10 #define SBid(i) ((i)>>logSB)
11 #define SBfirst(i) ((i) & (~(SB-1)))
12 #define SBlast(i) ((i) | (SB-1))
14 #define MBid(i) ((i)>>logMB)
15 #define MBfirst(i) ((i) & (~(MB-1)))
16 #define MBlast(i) ((i) | (MB-1))
18 ({ __typeof__ (a) _a = (a); \
19 __typeof__ (b) _b = (b); \
24 ({ __typeof__ (a) _a = (a); \
25 __typeof__ (b) _b = (b); \
26 _a <= _b ? _a : _b; })
29 pb getpat_preorder(pb *b)
34 pb getpat_postorder(pb *b)
43 w2 = (w1 << 1) + (b[1] >> (D-1));
48 pb getpat_inorder(pb *b)
52 w2 = (w1 << 1) + (b[1] >> (D-1));
57 pb getpat_dfuds_leaf(pb *b)
61 w2 = (w1 << 1) + (b[1] >> (D-1));
68 ///////////////////////////////////////////
69 // bp_depth(bp *b, int s)
70 // returns the depth of s
71 // The root node has depth 1
72 ///////////////////////////////////////////
73 int bp_depth(bp *b, int s)
77 d = 2 * bp_darray_rank(b->da,s) - (s+1);
81 int fast_depth(bp *b, int s)
88 if ((s & (RRR-1)) != 0) {
89 //printf("fast_depth:warning s=%d\n",s);
93 //d = 2 * (da->rl[s>>logR] + da->rm[s>>logRR] + da->rs[s>>logRRR]) - s;
95 r = da->rl[s>>logR] + da->rm[s>>logRR];
96 j = (s>>logRRR) & (RR/RRR-1);
98 r += da->rs[((s>>logRR)<<(logRR-logRRR))+j-1];
106 int search_SB_r(bp *b, int i, int rel)
112 il = min((SBid(i) + 1) << logSB,n);
120 w = (x >> (D-ETW)) & ((1<<ETW)-1);
121 if (rel >= -ETW && rel <= ETW) {
122 r = fwdtbl[((rel+ETW)<<ETW)+w];
124 if (i+r >= n) return NOTFOUND;
129 rel -= 2*popcount(w)-r;
138 int search_MB_r(bp *b, int i, int td)
147 il = min((MBid(i) + 1) << logMB, n);
148 for ( ; i < il; i+=SB) {
152 d = fast_depth(b,i-1);
154 m = d + b->sm[SBid(i)] - SB;
155 M = d + b->sM[SBid(i)] - 1;
156 if (m <= td && td <= M) {
157 return search_SB_r(b,i,td-d);
160 if (i >= n) return NOTFOUND;
164 ///////////////////////////////////////////
165 // bp_fwd_excess(bp *b,int s, int rel)
166 // find the leftmost value depth(s)+rel to the right of s (exclusive)
167 ///////////////////////////////////////////
168 int bp_fwd_excess(bp *b,int s, int rel)
180 d = search_SB_r(b,i,rel);
181 if (d >= NOTFOUND) return d;
183 i = min((SBid(i) + 1) << logSB, n);
184 td = bp_depth(b,s) + rel;
185 d = search_MB_r(b,i,td);
186 if (d >= NOTFOUND) return d;
198 if (m <= td && td <= M) {
203 if (!(m <= td && td <= M)) i++;
208 return search_MB_r(b,i,td);
218 int degree_SB_slow(bp *b, int i, int t, int rel, int *ans, int ith)
225 il = min((SBid(i) + 1) << logSB,n);
229 if (bp_getbit(b->B,i)==OP) {
234 if (d < rel) { // reached the end
242 if (d == rel) { // found the same depth
255 int degree_SB(bp *b, int i, int t, int rel, int *ans, int ith)
264 il = min((SBid(i) + 1) << logSB,n);
274 w = (x >> (D-ETW)) & ((1<<ETW)-1);
276 if (j < ETW || il-i < ETW) {
277 r = max(ETW-j,ETW-(il-i));
280 v = minmaxtbl_v[0][w | w2];
284 for (r = 0; r < ETW; r++) {
302 r = childtbl2[rel-d+ETW][ith-1][w];
310 r = ETW-1-minmaxtbl_i[0][w | w2];
312 deg += degtbl2[((rel-d+ETW)<<ETW) + (w & (~w2))];
321 *ans = i + childtbl[((ith-1)<<ETW) + (w | w2)];
329 d += 2*popcount(w)-r;
340 int degree_MB(bp *b, int i, int t, int td, int *ans,int ith)
351 il = min((MBid(i) + 1) << logMB,n);
353 for ( ; i+SB-1 < il; i+=SB) {
357 d = fast_depth(b,i-1);
359 m = d + b->sm[SBid(i)] - SB;
361 r = degree_SB(b,i,n,td-d,°tmp,ith);
363 if (r == NOTFOUND) return NOTFOUND;
373 if (ith <= b->sd[SBid(i)]) break;
374 ith -= b->sd[SBid(i)];
376 deg += b->sd[SBid(i)];
383 d = fast_depth(b,i-1);
385 r = degree_SB(b,i,n,td-d,°tmp,ith);
387 if (r == NOTFOUND) return NOTFOUND;
397 if (i >= n) return END;
401 static void partition_range(int s,int t)
405 printf("partition [%d,%d] => ",s,t);
410 printf("[%d,%d] ",s,s+h-1);
420 printf("[%d,%d] ",s,s+h-1);
431 ///////////////////////////////////////////
432 // bp_fast_degree(bp *b,int s, int t, int ith)
433 // returns the number of children of s, to the left of t
434 // returns the position of (ith)-th child of s if (ith > 0)
435 ///////////////////////////////////////////
436 int bp_fast_degree(bp *b,int s, int t, int ith)
452 if (i != SBfirst(i)) {
453 d = degree_SB(b,i,n,0,°tmp,ith);
455 if (d == NOTFOUND) return -1;
456 if (d == FOUND) return degtmp;
459 if (d == END) return degtmp;
465 i = SBid(i+SB-1) << logSB;
467 if (i != MBfirst(i)) {
468 d = degree_MB(b,i,n,td,°tmp,ith);
470 if (d == NOTFOUND) return -1;
471 if (d == FOUND) return degtmp;
486 tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
489 sm += m_ofs; tm += m_ofs;
490 for (i=sm; i<=tm; i++) {
494 if (b->mm[i] == td) {
496 if (ith <= b->md[i]) break;
505 tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
508 sm += m_ofs; tm += m_ofs;
511 //partition_range(sm,tm);
513 //printf("partition [%d,%d] => ",sm,tm);
518 //printf("[%d,%d] ",sm,sm+h-1);
524 if (b->mm[j] == td) {
526 if (ith <= b->md[j]) {
537 if (sm+h > tm) break;
543 //printf("[%d,%d] ",sm,sm+h-1);
546 if (b->mm[j] >= td) {
547 if (b->mm[j] == td) {
548 if (ith > b->md[j]) {
559 if (b->mm[j] >= td) {
560 if (b->mm[j] == td) {
578 d = degree_MB(b,ss,n,td,°tmp,ith);
580 if (d == NOTFOUND) return -1;
581 if (d == FOUND) return degtmp;
584 if (d == END) return deg;
588 printf("degree: ???\n");
593 int search_SB_l(bp *b, int i, int rel)
606 w = x & ((1<<ETW)-1);
607 if (rel >= -ETW && rel <= ETW) {
608 r = bwdtbl[((rel+ETW)<<ETW)+w];
610 if (i-r < 0) return NOTFOUND;
615 rel += 2*popcount(w)-r;
622 if (i < 0) return NOTFOUND;
626 int search_MB_l(bp *b, int i, int td)
633 if (i % SB != SB-1) {
634 printf("search_MB_l:error!!! i=%d SB=%d\n",i,i%SB);
640 for ( ; i >= il; i-=SB) {
642 d = bp_depth(b,i-SB);
644 d = fast_depth(b,i-SB);
646 m = d + b->sm[SBid(i)] - SB;
647 M = d + b->sM[SBid(i)] - 1;
648 if (m <= td && td <= M) {
654 if (d == td) return i;
655 return search_SB_l(b,i,td-d);
661 ///////////////////////////////////////////
662 // bwd_excess(bp *b,int s, int rel)
663 // find the rightmost value depth(s)+rel to the left of s (exclusive)
664 ///////////////////////////////////////////
665 int bp_bwd_excess(bp *b,int s, int rel)
675 d = search_SB_l(b,i,rel);
676 if (d >= NOTFOUND) return d;
680 td = bp_depth(b,s) + rel;
682 d = search_MB_l(b,i,td);
683 if (d >= NOTFOUND) return d;
686 i = (s>>logMB) + m_ofs;
692 if (m <= td && td <= M) break;
697 if (td == 0) return -1;
698 else return NOTFOUND;
705 if (!(m <= td && td <= M)) i--;
708 i = ((i+1)<<logMB)-1;
710 d = search_MB_l(b,i,td);
711 if (d >= NOTFOUND) return d;
714 printf("bwd_excess: ???\n");
719 int rmq_SB(bp *b, int s, int t, int opt, int *dm)
728 lr = 0; if (opt & OPT_RIGHT) lr = 1;
729 op = opt & (OPT_RIGHT | OPT_MAX);
742 w = (x >> (D-ETW)) & ((1<<ETW)-1);
744 if (j < ETW || t-i < ETW-1) {
745 r = max(ETW-j,ETW-1-(t-i));
750 v = minmaxtbl_v[op][w & (~w2)];
751 if (d + v + lr > ds) {
753 is = i + minmaxtbl_i[op][w & (~w2)];
756 v = minmaxtbl_v[op][w | w2];
757 if (d + v < ds +lr) {
759 is = i + minmaxtbl_i[op][w | w2];
764 d += 2*popcount(w)-r;
772 if (bp_getbit(b->B,i)==OP) {
781 if (!(op & OPT_MAX)) {
794 int rmq_MB(bp *b, int s, int t, int opt, int *dm)
800 lr = 0; if (opt & OPT_RIGHT) lr = 1;
803 for (i = s; i <= t; i++) {
805 d = bp_depth(b,(i<<logSB)-1);
807 d = fast_depth(b,(i<<logSB)-1);
810 m = d + b->sM[i] - 1;
815 m = d + b->sm[i] - SB;
828 ///////////////////////////////////////////
829 // bp_rmq(bp *b, int s, int t, int opt)
830 // returns the index of leftmost/rightmost minimum/maximum value
831 // in the range [s,t] (inclusive)
832 // returns the leftmost minimum if opt == 0
833 // the rightmost one if (opt & OPT_RIGHT)
834 // the maximum if (opt & OPT_MAX)
835 ///////////////////////////////////////////
836 int bp_rmq(bp *b, int s, int t, int opt)
838 int ss, ts; // SB index of s and t
839 int sm, tm; // MB index of s and t
840 int ds; // current min value
841 int is; // current min index
842 int ys; // type of current min index
843 // 0: is is the index of min
844 // 1: is is the SB index
845 // 2: is is the MB index
851 lr = 0; if (opt & OPT_RIGHT) lr = 1;
854 if (s < 0 || t >= n || s > t) {
855 printf("rmq: error s=%d t=%d n=%d\n",s,t,n);
858 if (s == t) return s;
861 ////////////////////////////////////////////////////////////
862 // search the SB of s
863 ////////////////////////////////////////////////////////////
865 il = min(SBlast(s),t);
866 is = rmq_SB(b,s,il,opt,&dm);
867 if (il == t) { // scan reached the end of the range
870 ds = bp_depth(b,s) + dm; ys = 0;
872 ////////////////////////////////////////////////////////////
873 // search the MB of s
874 ////////////////////////////////////////////////////////////
877 il = min(SBid(MBlast(s)),SBid(t)-1);
879 j = rmq_MB(b,ss,il,opt,&dm);
880 //if (dm < ds + lr) {
882 ds = dm; is = j; ys = 1;
885 ////////////////////////////////////////////////////////////
887 ////////////////////////////////////////////////////////////
895 sm += m_ofs; tm += m_ofs;
896 for (i=sm; i<=tm; i++) {
898 if (b->mM[i] + lr > ds) {
899 ds = b->mM[i]; is = i - m_ofs; ys = 2;
902 if (b->mm[i] < ds + lr) {
903 ds = b->mm[i]; is = i - m_ofs; ys = 2;
914 sm += m_ofs; tm += m_ofs;
917 if (b->mM[sm] + lr > ds) {
918 ds = b->mM[sm]; is = sm; ys = 2;
920 for (i=0; i<=h-2; i++) {
923 if (b->mM[j+1] + lr > ds) {
924 ds = b->mM[j+1]; is = j+1; ys = 2;
928 for (i=h-2; i>=0; i--) {
931 if (b->mM[j-1] + lr > ds) {
932 ds = b->mM[j-1]; is = j-1; ys = 2;
936 if (b->mM[tm] + lr > ds) {
937 ds = b->mM[tm]; is = tm; ys = 2;
942 if (b->mM[is+1] + lr > b->mM[is]) is++;
947 if (b->mm[sm] < ds + lr) {
948 ds = b->mm[sm]; is = sm; ys = 2;
950 for (i=0; i<=h-2; i++) {
953 if (b->mm[j+1] < ds + lr) {
954 ds = b->mm[j+1]; is = j+1; ys = 2;
958 for (i=h-2; i>=0; i--) {
961 if (b->mm[j-1] < ds + lr) {
962 ds = b->mm[j-1]; is = j-1; ys = 2;
966 if (b->mm[tm] < ds + lr) {
967 ds = b->mm[tm]; is = tm; ys = 2;
972 if (b->mm[is+1] < b->mm[is] + lr) is++;
981 ////////////////////////////////////////////////////////////
982 // search the MB of t
983 ////////////////////////////////////////////////////////////
986 ts = max(SBid(MBfirst(t)),SBid(s)+1);
987 il = SBid(SBfirst(t)-1);
989 j = rmq_MB(b,ts,il,opt,&dm);
990 //if (dm < ds + lr) {
992 ds = dm; is = j; ys = 1;
995 ////////////////////////////////////////////////////////////
996 // search the SB of t
997 ////////////////////////////////////////////////////////////
1000 j = rmq_SB(b,i,t,opt,&dm);
1001 d = bp_depth(b,i) + dm;
1002 if (opt & OPT_MAX) {
1004 ds = d; is = j; ys = 0;
1008 ds = d; is = j; ys = 0;
1012 ////////////////////////////////////////////////////////////
1014 ////////////////////////////////////////////////////////////
1017 ss = SBid(is << logMB);
1018 il = SBid(MBlast(is << logMB));
1019 if (opt & OPT_MAX) {
1024 j = rmq_MB(b,ss,il,opt,&dm);
1025 ds = dm; is = j; ys = 1;
1030 il = SBlast(is << logSB);
1031 is = rmq_SB(b,ss,il,opt,&dm);
1032 //ds = bp_depth(b,ss) + dm; ys = 0;