5 #define min(x,y) ((x)<(y)?(x):(y))
8 #define max(x,y) ((x)>(y)?(x):(y))
16 #define SBid(i) ((i)>>logSB)
17 #define SBfirst(i) ((i) & (~(SB-1)))
18 #define SBlast(i) ((i) | (SB-1))
20 #define MBid(i) ((i)>>logMB)
21 #define MBfirst(i) ((i) & (~(MB-1)))
22 #define MBlast(i) ((i) | (MB-1))
24 pb getpat_preorder(pb *b)
29 pb getpat_postorder(pb *b)
38 w2 = (w1 << 1) + (b[1] >> (D-1));
43 pb getpat_inorder(pb *b)
47 w2 = (w1 << 1) + (b[1] >> (D-1));
52 pb getpat_dfuds_leaf(pb *b)
56 w2 = (w1 << 1) + (b[1] >> (D-1));
63 ///////////////////////////////////////////
64 // depth(bp *b, int s)
65 // returns the depth of s
66 // The root node has depth 1
67 ///////////////////////////////////////////
68 int depth(bp *b, int s)
72 d = 2 * darray_rank(b->da,s) - (s+1);
74 if (d != naive_depth(b,s)) {
78 //printf("depth(%d)=%d\n",s,d);
83 int fast_depth(bp *b, int s)
90 if ((s & (RRR-1)) != 0) {
91 //printf("fast_depth:warning s=%d\n",s);
95 //d = 2 * (da->rl[s>>logR] + da->rm[s>>logRR] + da->rs[s>>logRRR]) - s;
97 r = da->rl[s>>logR] + da->rm[s>>logRR];
98 j = (s>>logRRR) & (RR/RRR-1);
100 r += da->rs[((s>>logRR)<<(logRR-logRRR))+j-1];
108 int search_SB_r(bp *b, int i, int rel)
114 il = min((SBid(i) + 1) << logSB,n);
122 w = (x >> (D-ETW)) & ((1<<ETW)-1);
123 if (rel >= -ETW && rel <= ETW) {
124 r = fwdtbl[((rel+ETW)<<ETW)+w];
126 if (i+r >= n) return NOTFOUND;
131 rel -= 2*popCount[w]-r;
140 int search_MB_r(bp *b, int i, int td)
149 il = min((MBid(i) + 1) << logMB,n);
150 for ( ; i < il; i+=SB) {
154 d = fast_depth(b,i-1);
156 m = d + b->sm[SBid(i)] - SB;
157 M = d + b->sM[SBid(i)] - 1;
158 if (m <= td && td <= M) {
159 return search_SB_r(b,i,td-d);
162 if (i >= n) return NOTFOUND;
166 ///////////////////////////////////////////
167 // fwd_excess(bp *b,int s, int rel)
168 // find the leftmost value depth(s)+rel to the right of s (exclusive)
169 ///////////////////////////////////////////
170 int fwd_excess(bp *b,int s, int rel)
181 d = search_SB_r(b,i,rel);
182 if (d >= NOTFOUND) return d;
184 i = min((SBid(i) + 1) << logSB,n);
185 td = depth(b,s) + rel;
186 d = search_MB_r(b,i,td);
187 if (d >= NOTFOUND) return d;
189 if (i != SBfirst(i)) {
190 d = search_SB_r(b,i,rel);
191 if (d >= NOTFOUND) return d;
194 td = depth(b,s) + rel;
196 i = SBid(i+SB-1) << logSB;
198 if (i != MBfirst(i)) {
199 d = search_MB_r(b,i,td);
200 if (d >= NOTFOUND) return d;
211 if (m <= td && td <= M) break;
215 if (i == 0) return NOTFOUND;
221 if (!(m <= td && td <= M)) i++;
226 d = search_MB_r(b,i,td);
227 if (d >= NOTFOUND) return d;
230 printf("fwd_excess: ???\n");
236 int degree_SB_slow(bp *b, int i, int t, int rel, int *ans, int ith)
243 il = min((SBid(i) + 1) << logSB,n);
247 if (getbit(b->B,i)==OP) {
252 if (d < rel) { // reached the end
260 if (d == rel) { // found the same depth
273 int degree_SB(bp *b, int i, int t, int rel, int *ans, int ith)
282 il = min((SBid(i) + 1) << logSB,n);
292 w = (x >> (D-ETW)) & ((1<<ETW)-1);
294 if (j < ETW || il-i < ETW) {
295 r = max(ETW-j,ETW-(il-i));
298 v = minmaxtbl_v[0][w | w2];
302 for (r = 0; r < ETW; r++) {
320 r = childtbl2[rel-d+ETW][ith-1][w];
328 r = ETW-1-minmaxtbl_i[0][w | w2];
330 deg += degtbl2[((rel-d+ETW)<<ETW) + (w & (~w2))];
339 *ans = i + childtbl[((ith-1)<<ETW) + (w | w2)];
347 d += 2*popCount[w]-r;
358 int degree_MB(bp *b, int i, int t, int td, int *ans,int ith)
369 il = min((MBid(i) + 1) << logMB,n);
371 for ( ; i+SB-1 < il; i+=SB) {
375 d = fast_depth(b,i-1);
377 m = d + b->sm[SBid(i)] - SB;
379 r = degree_SB(b,i,n,td-d,°tmp,ith);
381 if (r == NOTFOUND) return NOTFOUND;
391 if (ith <= b->sd[SBid(i)]) break;
392 ith -= b->sd[SBid(i)];
394 deg += b->sd[SBid(i)];
401 d = fast_depth(b,i-1);
403 r = degree_SB(b,i,n,td-d,°tmp,ith);
405 if (r == NOTFOUND) return NOTFOUND;
415 if (i >= n) return END;
419 static int partition_range(int s,int t)
423 printf("partition [%d,%d] => ",s,t);
428 printf("[%d,%d] ",s,s+h-1);
438 printf("[%d,%d] ",s,s+h-1);
449 ///////////////////////////////////////////
450 // fast_degree(bp *b,int s, int t, int ith)
451 // returns the number of children of s, to the left of t
452 // returns the position of (ith)-th child of s if (ith > 0)
453 ///////////////////////////////////////////
454 int fast_degree(bp *b,int s, int t, int ith)
470 if (i != SBfirst(i)) {
471 d = degree_SB(b,i,n,0,°tmp,ith);
473 if (d == NOTFOUND) return -1;
474 if (d == FOUND) return degtmp;
477 if (d == END) return degtmp;
483 i = SBid(i+SB-1) << logSB;
485 if (i != MBfirst(i)) {
486 d = degree_MB(b,i,n,td,°tmp,ith);
488 if (d == NOTFOUND) return -1;
489 if (d == FOUND) return degtmp;
504 tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
507 sm += m_ofs; tm += m_ofs;
508 for (i=sm; i<=tm; i++) {
512 if (b->mm[i] == td) {
514 if (ith <= b->md[i]) break;
523 tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
526 sm += m_ofs; tm += m_ofs;
529 //partition_range(sm,tm);
531 //printf("partition [%d,%d] => ",sm,tm);
536 //printf("[%d,%d] ",sm,sm+h-1);
542 if (b->mm[j] == td) {
544 if (ith <= b->md[j]) {
555 if (sm+h > tm) break;
561 //printf("[%d,%d] ",sm,sm+h-1);
564 if (b->mm[j] >= td) {
565 if (b->mm[j] == td) {
566 if (ith > b->md[j]) {
577 if (b->mm[j] >= td) {
578 if (b->mm[j] == td) {
596 d = degree_MB(b,ss,n,td,°tmp,ith);
598 if (d == NOTFOUND) return -1;
599 if (d == FOUND) return degtmp;
602 if (d == END) return deg;
606 printf("degree: ???\n");
611 int search_SB_l(bp *b, int i, int rel)
624 w = x & ((1<<ETW)-1);
625 if (rel >= -ETW && rel <= ETW) {
626 r = bwdtbl[((rel+ETW)<<ETW)+w];
628 if (i-r < 0) return NOTFOUND;
633 rel += 2*popCount[w]-r;
640 if (i < 0) return NOTFOUND;
644 int search_MB_l(bp *b, int i, int td)
651 if (i % SB != SB-1) {
652 printf("search_MB_l:error!!! i=%d SB=%d\n",i,i%SB);
658 for ( ; i >= il; i-=SB) {
662 d = fast_depth(b,i-SB);
664 m = d + b->sm[SBid(i)] - SB;
665 M = d + b->sM[SBid(i)] - 1;
666 if (m <= td && td <= M) {
672 if (d == td) return i;
673 return search_SB_l(b,i,td-d);
679 ///////////////////////////////////////////
680 // bwd_excess(bp *b,int s, int rel)
681 // find the rightmost value depth(s)+rel to the left of s (exclusive)
682 ///////////////////////////////////////////
683 int bwd_excess(bp *b,int s, int rel)
693 d = search_SB_l(b,i,rel);
694 if (d >= NOTFOUND) return d;
698 td = depth(b,s) + rel;
700 d = search_MB_l(b,i,td);
701 if (d >= NOTFOUND) return d;
704 i = (s>>logMB) + m_ofs;
710 if (m <= td && td <= M) break;
715 if (td == 0) return -1;
716 else return NOTFOUND;
723 if (!(m <= td && td <= M)) i--;
726 i = ((i+1)<<logMB)-1;
728 d = search_MB_l(b,i,td);
729 if (d >= NOTFOUND) return d;
732 printf("bwd_excess: ???\n");
737 int rmq_SB(bp *b, int s, int t, int opt, int *dm)
746 lr = 0; if (opt & OPT_RIGHT) lr = 1;
747 op = opt & (OPT_RIGHT | OPT_MAX);
760 w = (x >> (D-ETW)) & ((1<<ETW)-1);
762 if (j < ETW || t-i < ETW-1) {
763 r = max(ETW-j,ETW-1-(t-i));
768 v = minmaxtbl_v[op][w & (~w2)];
769 if (d + v + lr > ds) {
771 is = i + minmaxtbl_i[op][w & (~w2)];
774 v = minmaxtbl_v[op][w | w2];
775 if (d + v < ds +lr) {
777 is = i + minmaxtbl_i[op][w | w2];
782 d += 2*popCount[w]-r;
790 if (getbit(b->B,i)==OP) {
799 if (!(op & OPT_MAX)) {
812 int rmq_MB(bp *b, int s, int t, int opt, int *dm)
818 lr = 0; if (opt & OPT_RIGHT) lr = 1;
821 for (i = s; i <= t; i++) {
823 d = depth(b,(i<<logSB)-1);
825 d = fast_depth(b,(i<<logSB)-1);
828 m = d + b->sM[i] - 1;
833 m = d + b->sm[i] - SB;
846 ///////////////////////////////////////////
847 // rmq(bp *b, int s, int t, int opt)
848 // returns the index of leftmost/rightmost minimum/maximum value
849 // in the range [s,t] (inclusive)
850 // returns the leftmost minimum if opt == 0
851 // the rightmost one if (opt & OPT_RIGHT)
852 // the maximum if (opt & OPT_MAX)
853 ///////////////////////////////////////////
854 int rmq(bp *b, int s, int t, int opt)
856 int ss, ts; // SB index of s and t
857 int sm, tm; // MB index of s and t
858 int ds; // current min value
859 int is; // current min index
860 int ys; // type of current min index
861 // 0: is is the index of min
862 // 1: is is the SB index
863 // 2: is is the MB index
869 lr = 0; if (opt & OPT_RIGHT) lr = 1;
872 if (s < 0 || t >= n || s > t) {
873 printf("rmq: error s=%d t=%d n=%d\n",s,t,n);
876 if (s == t) return s;
879 ////////////////////////////////////////////////////////////
880 // search the SB of s
881 ////////////////////////////////////////////////////////////
883 il = min(SBlast(s),t);
884 is = rmq_SB(b,s,il,opt,&dm);
885 if (il == t) { // scan reached the end of the range
888 ds = depth(b,s) + dm; ys = 0;
890 ////////////////////////////////////////////////////////////
891 // search the MB of s
892 ////////////////////////////////////////////////////////////
895 il = min(SBid(MBlast(s)),SBid(t)-1);
897 j = rmq_MB(b,ss,il,opt,&dm);
898 //if (dm < ds + lr) {
900 ds = dm; is = j; ys = 1;
903 ////////////////////////////////////////////////////////////
905 ////////////////////////////////////////////////////////////
913 sm += m_ofs; tm += m_ofs;
914 for (i=sm; i<=tm; i++) {
916 if (b->mM[i] + lr > ds) {
917 ds = b->mM[i]; is = i - m_ofs; ys = 2;
920 if (b->mm[i] < ds + lr) {
921 ds = b->mm[i]; is = i - m_ofs; ys = 2;
932 sm += m_ofs; tm += m_ofs;
935 if (b->mM[sm] + lr > ds) {
936 ds = b->mM[sm]; is = sm; ys = 2;
938 for (i=0; i<=h-2; i++) {
941 if (b->mM[j+1] + lr > ds) {
942 ds = b->mM[j+1]; is = j+1; ys = 2;
946 for (i=h-2; i>=0; i--) {
949 if (b->mM[j-1] + lr > ds) {
950 ds = b->mM[j-1]; is = j-1; ys = 2;
954 if (b->mM[tm] + lr > ds) {
955 ds = b->mM[tm]; is = tm; ys = 2;
960 if (b->mM[is+1] + lr > b->mM[is]) is++;
965 if (b->mm[sm] < ds + lr) {
966 ds = b->mm[sm]; is = sm; ys = 2;
968 for (i=0; i<=h-2; i++) {
971 if (b->mm[j+1] < ds + lr) {
972 ds = b->mm[j+1]; is = j+1; ys = 2;
976 for (i=h-2; i>=0; i--) {
979 if (b->mm[j-1] < ds + lr) {
980 ds = b->mm[j-1]; is = j-1; ys = 2;
984 if (b->mm[tm] < ds + lr) {
985 ds = b->mm[tm]; is = tm; ys = 2;
990 if (b->mm[is+1] < b->mm[is] + lr) is++;
999 ////////////////////////////////////////////////////////////
1000 // search the MB of t
1001 ////////////////////////////////////////////////////////////
1004 ts = max(SBid(MBfirst(t)),SBid(s)+1);
1005 il = SBid(SBfirst(t)-1);
1007 j = rmq_MB(b,ts,il,opt,&dm);
1008 //if (dm < ds + lr) {
1010 ds = dm; is = j; ys = 1;
1013 ////////////////////////////////////////////////////////////
1014 // search the SB of t
1015 ////////////////////////////////////////////////////////////
1018 j = rmq_SB(b,i,t,opt,&dm);
1019 d = depth(b,i) + dm;
1020 if (opt & OPT_MAX) {
1022 ds = d; is = j; ys = 0;
1026 ds = d; is = j; ys = 0;
1030 ////////////////////////////////////////////////////////////
1032 ////////////////////////////////////////////////////////////
1035 ss = SBid(is << logMB);
1036 il = SBid(MBlast(is << logMB));
1037 if (opt & OPT_MAX) {
1042 j = rmq_MB(b,ss,il,opt,&dm);
1043 ds = dm; is = j; ys = 1;
1048 il = SBlast(is << logSB);
1049 is = rmq_SB(b,ss,il,opt,&dm);
1050 //ds = depth(b,ss) + dm; ys = 0;