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))
31 pb getpat_preorder(pb *b)
36 pb getpat_postorder(pb *b)
45 w2 = (w1 << 1) + (b[1] >> (D-1));
50 pb getpat_inorder(pb *b)
54 w2 = (w1 << 1) + (b[1] >> (D-1));
59 pb getpat_dfuds_leaf(pb *b)
63 w2 = (w1 << 1) + (b[1] >> (D-1));
70 ///////////////////////////////////////////
71 // depth(bp *b, int s)
72 // returns the depth of s
73 // The root node has depth 1
74 ///////////////////////////////////////////
75 int depth(bp *b, int s)
79 d = 2 * darray_rank(b->da,s) - (s+1);
81 if (d != naive_depth(b,s)) {
85 //printf("depth(%d)=%d\n",s,d);
90 int fast_depth(bp *b, int s)
97 if ((s & (RRR-1)) != 0) {
98 //printf("fast_depth:warning s=%d\n",s);
102 //d = 2 * (da->rl[s>>logR] + da->rm[s>>logRR] + da->rs[s>>logRRR]) - s;
104 r = da->rl[s>>logR] + da->rm[s>>logRR];
105 j = (s>>logRRR) & (RR/RRR-1);
107 r += da->rs[((s>>logRR)<<(logRR-logRRR))+j-1];
115 int search_SB_r(bp *b, int i, int rel)
121 il = min((SBid(i) + 1) << logSB,n);
129 w = (x >> (D-ETW)) & ((1<<ETW)-1);
130 if (rel >= -ETW && rel <= ETW) {
131 r = fwdtbl[((rel+ETW)<<ETW)+w];
133 if (i+r >= n) return NOTFOUND;
138 rel -= 2*popCount[w]-r;
147 int search_MB_r(bp *b, int i, int td)
156 il = min((MBid(i) + 1) << logMB,n);
157 for ( ; i < il; i+=SB) {
161 d = fast_depth(b,i-1);
163 m = d + b->sm[SBid(i)] - SB;
164 M = d + b->sM[SBid(i)] - 1;
165 if (m <= td && td <= M) {
166 return search_SB_r(b,i,td-d);
169 if (i >= n) return NOTFOUND;
173 ///////////////////////////////////////////
174 // fwd_excess(bp *b,int s, int rel)
175 // find the leftmost value depth(s)+rel to the right of s (exclusive)
176 ///////////////////////////////////////////
177 int fwd_excess(bp *b,int s, int rel)
188 d = search_SB_r(b,i,rel);
189 if (d >= NOTFOUND) return d;
191 i = min((SBid(i) + 1) << logSB,n);
192 td = depth(b,s) + rel;
193 d = search_MB_r(b,i,td);
194 if (d >= NOTFOUND) return d;
196 if (i != SBfirst(i)) {
197 d = search_SB_r(b,i,rel);
198 if (d >= NOTFOUND) return d;
201 td = depth(b,s) + rel;
203 i = SBid(i+SB-1) << logSB;
205 if (i != MBfirst(i)) {
206 d = search_MB_r(b,i,td);
207 if (d >= NOTFOUND) return d;
218 if (m <= td && td <= M) break;
222 if (i == 0) return NOTFOUND;
228 if (!(m <= td && td <= M)) i++;
233 d = search_MB_r(b,i,td);
234 if (d >= NOTFOUND) return d;
237 printf("fwd_excess: ???\n");
243 int degree_SB_slow(bp *b, int i, int t, int rel, int *ans, int ith)
250 il = min((SBid(i) + 1) << logSB,n);
254 if (getbit(b->B,i)==OP) {
259 if (d < rel) { // reached the end
267 if (d == rel) { // found the same depth
280 int degree_SB(bp *b, int i, int t, int rel, int *ans, int ith)
289 il = min((SBid(i) + 1) << logSB,n);
299 w = (x >> (D-ETW)) & ((1<<ETW)-1);
301 if (j < ETW || il-i < ETW) {
302 r = max(ETW-j,ETW-(il-i));
305 v = minmaxtbl_v[0][w | w2];
309 for (r = 0; r < ETW; r++) {
327 r = childtbl2[rel-d+ETW][ith-1][w];
335 r = ETW-1-minmaxtbl_i[0][w | w2];
337 deg += degtbl2[((rel-d+ETW)<<ETW) + (w & (~w2))];
346 *ans = i + childtbl[((ith-1)<<ETW) + (w | w2)];
354 d += 2*popCount[w]-r;
365 int degree_MB(bp *b, int i, int t, int td, int *ans,int ith)
376 il = min((MBid(i) + 1) << logMB,n);
378 for ( ; i+SB-1 < il; i+=SB) {
382 d = fast_depth(b,i-1);
384 m = d + b->sm[SBid(i)] - SB;
386 r = degree_SB(b,i,n,td-d,°tmp,ith);
388 if (r == NOTFOUND) return NOTFOUND;
398 if (ith <= b->sd[SBid(i)]) break;
399 ith -= b->sd[SBid(i)];
401 deg += b->sd[SBid(i)];
408 d = fast_depth(b,i-1);
410 r = degree_SB(b,i,n,td-d,°tmp,ith);
412 if (r == NOTFOUND) return NOTFOUND;
422 if (i >= n) return END;
426 static int partition_range(int s,int t)
430 printf("partition [%d,%d] => ",s,t);
435 printf("[%d,%d] ",s,s+h-1);
445 printf("[%d,%d] ",s,s+h-1);
456 ///////////////////////////////////////////
457 // fast_degree(bp *b,int s, int t, int ith)
458 // returns the number of children of s, to the left of t
459 // returns the position of (ith)-th child of s if (ith > 0)
460 ///////////////////////////////////////////
461 int fast_degree(bp *b,int s, int t, int ith)
477 if (i != SBfirst(i)) {
478 d = degree_SB(b,i,n,0,°tmp,ith);
480 if (d == NOTFOUND) return -1;
481 if (d == FOUND) return degtmp;
484 if (d == END) return degtmp;
490 i = SBid(i+SB-1) << logSB;
492 if (i != MBfirst(i)) {
493 d = degree_MB(b,i,n,td,°tmp,ith);
495 if (d == NOTFOUND) return -1;
496 if (d == FOUND) return degtmp;
511 tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
514 sm += m_ofs; tm += m_ofs;
515 for (i=sm; i<=tm; i++) {
519 if (b->mm[i] == td) {
521 if (ith <= b->md[i]) break;
530 tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
533 sm += m_ofs; tm += m_ofs;
536 //partition_range(sm,tm);
538 //printf("partition [%d,%d] => ",sm,tm);
543 //printf("[%d,%d] ",sm,sm+h-1);
549 if (b->mm[j] == td) {
551 if (ith <= b->md[j]) {
562 if (sm+h > tm) break;
568 //printf("[%d,%d] ",sm,sm+h-1);
571 if (b->mm[j] >= td) {
572 if (b->mm[j] == td) {
573 if (ith > b->md[j]) {
584 if (b->mm[j] >= td) {
585 if (b->mm[j] == td) {
603 d = degree_MB(b,ss,n,td,°tmp,ith);
605 if (d == NOTFOUND) return -1;
606 if (d == FOUND) return degtmp;
609 if (d == END) return deg;
613 printf("degree: ???\n");
618 int search_SB_l(bp *b, int i, int rel)
631 w = x & ((1<<ETW)-1);
632 if (rel >= -ETW && rel <= ETW) {
633 r = bwdtbl[((rel+ETW)<<ETW)+w];
635 if (i-r < 0) return NOTFOUND;
640 rel += 2*popCount[w]-r;
647 if (i < 0) return NOTFOUND;
651 int search_MB_l(bp *b, int i, int td)
658 if (i % SB != SB-1) {
659 printf("search_MB_l:error!!! i=%d SB=%d\n",i,i%SB);
665 for ( ; i >= il; i-=SB) {
669 d = fast_depth(b,i-SB);
671 m = d + b->sm[SBid(i)] - SB;
672 M = d + b->sM[SBid(i)] - 1;
673 if (m <= td && td <= M) {
679 if (d == td) return i;
680 return search_SB_l(b,i,td-d);
686 ///////////////////////////////////////////
687 // bwd_excess(bp *b,int s, int rel)
688 // find the rightmost value depth(s)+rel to the left of s (exclusive)
689 ///////////////////////////////////////////
690 int bwd_excess(bp *b,int s, int rel)
700 d = search_SB_l(b,i,rel);
701 if (d >= NOTFOUND) return d;
705 td = depth(b,s) + rel;
707 d = search_MB_l(b,i,td);
708 if (d >= NOTFOUND) return d;
711 i = (s>>logMB) + m_ofs;
717 if (m <= td && td <= M) break;
722 if (td == 0) return -1;
723 else return NOTFOUND;
730 if (!(m <= td && td <= M)) i--;
733 i = ((i+1)<<logMB)-1;
735 d = search_MB_l(b,i,td);
736 if (d >= NOTFOUND) return d;
739 printf("bwd_excess: ???\n");
744 int rmq_SB(bp *b, int s, int t, int opt, int *dm)
753 lr = 0; if (opt & OPT_RIGHT) lr = 1;
754 op = opt & (OPT_RIGHT | OPT_MAX);
767 w = (x >> (D-ETW)) & ((1<<ETW)-1);
769 if (j < ETW || t-i < ETW-1) {
770 r = max(ETW-j,ETW-1-(t-i));
775 v = minmaxtbl_v[op][w & (~w2)];
776 if (d + v + lr > ds) {
778 is = i + minmaxtbl_i[op][w & (~w2)];
781 v = minmaxtbl_v[op][w | w2];
782 if (d + v < ds +lr) {
784 is = i + minmaxtbl_i[op][w | w2];
789 d += 2*popCount[w]-r;
797 if (getbit(b->B,i)==OP) {
806 if (!(op & OPT_MAX)) {
819 int rmq_MB(bp *b, int s, int t, int opt, int *dm)
825 lr = 0; if (opt & OPT_RIGHT) lr = 1;
828 for (i = s; i <= t; i++) {
830 d = depth(b,(i<<logSB)-1);
832 d = fast_depth(b,(i<<logSB)-1);
835 m = d + b->sM[i] - 1;
840 m = d + b->sm[i] - SB;
853 ///////////////////////////////////////////
854 // rmq(bp *b, int s, int t, int opt)
855 // returns the index of leftmost/rightmost minimum/maximum value
856 // in the range [s,t] (inclusive)
857 // returns the leftmost minimum if opt == 0
858 // the rightmost one if (opt & OPT_RIGHT)
859 // the maximum if (opt & OPT_MAX)
860 ///////////////////////////////////////////
861 int rmq(bp *b, int s, int t, int opt)
863 int ss, ts; // SB index of s and t
864 int sm, tm; // MB index of s and t
865 int ds; // current min value
866 int is; // current min index
867 int ys; // type of current min index
868 // 0: is is the index of min
869 // 1: is is the SB index
870 // 2: is is the MB index
876 lr = 0; if (opt & OPT_RIGHT) lr = 1;
879 if (s < 0 || t >= n || s > t) {
880 printf("rmq: error s=%d t=%d n=%d\n",s,t,n);
883 if (s == t) return s;
886 ////////////////////////////////////////////////////////////
887 // search the SB of s
888 ////////////////////////////////////////////////////////////
890 il = min(SBlast(s),t);
891 is = rmq_SB(b,s,il,opt,&dm);
892 if (il == t) { // scan reached the end of the range
895 ds = depth(b,s) + dm; ys = 0;
897 ////////////////////////////////////////////////////////////
898 // search the MB of s
899 ////////////////////////////////////////////////////////////
902 il = min(SBid(MBlast(s)),SBid(t)-1);
904 j = rmq_MB(b,ss,il,opt,&dm);
905 //if (dm < ds + lr) {
907 ds = dm; is = j; ys = 1;
910 ////////////////////////////////////////////////////////////
912 ////////////////////////////////////////////////////////////
920 sm += m_ofs; tm += m_ofs;
921 for (i=sm; i<=tm; i++) {
923 if (b->mM[i] + lr > ds) {
924 ds = b->mM[i]; is = i - m_ofs; ys = 2;
927 if (b->mm[i] < ds + lr) {
928 ds = b->mm[i]; is = i - m_ofs; ys = 2;
939 sm += m_ofs; tm += m_ofs;
942 if (b->mM[sm] + lr > ds) {
943 ds = b->mM[sm]; is = sm; ys = 2;
945 for (i=0; i<=h-2; i++) {
948 if (b->mM[j+1] + lr > ds) {
949 ds = b->mM[j+1]; is = j+1; ys = 2;
953 for (i=h-2; i>=0; i--) {
956 if (b->mM[j-1] + lr > ds) {
957 ds = b->mM[j-1]; is = j-1; ys = 2;
961 if (b->mM[tm] + lr > ds) {
962 ds = b->mM[tm]; is = tm; ys = 2;
967 if (b->mM[is+1] + lr > b->mM[is]) is++;
972 if (b->mm[sm] < ds + lr) {
973 ds = b->mm[sm]; is = sm; ys = 2;
975 for (i=0; i<=h-2; i++) {
978 if (b->mm[j+1] < ds + lr) {
979 ds = b->mm[j+1]; is = j+1; ys = 2;
983 for (i=h-2; i>=0; i--) {
986 if (b->mm[j-1] < ds + lr) {
987 ds = b->mm[j-1]; is = j-1; ys = 2;
991 if (b->mm[tm] < ds + lr) {
992 ds = b->mm[tm]; is = tm; ys = 2;
997 if (b->mm[is+1] < b->mm[is] + lr) is++;
1006 ////////////////////////////////////////////////////////////
1007 // search the MB of t
1008 ////////////////////////////////////////////////////////////
1011 ts = max(SBid(MBfirst(t)),SBid(s)+1);
1012 il = SBid(SBfirst(t)-1);
1014 j = rmq_MB(b,ts,il,opt,&dm);
1015 //if (dm < ds + lr) {
1017 ds = dm; is = j; ys = 1;
1020 ////////////////////////////////////////////////////////////
1021 // search the SB of t
1022 ////////////////////////////////////////////////////////////
1025 j = rmq_SB(b,i,t,opt,&dm);
1026 d = depth(b,i) + dm;
1027 if (opt & OPT_MAX) {
1029 ds = d; is = j; ys = 0;
1033 ds = d; is = j; ys = 0;
1037 ////////////////////////////////////////////////////////////
1039 ////////////////////////////////////////////////////////////
1042 ss = SBid(is << logMB);
1043 il = SBid(MBlast(is << logMB));
1044 if (opt & OPT_MAX) {
1049 j = rmq_MB(b,ss,il,opt,&dm);
1050 ds = dm; is = j; ys = 1;
1055 il = SBlast(is << logSB);
1056 is = rmq_SB(b,ss,il,opt,&dm);
1057 //ds = depth(b,ss) + dm; ys = 0;