6 #define min(x,y) ((x)<(y)?(x):(y))
10 #define max(x,y) ((x)>(y)?(x):(y))
18 #define SBid(i) ((i)>>logSB)
19 #define SBfirst(i) ((i) & (~(SB-1)))
20 #define SBlast(i) ((i) | (SB-1))
22 #define MBid(i) ((i)>>logMB)
23 #define MBfirst(i) ((i) & (~(MB-1)))
24 #define MBlast(i) ((i) | (MB-1))
26 pb getpat_preorder(pb *b)
31 pb getpat_postorder(pb *b)
40 w2 = (w1 << 1) + (b[1] >> (D-1));
45 pb getpat_inorder(pb *b)
49 w2 = (w1 << 1) + (b[1] >> (D-1));
54 pb getpat_dfuds_leaf(pb *b)
58 w2 = (w1 << 1) + (b[1] >> (D-1));
65 ///////////////////////////////////////////
66 // depth(bp *b, int s)
67 // returns the depth of s
68 // The root node has depth 1
69 ///////////////////////////////////////////
70 int depth(bp *b, int s)
74 d = 2 * darray_rank(b->da,s) - (s+1);
76 if (d != naive_depth(b,s)) {
80 //printf("depth(%d)=%d\n",s,d);
85 int fast_depth(bp *b, int s)
92 if ((s & (RRR-1)) != 0) {
93 //printf("fast_depth:warning s=%d\n",s);
97 //d = 2 * (da->rl[s>>logR] + da->rm[s>>logRR] + da->rs[s>>logRRR]) - s;
99 r = da->rl[s>>logR] + da->rm[s>>logRR];
100 j = (s>>logRRR) & (RR/RRR-1);
102 r += da->rs[((s>>logRR)<<(logRR-logRRR))+j-1];
110 int search_SB_r(bp *b, int i, int rel)
116 il = min((SBid(i) + 1) << logSB,n);
124 w = (x >> (D-ETW)) & ((1<<ETW)-1);
125 if (rel >= -ETW && rel <= ETW) {
126 r = fwdtbl[((rel+ETW)<<ETW)+w];
128 if (i+r >= n) return NOTFOUND;
133 rel -= 2*popcount(w)-r;
142 int search_MB_r(bp *b, int i, int td)
151 il = min((MBid(i) + 1) << logMB,n);
152 for ( ; i < il; i+=SB) {
156 d = fast_depth(b,i-1);
158 m = d + b->sm[SBid(i)] - SB;
159 M = d + b->sM[SBid(i)] - 1;
160 if (m <= td && td <= M) {
161 return search_SB_r(b,i,td-d);
164 if (i >= n) return NOTFOUND;
168 ///////////////////////////////////////////
169 // fwd_excess(bp *b,int s, int rel)
170 // find the leftmost value depth(s)+rel to the right of s (exclusive)
171 ///////////////////////////////////////////
172 int fwd_excess(bp *b,int s, int rel)
183 d = search_SB_r(b,i,rel);
184 if (d >= NOTFOUND) return d;
186 i = min((SBid(i) + 1) << logSB, n);
187 td = depth(b,s) + rel;
188 d = search_MB_r(b,i,td);
189 if (d >= NOTFOUND) return d;
201 if (m <= td && td <= M) {
206 if (!(m <= td && td <= M)) i++;
211 return search_MB_r(b,i,td);
221 int degree_SB_slow(bp *b, int i, int t, int rel, int *ans, int ith)
228 il = min((SBid(i) + 1) << logSB,n);
232 if (getbit(b->B,i)==OP) {
237 if (d < rel) { // reached the end
245 if (d == rel) { // found the same depth
258 int degree_SB(bp *b, int i, int t, int rel, int *ans, int ith)
267 il = min((SBid(i) + 1) << logSB,n);
277 w = (x >> (D-ETW)) & ((1<<ETW)-1);
279 if (j < ETW || il-i < ETW) {
280 r = max(ETW-j,ETW-(il-i));
283 v = minmaxtbl_v[0][w | w2];
287 for (r = 0; r < ETW; r++) {
305 r = childtbl2[rel-d+ETW][ith-1][w];
313 r = ETW-1-minmaxtbl_i[0][w | w2];
315 deg += degtbl2[((rel-d+ETW)<<ETW) + (w & (~w2))];
324 *ans = i + childtbl[((ith-1)<<ETW) + (w | w2)];
332 d += 2*popcount(w)-r;
343 int degree_MB(bp *b, int i, int t, int td, int *ans,int ith)
354 il = min((MBid(i) + 1) << logMB,n);
356 for ( ; i+SB-1 < il; i+=SB) {
360 d = fast_depth(b,i-1);
362 m = d + b->sm[SBid(i)] - SB;
364 r = degree_SB(b,i,n,td-d,°tmp,ith);
366 if (r == NOTFOUND) return NOTFOUND;
376 if (ith <= b->sd[SBid(i)]) break;
377 ith -= b->sd[SBid(i)];
379 deg += b->sd[SBid(i)];
386 d = fast_depth(b,i-1);
388 r = degree_SB(b,i,n,td-d,°tmp,ith);
390 if (r == NOTFOUND) return NOTFOUND;
400 if (i >= n) return END;
404 static int partition_range(int s,int t)
408 printf("partition [%d,%d] => ",s,t);
413 printf("[%d,%d] ",s,s+h-1);
423 printf("[%d,%d] ",s,s+h-1);
434 ///////////////////////////////////////////
435 // fast_degree(bp *b,int s, int t, int ith)
436 // returns the number of children of s, to the left of t
437 // returns the position of (ith)-th child of s if (ith > 0)
438 ///////////////////////////////////////////
439 int fast_degree(bp *b,int s, int t, int ith)
455 if (i != SBfirst(i)) {
456 d = degree_SB(b,i,n,0,°tmp,ith);
458 if (d == NOTFOUND) return -1;
459 if (d == FOUND) return degtmp;
462 if (d == END) return degtmp;
468 i = SBid(i+SB-1) << logSB;
470 if (i != MBfirst(i)) {
471 d = degree_MB(b,i,n,td,°tmp,ith);
473 if (d == NOTFOUND) return -1;
474 if (d == FOUND) return degtmp;
489 tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
492 sm += m_ofs; tm += m_ofs;
493 for (i=sm; i<=tm; i++) {
497 if (b->mm[i] == td) {
499 if (ith <= b->md[i]) break;
508 tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
511 sm += m_ofs; tm += m_ofs;
514 //partition_range(sm,tm);
516 //printf("partition [%d,%d] => ",sm,tm);
521 //printf("[%d,%d] ",sm,sm+h-1);
527 if (b->mm[j] == td) {
529 if (ith <= b->md[j]) {
540 if (sm+h > tm) break;
546 //printf("[%d,%d] ",sm,sm+h-1);
549 if (b->mm[j] >= td) {
550 if (b->mm[j] == td) {
551 if (ith > b->md[j]) {
562 if (b->mm[j] >= td) {
563 if (b->mm[j] == td) {
581 d = degree_MB(b,ss,n,td,°tmp,ith);
583 if (d == NOTFOUND) return -1;
584 if (d == FOUND) return degtmp;
587 if (d == END) return deg;
591 printf("degree: ???\n");
596 int search_SB_l(bp *b, int i, int rel)
609 w = x & ((1<<ETW)-1);
610 if (rel >= -ETW && rel <= ETW) {
611 r = bwdtbl[((rel+ETW)<<ETW)+w];
613 if (i-r < 0) return NOTFOUND;
618 rel += 2*popcount(w)-r;
625 if (i < 0) return NOTFOUND;
629 int search_MB_l(bp *b, int i, int td)
636 if (i % SB != SB-1) {
637 printf("search_MB_l:error!!! i=%d SB=%d\n",i,i%SB);
643 for ( ; i >= il; i-=SB) {
647 d = fast_depth(b,i-SB);
649 m = d + b->sm[SBid(i)] - SB;
650 M = d + b->sM[SBid(i)] - 1;
651 if (m <= td && td <= M) {
657 if (d == td) return i;
658 return search_SB_l(b,i,td-d);
664 ///////////////////////////////////////////
665 // bwd_excess(bp *b,int s, int rel)
666 // find the rightmost value depth(s)+rel to the left of s (exclusive)
667 ///////////////////////////////////////////
668 int bwd_excess(bp *b,int s, int rel)
678 d = search_SB_l(b,i,rel);
679 if (d >= NOTFOUND) return d;
683 td = depth(b,s) + rel;
685 d = search_MB_l(b,i,td);
686 if (d >= NOTFOUND) return d;
689 i = (s>>logMB) + m_ofs;
695 if (m <= td && td <= M) break;
700 if (td == 0) return -1;
701 else return NOTFOUND;
708 if (!(m <= td && td <= M)) i--;
711 i = ((i+1)<<logMB)-1;
713 d = search_MB_l(b,i,td);
714 if (d >= NOTFOUND) return d;
717 printf("bwd_excess: ???\n");
722 int rmq_SB(bp *b, int s, int t, int opt, int *dm)
731 lr = 0; if (opt & OPT_RIGHT) lr = 1;
732 op = opt & (OPT_RIGHT | OPT_MAX);
745 w = (x >> (D-ETW)) & ((1<<ETW)-1);
747 if (j < ETW || t-i < ETW-1) {
748 r = max(ETW-j,ETW-1-(t-i));
753 v = minmaxtbl_v[op][w & (~w2)];
754 if (d + v + lr > ds) {
756 is = i + minmaxtbl_i[op][w & (~w2)];
759 v = minmaxtbl_v[op][w | w2];
760 if (d + v < ds +lr) {
762 is = i + minmaxtbl_i[op][w | w2];
767 d += 2*popcount(w)-r;
775 if (getbit(b->B,i)==OP) {
784 if (!(op & OPT_MAX)) {
797 int rmq_MB(bp *b, int s, int t, int opt, int *dm)
803 lr = 0; if (opt & OPT_RIGHT) lr = 1;
806 for (i = s; i <= t; i++) {
808 d = depth(b,(i<<logSB)-1);
810 d = fast_depth(b,(i<<logSB)-1);
813 m = d + b->sM[i] - 1;
818 m = d + b->sm[i] - SB;
831 ///////////////////////////////////////////
832 // rmq(bp *b, int s, int t, int opt)
833 // returns the index of leftmost/rightmost minimum/maximum value
834 // in the range [s,t] (inclusive)
835 // returns the leftmost minimum if opt == 0
836 // the rightmost one if (opt & OPT_RIGHT)
837 // the maximum if (opt & OPT_MAX)
838 ///////////////////////////////////////////
839 int rmq(bp *b, int s, int t, int opt)
841 int ss, ts; // SB index of s and t
842 int sm, tm; // MB index of s and t
843 int ds; // current min value
844 int is; // current min index
845 int ys; // type of current min index
846 // 0: is is the index of min
847 // 1: is is the SB index
848 // 2: is is the MB index
854 lr = 0; if (opt & OPT_RIGHT) lr = 1;
857 if (s < 0 || t >= n || s > t) {
858 printf("rmq: error s=%d t=%d n=%d\n",s,t,n);
861 if (s == t) return s;
864 ////////////////////////////////////////////////////////////
865 // search the SB of s
866 ////////////////////////////////////////////////////////////
868 il = min(SBlast(s),t);
869 is = rmq_SB(b,s,il,opt,&dm);
870 if (il == t) { // scan reached the end of the range
873 ds = depth(b,s) + dm; ys = 0;
875 ////////////////////////////////////////////////////////////
876 // search the MB of s
877 ////////////////////////////////////////////////////////////
880 il = min(SBid(MBlast(s)),SBid(t)-1);
882 j = rmq_MB(b,ss,il,opt,&dm);
883 //if (dm < ds + lr) {
885 ds = dm; is = j; ys = 1;
888 ////////////////////////////////////////////////////////////
890 ////////////////////////////////////////////////////////////
898 sm += m_ofs; tm += m_ofs;
899 for (i=sm; i<=tm; i++) {
901 if (b->mM[i] + lr > ds) {
902 ds = b->mM[i]; is = i - m_ofs; ys = 2;
905 if (b->mm[i] < ds + lr) {
906 ds = b->mm[i]; is = i - m_ofs; ys = 2;
917 sm += m_ofs; tm += m_ofs;
920 if (b->mM[sm] + lr > ds) {
921 ds = b->mM[sm]; is = sm; ys = 2;
923 for (i=0; i<=h-2; i++) {
926 if (b->mM[j+1] + lr > ds) {
927 ds = b->mM[j+1]; is = j+1; ys = 2;
931 for (i=h-2; i>=0; i--) {
934 if (b->mM[j-1] + lr > ds) {
935 ds = b->mM[j-1]; is = j-1; ys = 2;
939 if (b->mM[tm] + lr > ds) {
940 ds = b->mM[tm]; is = tm; ys = 2;
945 if (b->mM[is+1] + lr > b->mM[is]) is++;
950 if (b->mm[sm] < ds + lr) {
951 ds = b->mm[sm]; is = sm; ys = 2;
953 for (i=0; i<=h-2; i++) {
956 if (b->mm[j+1] < ds + lr) {
957 ds = b->mm[j+1]; is = j+1; ys = 2;
961 for (i=h-2; i>=0; i--) {
964 if (b->mm[j-1] < ds + lr) {
965 ds = b->mm[j-1]; is = j-1; ys = 2;
969 if (b->mm[tm] < ds + lr) {
970 ds = b->mm[tm]; is = tm; ys = 2;
975 if (b->mm[is+1] < b->mm[is] + lr) is++;
984 ////////////////////////////////////////////////////////////
985 // search the MB of t
986 ////////////////////////////////////////////////////////////
989 ts = max(SBid(MBfirst(t)),SBid(s)+1);
990 il = SBid(SBfirst(t)-1);
992 j = rmq_MB(b,ts,il,opt,&dm);
993 //if (dm < ds + lr) {
995 ds = dm; is = j; ys = 1;
998 ////////////////////////////////////////////////////////////
999 // search the SB of t
1000 ////////////////////////////////////////////////////////////
1003 j = rmq_SB(b,i,t,opt,&dm);
1004 d = depth(b,i) + dm;
1005 if (opt & OPT_MAX) {
1007 ds = d; is = j; ys = 0;
1011 ds = d; is = j; ys = 0;
1015 ////////////////////////////////////////////////////////////
1017 ////////////////////////////////////////////////////////////
1020 ss = SBid(is << logMB);
1021 il = SBid(MBlast(is << logMB));
1022 if (opt & OPT_MAX) {
1027 j = rmq_MB(b,ss,il,opt,&dm);
1028 ds = dm; is = j; ys = 1;
1033 il = SBlast(is << logSB);
1034 is = rmq_SB(b,ss,il,opt,&dm);
1035 //ds = depth(b,ss) + dm; ys = 0;