CC=gcc
-OBJECTS_BP=bp.o darray.o bpcore.o
+OBJECTS_BP=bp.o bp-utils.o bp-darray.o bp-core.o
LIB_BP=libbp.a
all: depend $(LIB_BP)
--- /dev/null
+#include "bp.h"
+#include "bp-utils.h"
+
+
+#define NOTFOUND -2
+#define CONTINUE -3
+#define END -4
+#define FOUND -5
+
+#define SBid(i) ((i)>>logSB)
+#define SBfirst(i) ((i) & (~(SB-1)))
+#define SBlast(i) ((i) | (SB-1))
+
+#define MBid(i) ((i)>>logMB)
+#define MBfirst(i) ((i) & (~(MB-1)))
+#define MBlast(i) ((i) | (MB-1))
+#define max(a,b) \
+ ({ __typeof__ (a) _a = (a); \
+ __typeof__ (b) _b = (b); \
+ _a > _b ? _a : _b; })
+
+
+#define min(a,b) \
+ ({ __typeof__ (a) _a = (a); \
+ __typeof__ (b) _b = (b); \
+ _a <= _b ? _a : _b; })
+
+
+pb getpat_preorder(pb *b)
+{
+ return *b;
+}
+
+pb getpat_postorder(pb *b)
+{
+ return ~(*b);
+}
+
+pb getpat_leaf(pb *b)
+{
+ pb w1,w2,w;
+ w1 = b[0];
+ w2 = (w1 << 1) + (b[1] >> (D-1));
+ w = w1 & (~w2);
+ return w;
+}
+
+pb getpat_inorder(pb *b)
+{
+ pb w1,w2,w;
+ w1 = b[0];
+ w2 = (w1 << 1) + (b[1] >> (D-1));
+ w = (~w1) & w2;
+ return w;
+}
+
+pb getpat_dfuds_leaf(pb *b)
+{
+ pb w1,w2,w;
+ w1 = b[0];
+ w2 = (w1 << 1) + (b[1] >> (D-1));
+ w = (~w1) & (~w2);
+ return w;
+}
+
+
+
+///////////////////////////////////////////
+// bp_depth(bp *b, int s)
+// returns the depth of s
+// The root node has depth 1
+///////////////////////////////////////////
+int bp_depth(bp *b, int s)
+{
+ int d;
+ if (s < 0) return 0;
+ d = 2 * bp_darray_rank(b->da,s) - (s+1);
+ return d;
+}
+
+int fast_depth(bp *b, int s)
+{
+ int d;
+ darray *da;
+ int r,j;
+
+ s++;
+ if ((s & (RRR-1)) != 0) {
+ //printf("fast_depth:warning s=%d\n",s);
+ return bp_depth(b,s);
+ }
+ da = b->da;
+ //d = 2 * (da->rl[s>>logR] + da->rm[s>>logRR] + da->rs[s>>logRRR]) - s;
+
+ r = da->rl[s>>logR] + da->rm[s>>logRR];
+ j = (s>>logRRR) & (RR/RRR-1);
+ while (j > 0) {
+ r += da->rs[((s>>logRR)<<(logRR-logRRR))+j-1];
+ j--;
+ }
+ d = 2 * r - s;
+
+ return d;
+}
+
+int search_SB_r(bp *b, int i, int rel)
+{
+ int j,r,n,il;
+ pb *p,x,w;
+
+ n = b->n;
+ il = min((SBid(i) + 1) << logSB,n);
+ p = &b->B[i>>logD];
+ while (i<il) {
+ x = *p++;
+ j = i & (D-1);
+ x <<= j;
+ j = D-j;
+ while (j>0) {
+ w = (x >> (D-ETW)) & ((1<<ETW)-1);
+ if (rel >= -ETW && rel <= ETW) {
+ r = fwdtbl[((rel+ETW)<<ETW)+w];
+ if (r<ETW && r<j) {
+ if (i+r >= n) return NOTFOUND;
+ return i+r;
+ }
+ }
+ r = min(j,ETW);
+ rel -= 2*popcount(w)-r;
+ x <<= r;
+ i += r;
+ j -= r;
+ }
+ }
+ return CONTINUE;
+}
+
+int search_MB_r(bp *b, int i, int td)
+{
+ int il,d;
+ int m,M,n;
+ pb *B;
+
+ B = b->B;
+ n = b->n;
+
+ il = min((MBid(i) + 1) << logMB,n);
+ for ( ; i < il; i+=SB) {
+#if (SB % RRR != 0)
+ d = bp_depth(b,i-1);
+#else
+ d = fast_depth(b,i-1);
+#endif
+ m = d + b->sm[SBid(i)] - SB;
+ M = d + b->sM[SBid(i)] - 1;
+ if (m <= td && td <= M) {
+ return search_SB_r(b,i,td-d);
+ }
+ }
+ if (i >= n) return NOTFOUND;
+ return CONTINUE;
+}
+
+///////////////////////////////////////////
+// bp_fwd_excess(bp *b,int s, int rel)
+// find the leftmost value depth(s)+rel to the right of s (exclusive)
+///////////////////////////////////////////
+int bp_fwd_excess(bp *b,int s, int rel)
+{
+ int i,n;
+ int d,td;
+ int m,M;
+ int m_ofs;
+ pb *B;
+ n = b->n; B = b->B;
+
+ i = s+1;
+
+ d = search_SB_r(b,i,rel);
+ if (d >= NOTFOUND) return d;
+
+ i = min((SBid(i) + 1) << logSB, n);
+ td = bp_depth(b,s) + rel;
+ d = search_MB_r(b,i,td);
+ if (d >= NOTFOUND) return d;
+
+ m_ofs = b->m_ofs;
+
+ i = MBid(s) + m_ofs;
+
+ while (i > 0) {
+ if ((i&1) == 0) {
+ i++;
+ m = b->mm[i];
+ M = b->mM[i];
+
+ if (m <= td && td <= M) {
+ while (i < m_ofs) {
+ i <<= 1;
+ m = b->mm[i];
+ M = b->mM[i];
+ if (!(m <= td && td <= M)) i++;
+ }
+ i -= m_ofs;
+ i <<= logMB;
+
+ return search_MB_r(b,i,td);
+ };
+
+ }
+ i >>= 1;
+ }
+ return NOTFOUND;
+}
+
+
+int degree_SB_slow(bp *b, int i, int t, int rel, int *ans, int ith)
+{
+ int j,r,n,il;
+ pb *p,x,w,w2;
+ int d, deg, v;
+
+ n = t;
+ il = min((SBid(i) + 1) << logSB,n);
+ d = deg = 0;
+
+ while (i < il) {
+ if (bp_getbit(b->B,i)==OP) {
+ d++;
+ } else {
+ d--;
+ }
+ if (d < rel) { // reached the end
+ if (ith > 0) {
+ return NOTFOUND;
+ } else {
+ *ans = deg;
+ return END;
+ }
+ }
+ if (d == rel) { // found the same depth
+ deg++;
+ if (deg == ith) {
+ *ans = i;
+ return FOUND;
+ }
+ }
+ i++;
+ }
+ *ans = deg;
+ return CONTINUE;
+}
+
+int degree_SB(bp *b, int i, int t, int rel, int *ans, int ith)
+{
+ int j,r,n,il;
+ pb *p,x,w,w2;
+ int d, deg, v;
+ int degtmp,itmp;
+ int ith2,d2;
+
+ n = t;
+ il = min((SBid(i) + 1) << logSB,n);
+ d = deg = 0;
+
+ p = &b->B[i>>logD];
+ while (i < il) {
+ x = *p++;
+ j = i & (D-1);
+ x <<= j;
+ j = min(D-j,il-i);
+ while (j>0) {
+ w = (x >> (D-ETW)) & ((1<<ETW)-1);
+ w2 = 0;
+ if (j < ETW || il-i < ETW) {
+ r = max(ETW-j,ETW-(il-i));
+ w2 = (1<<r)-1;
+ }
+ v = minmaxtbl_v[0][w | w2];
+ if (d + v < rel) {
+ if (ith > 0) {
+#if 0
+ for (r = 0; r < ETW; r++) {
+ if (w & 0x80) {
+ d++;
+ } else {
+ d--;
+ if (d < rel) break;
+ if (d == rel) {
+ ith--;
+ if (ith == 0) {
+ *ans = i + r;
+ return FOUND;
+ }
+ }
+ }
+ w <<= 1;
+ }
+ return NOTFOUND;
+#else
+ r = childtbl2[rel-d+ETW][ith-1][w];
+ if (r >= 0) {
+ *ans = i + r;
+ return FOUND;
+ }
+ return NOTFOUND;
+#endif
+ }
+ r = ETW-1-minmaxtbl_i[0][w | w2];
+ w2 = (1<<r)-1;
+ deg += degtbl2[((rel-d+ETW)<<ETW) + (w & (~w2))];
+ *ans = deg;
+ return END;
+ }
+ if (d + v == rel) {
+ r = degtbl[w | w2];
+ deg += r;
+ if (ith > 0) {
+ if (ith <= r) {
+ *ans = i + childtbl[((ith-1)<<ETW) + (w | w2)];
+ return FOUND;
+ }
+ ith -= r;
+ }
+ }
+
+ r = min(j,ETW);
+ d += 2*popcount(w)-r;
+ x <<= r;
+ i += r;
+ j -= r;
+ }
+ }
+
+ *ans = deg;
+ return CONTINUE;
+}
+
+int degree_MB(bp *b, int i, int t, int td, int *ans,int ith)
+{
+ int il,d;
+ int m,M,n,r;
+ pb *B;
+ int deg,degtmp;
+
+ d = 0;
+ B = b->B;
+ n = t;
+
+ il = min((MBid(i) + 1) << logMB,n);
+ deg = 0;
+ for ( ; i+SB-1 < il; i+=SB) {
+#if (SB % RRR != 0)
+ d = bp_depth(b,i-1);
+#else
+ d = fast_depth(b,i-1);
+#endif
+ m = d + b->sm[SBid(i)] - SB;
+ if (m < td) {
+ r = degree_SB(b,i,n,td-d,°tmp,ith);
+ if (ith > 0) {
+ if (r == NOTFOUND) return NOTFOUND;
+ *ans = degtmp;
+ return FOUND;
+ } else {
+ *ans = deg + degtmp;
+ return END;
+ }
+ }
+ if (m == td) {
+ if (ith > 0) {
+ if (ith <= b->sd[SBid(i)]) break;
+ ith -= b->sd[SBid(i)];
+ }
+ deg += b->sd[SBid(i)];
+ }
+ }
+ if (i < il) {
+#if (SB % RRR != 0)
+ d = bp_depth(b,i-1);
+#else
+ d = fast_depth(b,i-1);
+#endif
+ r = degree_SB(b,i,n,td-d,°tmp,ith);
+ if (ith > 0) {
+ if (r == NOTFOUND) return NOTFOUND;
+ if (r == FOUND) {
+ *ans = degtmp;
+ return FOUND;
+ }
+ } else {
+ deg += degtmp;
+ }
+ }
+ *ans = deg;
+ if (i >= n) return END;
+ return CONTINUE;
+}
+
+static int partition_range(int s,int t)
+{
+ int i,j,h;
+
+ printf("partition [%d,%d] => ",s,t);
+ h = 1;
+ while (s <= t) {
+ if (s & h) {
+ if (s+h-1 <= t) {
+ printf("[%d,%d] ",s,s+h-1);
+ s += h;
+ }
+ } else {
+ if (s+h > t) break;
+ }
+ h <<= 1;
+ }
+ while (h > 0) {
+ if (s+h-1 <= t) {
+ printf("[%d,%d] ",s,s+h-1);
+ s += h;
+ }
+ h >>= 1;
+ }
+ printf("\n");
+}
+
+
+
+
+///////////////////////////////////////////
+// bp_fast_degree(bp *b,int s, int t, int ith)
+// returns the number of children of s, to the left of t
+// returns the position of (ith)-th child of s if (ith > 0)
+///////////////////////////////////////////
+int bp_fast_degree(bp *b,int s, int t, int ith)
+{
+ int i,j,n;
+ int d,td;
+ int m,M;
+ int m_ofs;
+ pb *B;
+ int deg,degtmp;
+ int sm,tm,ss,h;
+
+ n = t;
+ B = b->B;
+
+ deg = 0;
+
+ i = s+1;
+ if (i != SBfirst(i)) {
+ d = degree_SB(b,i,n,0,°tmp,ith);
+ if (ith > 0) {
+ if (d == NOTFOUND) return -1;
+ if (d == FOUND) return degtmp;
+ ith -= degtmp;
+ }
+ if (d == END) return degtmp;
+ deg += degtmp;
+ }
+
+ td = bp_depth(b,s);
+
+ i = SBid(i+SB-1) << logSB;
+
+ if (i != MBfirst(i)) {
+ d = degree_MB(b,i,n,td,°tmp,ith);
+ if (ith > 0) {
+ if (d == NOTFOUND) return -1;
+ if (d == FOUND) return degtmp;
+ ith -= degtmp;
+ deg += degtmp;
+ } else {
+ deg += degtmp;
+ if (d == END) {
+ return deg;
+ }
+ }
+ }
+
+#if 0
+ // sequential search
+
+ sm = MBid(i+MB-1);
+ tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
+
+ m_ofs = b->m_ofs;
+ sm += m_ofs; tm += m_ofs;
+ for (i=sm; i<=tm; i++) {
+ if (b->mm[i] < td) {
+ break;
+ }
+ if (b->mm[i] == td) {
+ if (ith > 0) {
+ if (ith <= b->md[i]) break;
+ ith -= b->md[i];
+ }
+ deg += b->md[i];
+ }
+ }
+ ss = i - m_ofs;
+#else
+ sm = MBid(i+MB-1);
+ tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
+
+ m_ofs = b->m_ofs;
+ sm += m_ofs; tm += m_ofs;
+ ss = sm;
+
+ //partition_range(sm,tm);
+
+ //printf("partition [%d,%d] => ",sm,tm);
+ h = 1;
+ while (sm <= tm) {
+ if (sm & h) {
+ if (sm+h-1 <= tm) {
+ //printf("[%d,%d] ",sm,sm+h-1);
+ j = sm / h;
+ if (b->mm[j] < td) {
+ h >>= 1;
+ break;
+ }
+ if (b->mm[j] == td) {
+ if (ith > 0) {
+ if (ith <= b->md[j]) {
+ h >>= 1;
+ break;
+ }
+ ith -= b->md[j];
+ }
+ deg += b->md[j];
+ }
+ sm += h;
+ }
+ } else {
+ if (sm+h > tm) break;
+ }
+ h <<= 1;
+ }
+ while (h > 0) {
+ if (sm+h-1 <= tm) {
+ //printf("[%d,%d] ",sm,sm+h-1);
+ j = sm / h;
+ if (ith > 0) {
+ if (b->mm[j] >= td) {
+ if (b->mm[j] == td) {
+ if (ith > b->md[j]) {
+ ith -= b->md[j];
+ sm += h;
+ } else {
+ deg += b->md[j];
+ }
+ } else {
+ sm += h;
+ }
+ }
+ } else {
+ if (b->mm[j] >= td) {
+ if (b->mm[j] == td) {
+ deg += b->md[j];
+ }
+ sm += h;
+ }
+ }
+ }
+ h >>= 1;
+ }
+ //printf("\n");
+ ss = sm;
+
+ ss -= m_ofs;
+
+#endif
+
+ ss <<= logMB;
+
+ d = degree_MB(b,ss,n,td,°tmp,ith);
+ if (ith > 0) {
+ if (d == NOTFOUND) return -1;
+ if (d == FOUND) return degtmp;
+ }
+ deg += degtmp;
+ if (d == END) return deg;
+ return deg;
+
+ // unexpected (bug)
+ printf("degree: ???\n");
+ return -99;
+
+}
+
+int search_SB_l(bp *b, int i, int rel)
+{
+ int j,r,il;
+ pb *p,x,w;
+
+ il = SBfirst(i);
+
+ p = &b->B[i>>logD];
+ while (i>=il) {
+ x = *p--;
+ j = (i & (D-1))+1;
+ x >>= D-j;
+ while (j>0) {
+ w = x & ((1<<ETW)-1);
+ if (rel >= -ETW && rel <= ETW) {
+ r = bwdtbl[((rel+ETW)<<ETW)+w];
+ if (r<ETW && r<j) {
+ if (i-r < 0) return NOTFOUND;
+ return i-r-1;
+ }
+ }
+ r = min(j,ETW);
+ rel += 2*popcount(w)-r;
+ x >>= r;
+ i -= r;
+ j -= r;
+ }
+
+ }
+ if (i < 0) return NOTFOUND;
+ return CONTINUE;
+}
+
+int search_MB_l(bp *b, int i, int td)
+{
+ int il,d;
+ int m,M;
+ pb *B;
+
+#if 0
+ if (i % SB != SB-1) {
+ printf("search_MB_l:error!!! i=%d SB=%d\n",i,i%SB);
+ }
+#endif
+ B = b->B;
+
+ il = MBfirst(i);
+ for ( ; i >= il; i-=SB) {
+#if (SB % RRR != 0)
+ d = bp_depth(b,i-SB);
+#else
+ d = fast_depth(b,i-SB);
+#endif
+ m = d + b->sm[SBid(i)] - SB;
+ M = d + b->sM[SBid(i)] - 1;
+ if (m <= td && td <= M) {
+#if (SB % RRR != 0)
+ d = bp_depth(b,i);
+#else
+ d = fast_depth(b,i);
+#endif
+ if (d == td) return i;
+ return search_SB_l(b,i,td-d);
+ }
+ }
+ return CONTINUE;
+}
+
+///////////////////////////////////////////
+// bwd_excess(bp *b,int s, int rel)
+// find the rightmost value depth(s)+rel to the left of s (exclusive)
+///////////////////////////////////////////
+int bp_bwd_excess(bp *b,int s, int rel)
+{
+ int i,n;
+ int d,td;
+ int m,M;
+ int m_ofs;
+ pb *B;
+ n = b->n; B = b->B;
+
+ i = s;
+ d = search_SB_l(b,i,rel);
+ if (d >= NOTFOUND) return d;
+
+ i = SBfirst(i) -1;
+
+ td = bp_depth(b,s) + rel;
+
+ d = search_MB_l(b,i,td);
+ if (d >= NOTFOUND) return d;
+
+ m_ofs = b->m_ofs;
+ i = (s>>logMB) + m_ofs;
+ while (i > 0) {
+ if ((i&1) == 1) {
+ i--;
+ m = b->mm[i];
+ M = b->mM[i];
+ if (m <= td && td <= M) break;
+ }
+ i >>= 1;
+ }
+ if (i == 0) {
+ if (td == 0) return -1;
+ else return NOTFOUND;
+ }
+
+ while (i < m_ofs) {
+ i = i*2 + 1;
+ m = b->mm[i];
+ M = b->mM[i];
+ if (!(m <= td && td <= M)) i--;
+ }
+ i -= m_ofs;
+ i = ((i+1)<<logMB)-1;
+
+ d = search_MB_l(b,i,td);
+ if (d >= NOTFOUND) return d;
+
+ // unexpected (bug)
+ printf("bwd_excess: ???\n");
+ return -99;
+
+}
+
+int rmq_SB(bp *b, int s, int t, int opt, int *dm)
+{
+ int i,d;
+ int is,ds;
+ pb *p,x,w,w2;
+ int j,v,r;
+ int lr;
+ int op;
+
+ lr = 0; if (opt & OPT_RIGHT) lr = 1;
+ op = opt & (OPT_RIGHT | OPT_MAX);
+
+ is = s; ds = d = 0;
+ i = s+1;
+
+#if SB >= ETW
+ p = &b->B[i>>logD];
+ while (i <= t) {
+ x = *p++;
+ j = i & (D-1);
+ x <<= j;
+ j = min(D-j,t-i+1);
+ while (j>0) {
+ w = (x >> (D-ETW)) & ((1<<ETW)-1);
+ w2 = 0;
+ if (j < ETW || t-i < ETW-1) {
+ r = max(ETW-j,ETW-1-(t-i));
+ w2 = (1<<r)-1;
+ }
+
+ if (op & OPT_MAX) {
+ v = minmaxtbl_v[op][w & (~w2)];
+ if (d + v + lr > ds) {
+ ds = d + v;
+ is = i + minmaxtbl_i[op][w & (~w2)];
+ }
+ } else {
+ v = minmaxtbl_v[op][w | w2];
+ if (d + v < ds +lr) {
+ ds = d + v;
+ is = i + minmaxtbl_i[op][w | w2];
+ }
+ }
+
+ r = min(j,ETW);
+ d += 2*popcount(w)-r;
+ x <<= r;
+ i += r;
+ j -= r;
+ }
+ }
+#else
+ while (i <= t) {
+ if (bp_getbit(b->B,i)==OP) {
+ d++;
+ if (op & OPT_MAX) {
+ if (d + lr > ds) {
+ ds = d; is = i;
+ }
+ }
+ } else {
+ d--;
+ if (!(op & OPT_MAX)) {
+ if (d < ds + lr) {
+ ds = d; is = i;
+ }
+ }
+ }
+ i++;
+ }
+#endif
+ *dm = ds;
+ return is;
+}
+
+int rmq_MB(bp *b, int s, int t, int opt, int *dm)
+{
+ int i,d,m;
+ int mi,md;
+ int lr;
+
+ lr = 0; if (opt & OPT_RIGHT) lr = 1;
+
+ md = *dm; mi = -1;
+ for (i = s; i <= t; i++) {
+#if (SB % RRR != 0)
+ d = bp_depth(b,(i<<logSB)-1);
+#else
+ d = fast_depth(b,(i<<logSB)-1);
+#endif
+ if (opt & OPT_MAX) {
+ m = d + b->sM[i] - 1;
+ if (m + lr > md) {
+ md = m; mi = i;
+ }
+ } else {
+ m = d + b->sm[i] - SB;
+ if (m < md + lr) {
+ md = m; mi = i;
+ }
+ }
+ }
+ *dm = md;
+ return mi;
+}
+
+
+
+
+///////////////////////////////////////////
+// bp_rmq(bp *b, int s, int t, int opt)
+// returns the index of leftmost/rightmost minimum/maximum value
+// in the range [s,t] (inclusive)
+// returns the leftmost minimum if opt == 0
+// the rightmost one if (opt & OPT_RIGHT)
+// the maximum if (opt & OPT_MAX)
+///////////////////////////////////////////
+int bp_rmq(bp *b, int s, int t, int opt)
+{
+ int ss, ts; // SB index of s and t
+ int sm, tm; // MB index of s and t
+ int ds; // current min value
+ int is; // current min index
+ int ys; // type of current min index
+ // 0: is is the index of min
+ // 1: is is the SB index
+ // 2: is is the MB index
+ int m_ofs;
+ int i,j,il,d,n;
+ int dm;
+ int lr;
+
+ lr = 0; if (opt & OPT_RIGHT) lr = 1;
+
+ n = b->n;
+ if (s < 0 || t >= n || s > t) {
+ printf("rmq: error s=%d t=%d n=%d\n",s,t,n);
+ return -1;
+ }
+ if (s == t) return s;
+
+
+ ////////////////////////////////////////////////////////////
+ // search the SB of s
+ ////////////////////////////////////////////////////////////
+
+ il = min(SBlast(s),t);
+ is = rmq_SB(b,s,il,opt,&dm);
+ if (il == t) { // scan reached the end of the range
+ return is;
+ }
+ ds = bp_depth(b,s) + dm; ys = 0;
+
+ ////////////////////////////////////////////////////////////
+ // search the MB of s
+ ////////////////////////////////////////////////////////////
+
+ ss = SBid(s) + 1;
+ il = min(SBid(MBlast(s)),SBid(t)-1);
+ dm = ds;
+ j = rmq_MB(b,ss,il,opt,&dm);
+ //if (dm < ds + lr) {
+ if (j >= 0) {
+ ds = dm; is = j; ys = 1;
+ }
+
+ ////////////////////////////////////////////////////////////
+ // search the tree
+ ////////////////////////////////////////////////////////////
+
+ sm = MBid(s) + 1;
+ tm = MBid(t) - 1;
+
+#if 0
+ // sequential search
+ m_ofs = b->m_ofs;
+ sm += m_ofs; tm += m_ofs;
+ for (i=sm; i<=tm; i++) {
+ if (opt & OPT_MAX) {
+ if (b->mM[i] + lr > ds) {
+ ds = b->mM[i]; is = i - m_ofs; ys = 2;
+ }
+ } else {
+ if (b->mm[i] < ds + lr) {
+ ds = b->mm[i]; is = i - m_ofs; ys = 2;
+ }
+ }
+ }
+
+#else
+ if (sm <= tm) {
+ int h;
+ h = blog(sm ^ tm);
+
+ m_ofs = b->m_ofs;
+ sm += m_ofs; tm += m_ofs;
+
+ if (opt & OPT_MAX) {
+ if (b->mM[sm] + lr > ds) {
+ ds = b->mM[sm]; is = sm; ys = 2;
+ }
+ for (i=0; i<=h-2; i++) {
+ j = sm>>i;
+ if ((j&1) == 0) {
+ if (b->mM[j+1] + lr > ds) {
+ ds = b->mM[j+1]; is = j+1; ys = 2;
+ }
+ }
+ }
+ for (i=h-2; i>=0; i--) {
+ j = tm>>i;
+ if ((j&1)==1) {
+ if (b->mM[j-1] + lr > ds) {
+ ds = b->mM[j-1]; is = j-1; ys = 2;
+ }
+ }
+ }
+ if (b->mM[tm] + lr > ds) {
+ ds = b->mM[tm]; is = tm; ys = 2;
+ }
+ if (ys == 2) {
+ while (is < m_ofs) {
+ is <<= 1;
+ if (b->mM[is+1] + lr > b->mM[is]) is++;
+ }
+ is -= m_ofs;
+ }
+ } else { // MIN
+ if (b->mm[sm] < ds + lr) {
+ ds = b->mm[sm]; is = sm; ys = 2;
+ }
+ for (i=0; i<=h-2; i++) {
+ j = sm>>i;
+ if ((j&1) == 0) {
+ if (b->mm[j+1] < ds + lr) {
+ ds = b->mm[j+1]; is = j+1; ys = 2;
+ }
+ }
+ }
+ for (i=h-2; i>=0; i--) {
+ j = tm>>i;
+ if ((j&1)==1) {
+ if (b->mm[j-1] < ds + lr) {
+ ds = b->mm[j-1]; is = j-1; ys = 2;
+ }
+ }
+ }
+ if (b->mm[tm] < ds + lr) {
+ ds = b->mm[tm]; is = tm; ys = 2;
+ }
+ if (ys == 2) {
+ while (is < m_ofs) {
+ is <<= 1;
+ if (b->mm[is+1] < b->mm[is] + lr) is++;
+ }
+ is -= m_ofs;
+ }
+ }
+ }
+
+#endif
+
+ ////////////////////////////////////////////////////////////
+ // search the MB of t
+ ////////////////////////////////////////////////////////////
+
+
+ ts = max(SBid(MBfirst(t)),SBid(s)+1);
+ il = SBid(SBfirst(t)-1);
+ dm = ds;
+ j = rmq_MB(b,ts,il,opt,&dm);
+ //if (dm < ds + lr) {
+ if (j >= 0) {
+ ds = dm; is = j; ys = 1;
+ }
+
+ ////////////////////////////////////////////////////////////
+ // search the SB of t
+ ////////////////////////////////////////////////////////////
+
+ i = SBfirst(t);
+ j = rmq_SB(b,i,t,opt,&dm);
+ d = bp_depth(b,i) + dm;
+ if (opt & OPT_MAX) {
+ if (d + lr > ds) {
+ ds = d; is = j; ys = 0;
+ }
+ } else {
+ if (d < ds + lr) {
+ ds = d; is = j; ys = 0;
+ }
+ }
+
+ ////////////////////////////////////////////////////////////
+ // search the rest
+ ////////////////////////////////////////////////////////////
+
+ if (ys == 2) {
+ ss = SBid(is << logMB);
+ il = SBid(MBlast(is << logMB));
+ if (opt & OPT_MAX) {
+ dm = -n-1;
+ } else {
+ dm = n+1;
+ }
+ j = rmq_MB(b,ss,il,opt,&dm);
+ ds = dm; is = j; ys = 1;
+ }
+
+ if (ys == 1) {
+ ss = is << logSB;
+ il = SBlast(is << logSB);
+ is = rmq_SB(b,ss,il,opt,&dm);
+ //ds = bp_depth(b,ss) + dm; ys = 0;
+ }
+
+ return is;
+}
+
--- /dev/null
+#include <stdio.h>\r
+#include <stdlib.h>\r
+#include "bp-darray.h"\r
+#include "bp-utils.h"\r
+\r
+#define PBS (sizeof(pb)*8)\r
+#define D (1<<logD)\r
+#define logM 5\r
+#define M (1<<logM)\r
+#define logP 8\r
+#define P (1<<logP)\r
+#define logLL 16 // size of word\r
+#define LL (1<<logLL)\r
+#define logLLL 7\r
+//#define logLLL 5\r
+#define LLL (1<<logLLL)\r
+//#define logL 10\r
+//#define logL (logLL-3)\r
+#define logL (logLL-1-5)\r
+#define L (1<<logL)\r
+\r
+#define max(a,b) \\r
+ ({ __typeof__ (a) _a = (a); \\r
+ __typeof__ (b) _b = (b); \\r
+ _a > _b ? _a : _b; })\r
+\r
+\r
+#define min(a,b) \\r
+ ({ __typeof__ (a) _a = (a); \\r
+ __typeof__ (b) _b = (b); \\r
+ _a <= _b ? _a : _b; })\r
+\r
+\r
+\r
+\r
+static dword getbits(pb *B, int i, int d)\r
+{\r
+ dword j,x;\r
+\r
+ x = 0;\r
+ for (j=0; j<d; j++) {\r
+ x <<= 1;\r
+ x += bp_getbit(B,i+j);\r
+ }\r
+ return x;\r
+}\r
+\r
+static int getpattern(pb *B, int i, int k, pb pat)\r
+{\r
+ int j;\r
+ int x;\r
+ x = 1;\r
+ for (j=0; j<k; j++) {\r
+ x &= bp_getbit(B,i+j) ^ (~(pat>>(k-1-j)));\r
+ }\r
+ //printf("getpattern(%d) = %d\n",i,x);\r
+ return x;\r
+}\r
+\r
+static int selecttbl[8*256];\r
+static int selecttbl_init = 0;\r
+\r
+static void make_selecttbl(void)\r
+{\r
+ int i,x,r;\r
+ pb buf[1];\r
+ if (selecttbl_init) return;\r
+\r
+ selecttbl_init = 1;\r
+\r
+ for (x = 0; x < 256; x++) {\r
+ bp_setbits(buf,0,8,x);\r
+ for (r=0; r<8; r++) selecttbl[(r<<8)+x] = -1;\r
+ r = 0;\r
+ for (i=0; i<8; i++) {\r
+ if (bp_getbit(buf,i)) {\r
+ selecttbl[(r<<8)+x] = i;\r
+ r++;\r
+ }\r
+ }\r
+ }\r
+}\r
+\r
+\r
+darray * bp_darray_construct(int n, pb *buf, int opt)\r
+{\r
+ int i,j,k,m;\r
+ int nl;\r
+ int p,pp;\r
+ int il,is,ml,ms;\r
+ int r,m2;\r
+ int p1,p2,p3,p4,s1,s2,s3,s4;\r
+ darray * da;\r
+\r
+ bp_xmalloc(da, 1);\r
+\r
+ da->idx_size = 0;\r
+\r
+ make_selecttbl();\r
+\r
+\r
+ if (L/LLL == 0) {\r
+ printf("ERROR: L=%d LLL=%d\n",L,LLL);\r
+ exit(1);\r
+ }\r
+\r
+ m = 0;\r
+ for (i=0; i<n; i++) m += bp_getbit(buf,i);\r
+ da->n = n;\r
+ da->m = m;\r
+ //printf("n=%d m=%d\n",n,m);\r
+\r
+ da->buf = buf;\r
+\r
+ if (opt & (~OPT_NO_RANK)) { // construct select table\r
+\r
+ nl = (m-1) / L + 1;\r
+ bp_xmalloc(da->lp, nl+1);\r
+ da->idx_size += (nl+1) * sizeof(*da->lp);\r
+ bp_xmalloc(da->p, nl+1);\r
+ da->idx_size += (nl+1) * sizeof(*da->p);\r
+\r
+ for (r = 0; r < 2; r++) {\r
+ s1 = s2 = s3 = s4 = 0;\r
+ p1 = p2 = p3 = p4 = -1;\r
+\r
+ ml = ms = 0;\r
+ for (il = 0; il < nl; il++) {\r
+\r
+ while (s1 <= il*L) {\r
+ if (bp_getbit(buf,p1+1)) s1++;\r
+ p1++;\r
+ }\r
+ pp = p1;\r
+ da->lp[il] = pp;\r
+ i = min((il+1)*L-1,m-1);\r
+\r
+ while (s2 <= i) {\r
+ if (bp_getbit(buf,p2+1)) s2++;\r
+ p2++;\r
+ }\r
+ p = p2;\r
+\r
+ if (p - pp >= LL) {\r
+ if (r == 1) {\r
+ for (is = 0; is < L; is++) {\r
+ if (il*L+is >= m) break;\r
+\r
+ while (s3 <= il*L+is) {\r
+ if (bp_getbit(buf,p3+1)) s3++;\r
+ p3++;\r
+ }\r
+ da->sl[ml*L+is] = p3;\r
+ }\r
+ }\r
+ da->p[il] = -(ml+1);\r
+ ml++;\r
+ } else {\r
+ if (r == 1) {\r
+ for (is = 0; is < L/LLL; is++) {\r
+ if (il*L+is*LLL >= m) break;\r
+ while (s4 <= il*L+is*LLL) {\r
+ if (bp_getbit(buf,p4+1)) s4++;\r
+ p4++;\r
+ }\r
+\r
+ da->ss[ms*(L/LLL)+is] = p4 - pp;\r
+ }\r
+ }\r
+ da->p[il] = ms;\r
+ ms++;\r
+ }\r
+ }\r
+ if (r == 0) {\r
+ bp_xmalloc(da->sl,ml*L+1);\r
+ da->idx_size += (ml*L+1) * sizeof(*da->sl);\r
+ bp_xmalloc(da->ss, ms*(L/LLL)+1);\r
+ da->idx_size += (ms*(L/LLL)+1) * sizeof(*da->ss);\r
+ }\r
+ }\r
+\r
+ } else { // no select table\r
+ da->lp = NULL;\r
+ da->p = da->sl = NULL;\r
+ da->ss = NULL;\r
+ }\r
+\r
+ // construct rank table\r
+\r
+ if ((opt & OPT_NO_RANK) == 0) {\r
+ bp_xmalloc(da->rl,n/R1+2);\r
+ da->idx_size += (n/R1+2) * sizeof(*da->rl);\r
+ bp_xmalloc(da->rm,n/RR+2);\r
+\r
+ da->idx_size += (n/RR+2) * sizeof(*da->rm);\r
+\r
+ bp_xmalloc(da->rs,n/RRR+2);\r
+ da->idx_size += (n/RRR+2) * sizeof(*da->rs);\r
+\r
+ r = 0;\r
+ for (k=0; k<=n; k+=R1) {\r
+ da->rl[k/R1] = r;\r
+ m2 = 0;\r
+ for (i=0; i<R1; i+=RR) {\r
+ if (k+i <= n) da->rm[(k+i)/RR] = m2;\r
+ m = 0;\r
+ for (j=0; j<RR; j++) {\r
+ if (k+i+j < n && bp_getbit(buf,k+i+j)==1) m++;\r
+ if (j % RRR == RRR-1) {\r
+ if (k+i+j <= n) da->rs[(k+i+j)/RRR] = m;\r
+ m2 += m;\r
+ m = 0;\r
+ }\r
+ }\r
+ if (m != 0) {\r
+ printf("???\n");\r
+ }\r
+ }\r
+ r += m2;\r
+ }\r
+ }\r
+\r
+ return da;\r
+}\r
+\r
+\r
+void bp_darray_free(darray *da) {\r
+\r
+ if (!da) return;\r
+ if (da->buf) bp_free(da->buf);\r
+ if (da->lp) bp_free(da->lp);\r
+ if (da->sl) bp_free(da->sl);\r
+ if (da->ss) bp_free(da->ss);\r
+ if (da->p) bp_free(da->p);\r
+ if (da->rl) bp_free(da->rl);\r
+ if (da->rm) bp_free(da->rm);\r
+ if (da->rs) bp_free(da->rs);\r
+ bp_free(da);\r
+\r
+}\r
+\r
+static int darray_rank0(darray *da, int i)\r
+{\r
+ int r,j;\r
+ pb *p;\r
+\r
+#if (RRR == D*2)\r
+ r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];\r
+ p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
+ j = i & (RRR-1);\r
+ if (j < D) r += popcount(*p >> (D-1-j));\r
+ else r += popcount(*p) + popcount(p[1] >> (D-1-(j-D)));\r
+#else\r
+\r
+ j = i & (RRR-1);\r
+ if (j < RRR/2) {\r
+ r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];\r
+ p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
+ while (j >= D) {\r
+ r += popcount(*p++);\r
+ j -= D;\r
+ }\r
+ r += popcount(*p >> (D-1-j));\r
+ } else {\r
+ j = RRR-1 - (i & (RRR-1));\r
+ i += j+1;\r
+ r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];\r
+ p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
+ while (j >= D) {\r
+ r -= popcount(*--p);\r
+ j -= D;\r
+ }\r
+ if (j > 0) r -= popcount(*--p << (D-j));\r
+ }\r
+\r
+#endif\r
+\r
+ return r;\r
+}\r
+\r
+int bp_darray_rank(darray *da, int i)\r
+{\r
+ int r,j,i_rr, i_rrr;\r
+ pb *p;\r
+ i_rr = i >> logRR;\r
+ i_rrr = i >> logRRR;\r
+ r = da->rl[i>>logR] + da->rm[i_rr];\r
+\r
+ j = (i_rrr) & (RR/RRR-1);\r
+ while (j > 0) {\r
+ r += da->rs[((i_rr)<<(logRR-logRRR))+j-1];\r
+ j--;\r
+ }\r
+\r
+ p = da->buf + ((i_rrr)<<(logRRR-logD));\r
+ j = i & (RRR-1);\r
+ while (j >= D) {\r
+ r += popcount(*p++);\r
+ j -= D;\r
+ }\r
+ r += popcount(*p >> (D-1-j));\r
+\r
+ return r;\r
+}\r
+\r
+int bp_darray_select_bsearch(darray *da, int i, pb (*getpat)(pb *))\r
+{\r
+ int j;\r
+ int l,r,m,n;\r
+ int s;\r
+ int t,x,rr;\r
+ pb *p,w;\r
+\r
+\r
+ if (i == 0) return -1;\r
+\r
+ if (i > da->m) {\r
+ return -1;\r
+ }\r
+\r
+ i--;\r
+\r
+\r
+\r
+ n = da->m;\r
+\r
+ t = i;\r
+\r
+ l = 0; r = (n-1)>>logR;\r
+ // find the smallest index x s.t. rl[x] >= t\r
+ while (l < r) {\r
+ m = (l+r)/2;\r
+ //printf("m=%d rl[m+1]=%d t=%d\n",m,da->rl[m+1],t);\r
+ if (da->rl[m+1] > t) { // m+1 is out of range\r
+ r = m; // new r = m >= l\r
+ } else {\r
+ l = m+1; // new l = m+1 <= r\r
+ }\r
+ }\r
+ x = l;\r
+ t -= da->rl[x];\r
+\r
+ x <<= logR;\r
+\r
+ l = x >> logRR; r = (min(x+R1-1,n))>>logRR;\r
+ while (l < r) {\r
+ m = (l+r)/2;\r
+ //printf("m=%d rm[m+1]=%d t=%d\n",m,da->rm[m+1],t);\r
+ if (da->rm[m+1] > t) { // m+1 is out of range\r
+ r = m;\r
+ } else {\r
+ l = m+1; // new l = m+1 <= r\r
+ }\r
+ }\r
+ x = l;\r
+ t -= da->rm[x];\r
+\r
+ x <<= logRR;\r
+\r
+#if 0\r
+ l = x >> logRRR; r = (min(x+RR-1,n))>>logRRR;\r
+ while (l < r) {\r
+ m = (l+r)/2;\r
+ //printf("m=%d rs[m+1]=%d t=%d\n",m,da->rs[m+1],t);\r
+ if (da->rs[m+1] > t) { // m+1 is out of range\r
+ r = m;\r
+ } else {\r
+ l = m+1; // new l = m+1 <= r\r
+ }\r
+ }\r
+ x = l;\r
+ t -= da->rs[x];\r
+#else\r
+ l = x >> logRRR;\r
+ while (t > da->rs[l]) {\r
+ t -= da->rs[l];\r
+ l++;\r
+ }\r
+ x = l;\r
+#endif\r
+\r
+ x <<= logRRR;\r
+\r
+ p = &da->buf[x >> logD];\r
+ while (1) {\r
+ m = popcount(getpat(p));\r
+ if (m > t) break;\r
+ t -= m;\r
+ x += D;\r
+ p++;\r
+ }\r
+\r
+ w = getpat(p);\r
+ while (1) {\r
+ rr = popcount8(w >> (D-8));\r
+ if (rr > t) break;\r
+ t -= rr;\r
+ x += 8;\r
+ w <<= 8;\r
+ }\r
+ x += selecttbl[((t-0)<<8)+(w>>(D-8))];\r
+\r
+#if 0\r
+ if (x != s) {\r
+ printf("error x=%d s=%d\n",x,s);\r
+ }\r
+#endif\r
+ return x;\r
+}\r
+\r
+int bp_darray_pat_rank(darray *da, int i, pb (*getpat)(pb *))\r
+{\r
+ int r,j;\r
+ pb *p;\r
+\r
+ r = da->rl[i>>logR] + da->rm[i>>logRR];\r
+ j = (i>>logRRR) & (RR/RRR-1);\r
+ while (j > 0) {\r
+ r += da->rs[((i>>logRR)<<(logRR-logRRR))+j-1];\r
+ j--;\r
+ }\r
+\r
+ p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
+ j = i & (RRR-1);\r
+ while (j >= D) {\r
+ r += popcount(getpat(p));\r
+ p++;\r
+ j -= D;\r
+ }\r
+ r += popcount(getpat(p) >> (D-1-j));\r
+\r
+ return r;\r
+}\r
+\r
+\r
+int bp_darray_select(darray *da, int i,int f)\r
+{\r
+ int p,r;\r
+ int il;\r
+ int rr;\r
+ pb x;\r
+ pb *q;\r
+\r
+ if (i == 0) return -1;\r
+\r
+ if (i > da->m) {\r
+ return -1;\r
+ //printf("ERROR: m=%d i=%d\n",da->m,i);\r
+ //exit(1);\r
+ }\r
+\r
+ i--;\r
+\r
+ il = da->p[i>>logL];\r
+ if (il < 0) {\r
+ il = -il-1;\r
+ p = da->sl[(il<<logL)+(i & (L-1))];\r
+ } else {\r
+ p = da->lp[i>>logL];\r
+ p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];\r
+ r = i - (i & (LLL-1));\r
+\r
+ q = &(da->buf[p>>logD]);\r
+\r
+ if (f == 1) {\r
+ rr = p & (D-1);\r
+ r -= popcount(*q >> (D-1-rr));\r
+ p = p - rr;\r
+\r
+ while (1) {\r
+ rr = popcount(*q);\r
+ if (r + rr >= i) break;\r
+ r += rr;\r
+ p += D;\r
+ q++;\r
+ }\r
+\r
+ x = *q;\r
+ while (1) {\r
+ //rr = popcount(x >> (D-8));\r
+ //rr = popcount(x >> (D-8));\r
+ rr = popcount8(x >> (D-8));\r
+ if (r + rr >= i) break;\r
+ r += rr;\r
+ p += 8;\r
+ x <<= 8;\r
+ }\r
+ p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];\r
+ } else {\r
+ rr = p & (D-1);\r
+ r -= popcount((~(*q)) >> (D-1-rr));\r
+ p = p - rr;\r
+\r
+ while (1) {\r
+ rr = popcount(~(*q));\r
+ if (r + rr >= i) break;\r
+ r += rr;\r
+ p += D;\r
+ q++;\r
+ }\r
+\r
+ x = ~(*q);\r
+\r
+ while (1) {\r
+ rr = popcount8(x >> (D-8));\r
+ if (r + rr >= i) break;\r
+ r += rr;\r
+ p += 8;\r
+ x <<= 8;\r
+ }\r
+ p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];\r
+ }\r
+ }\r
+ return p;\r
+}\r
+\r
+int bp_darray_pat_select(darray *da, int i, pb (*getpat)(pb *))\r
+{\r
+ int p,r;\r
+ int il;\r
+ int rr;\r
+ pb x;\r
+ pb *q;\r
+\r
+ if (i == 0) return -1;\r
+\r
+ if (i > da->m) {\r
+ return -1;\r
+ //printf("ERROR: m=%d i=%d\n",da->m,i);\r
+ //exit(1);\r
+ }\r
+\r
+ i--;\r
+\r
+ il = da->p[i>>logL];\r
+ if (il < 0) {\r
+ il = -il-1;\r
+ p = da->sl[(il<<logL)+(i & (L-1))];\r
+ } else {\r
+ p = da->lp[i>>logL];\r
+ p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];\r
+ r = i - (i & (LLL-1));\r
+\r
+ q = &(da->buf[p>>logD]);\r
+\r
+ rr = p & (D-1);\r
+ r -= popcount(getpat(q) >> (D-1-rr));\r
+ p = p - rr;\r
+\r
+ while (1) {\r
+ rr = popcount(getpat(q));\r
+ if (r + rr >= i) break;\r
+ r += rr;\r
+ p += D;\r
+ q++;\r
+ }\r
+\r
+ x = getpat(q);\r
+ while (1) {\r
+ rr = popcount8(x >> (D-8));\r
+ if (r + rr >= i) break;\r
+ r += rr;\r
+ p += 8;\r
+ x <<= 8;\r
+ }\r
+ p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];\r
+ }\r
+ return p;\r
+}\r
+\r
+\r
+darray * bp_darray_pat_construct(int n, pb *buf, int k, pb pat, int opt)\r
+{\r
+ int i;\r
+ pb *b;\r
+ darray *da;\r
+\r
+ bp_xmalloc(b, (n+D-1)/D);\r
+\r
+ for (i=0; i<n-k+1; i++) {\r
+ bp_setbit(b, i, getpattern(buf,i,k,pat));\r
+ }\r
+ for (i=n-k+1; i<n; i++) {\r
+ bp_setbit(b, i, 0);\r
+ }\r
+\r
+ da = bp_darray_construct(n, b, opt);\r
+ da->buf = buf;\r
+\r
+ bp_free(b);\r
+\r
+ return da;\r
+}\r
--- /dev/null
+#ifndef BP_DARRAY_H_
+#define BP_DARRAY_H_
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef unsigned char byte;
+typedef unsigned short word;
+typedef unsigned int dword;
+
+#define OPT_NO_RANK (1<<30)
+
+
+typedef dword pb;
+
+
+typedef struct {
+ int n,m;
+ int size;
+ pb *buf;
+ dword *lp;
+ dword *sl;
+ word *ss;
+ dword *p;
+
+ dword *rl;
+ word *rm;
+ byte *rs;
+
+ int idx_size;
+} darray;
+
+
+darray * bp_darray_construct(int n, pb *buf,int opt);
+void bp_darray_free(darray *da);
+
+int bp_darray_select(darray *da, int i,int f);
+int bp_darray_rank(darray *da, int i);
+darray * bp_darray_pat_construct(int n, pb *buf, int k, pb pat, int opt);
+int bp_darray_pat_select(darray *da, int i, pb (*getpat)(pb *));
+int bp_darray_pat_rank(darray *da, int i, pb (*getpat)(pb *));
+
+int bp_darray_select_bsearch(darray *da, int i, pb (*getpat)(pb *));
+
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif
--- /dev/null
+#include <stdio.h>
+#include "bp-utils.h"
+
+static size_t allocated = 0;
+
+void * bp_malloc(size_t n)
+{
+ void * res = malloc(n);
+ if (!res) {
+ fprintf(stderr, __FILE__ ": failure to allocate %lu bytes\n", n);
+ exit(1);
+ };
+ allocated += n;
+ return res;
+}
+
+void bp_free(void * p)
+{
+ if (p) free(p);
+}
+size_t bp_get_alloc_stats(void)
+{
+ return allocated;
+}
+void bp_reset_alloc_states(void)
+{
+ allocated = 0;
+}
+
+int bp_setbit(unsigned int *B, int i,int x)
+{
+ int j,l;
+
+ j = i / D;
+ l = i % D;
+ if (x==0) B[j] &= (~(1<<(D-1-l)));
+ else if (x==1) B[j] |= (1<<(D-1-l));
+ else {
+ printf("error setbit x=%d\n",x);
+ exit(1);
+ }
+ return x;
+}
+
+int bp_setbits(unsigned int *B, int i, int d, int x)
+{
+ int j;
+
+ for (j=0; j<d; j++) {
+ bp_setbit(B,i+j,(x>>(d-j-1))&1);
+ }
+ return x;
+}
--- /dev/null
+#ifndef BP_UTILS_H_
+#define BP_UTILS_H_
+
+#ifdef __cplusplus
+
+extern "C" {
+#endif
+
+#define logD 5
+#define D (1<<logD)
+
+#define logR 16
+#define R1 (1<<logR)
+#define logRR 10
+//#define logRR 8
+#define RR (1<<logRR)
+#define logRRR 7
+#define RRR (1<<logRRR)
+
+
+
+
+
+#include <stdlib.h>
+
+
+
+#ifdef HAS_NATIVE_POPCOUNT
+static inline unsigned int popcount(unsigned int n){
+ asm ("popcnt %1, %0" : "=r" (n) : "0" (n));
+ return n;
+}
+
+static inline unsigned int popcount8(unsigned int n) {
+ return popcount(n & 0xff);
+}
+
+#else
+
+static unsigned int popcount8(unsigned int x)
+{
+ unsigned int r;
+ r = x;
+ r = ((r & 0xaa)>>1) + (r & 0x55);
+ r = ((r & 0xcc)>>2) + (r & 0x33);
+ r = ((r>>4) + r) & 0x0f;
+ return r;
+}
+
+static inline unsigned int popcount(unsigned int x)
+{
+ unsigned int m1 = 0x55555555;
+ unsigned int m2 = 0xc30c30c3;
+ x -= (x >> 1) & m1;
+ x = (x & m2) + ((x >> 2) & m2) + ((x >> 4) & m2);
+ x += x >> 6;
+ return (x + (x >> 12) + (x >> 24)) & 0x3f;
+}
+
+#endif
+
+
+void * bp_malloc(size_t);
+
+#define bp_xmalloc(p, n) p = (__typeof__(p)) bp_malloc((n)*sizeof(*p))
+
+
+
+void bp_free(void *);
+
+size_t bp_alloc_stats(void);
+
+void bp_reset_alloc_states(void);
+
+int bp_setbit(unsigned int *B, int i,int x);
+
+int bp_setbits(unsigned int *B, int i, int d, int x);
+
+static inline int bp_getbit(unsigned int *B, int i)
+{
+ int j = i >> logD;
+ int l = i & (D-1);
+ return (B[j] >> (D-1-l)) & 1;
+}
+
+#ifdef __cplusplus
+ }
+#endif
+
+#endif
#include "bp.h"
+//#define CHECK
+#define RANDOM
+
+
#define max(a,b) \
({ __typeof__ (a) _a = (a); \
__typeof__ (b) _b = (b); \
__typeof__ (b) _b = (b); \
_a <= _b ? _a : _b; })
-//#define CHECK
-#define RANDOM
-int msize=0;
-#define mymalloc(p,n,f) { \
- p = (__typeof__(p)) malloc((n)*sizeof(*p)); \
-if ((p)==NULL) {printf("not enough memory (%d bytes) in line %d\n",msize,__LINE__); \
- exit(1);}; \
-msize += (f)*(n)*sizeof(*p); \
-;}
-int postorder_select_bsearch(bp *b,int s);
+static int postorder_select_bsearch(bp *b,int s);
-int naive_depth(bp *b, int s)
+static int naive_depth(bp *b, int s)
{
int i,d;
if (s < 0) return 0;
d = 0;
for (i=0; i<=s; i++) {
- if (getbit(b->B,i)==OP) {
+ if (bp_getbit(b->B,i)==OP) {
d++;
} else {
d--;
return d;
}
-void printbp(bp *b, int s, int t)
+void bp_print(bp *b, int s, int t)
{
int i,c,d;
d = 0;
for (i=s; i<=t; i++) {
- if (getbit(b->B,i)==OP) {
+ if (bp_getbit(b->B,i)==OP) {
c = '(';
d++;
} else {
int i,v;
int *stack,s;
- mymalloc(matchtbl,n,0);
- mymalloc(parenttbl,n,0);
- mymalloc(stack,n,0);
+ bp_xmalloc(matchtbl,n);
+ bp_xmalloc(parenttbl,n);
+ bp_xmalloc(stack,n);
for (i=0; i<n; i++) matchtbl[i] = parenttbl[i] = -1;
v = 0;
stack[0] = -1;
for (i=0; i<n; i++) {
- if (getbit(B,i)==OP) {
+ if (bp_getbit(B,i)==OP) {
v++;
if (v > 0) {
parenttbl[i] = stack[v-1];
free(stack);
}
-int popCount[1<<ETW];
+
int fwdtbl[(2*ETW+1)*(1<<ETW)];
int bwdtbl[(2*ETW+1)*(1<<ETW)];
#if 0
int childtbl2[2*ETW+1][ETW][(1<<ETW)];
int depthtbl[(2*ETW+1)*(1<<ETW)];
int inited = 0;
+
void make_matchtbl(void)
{
int i,j,x,r;
return;
inited = 1;
for (x = 0; x < (1<<ETW); x++) {
- setbits(buf,0,ETW,x);
+ bp_setbits(buf,0,ETW,x);
for (r=-ETW; r<=ETW; r++) fwdtbl[((r+ETW)<<ETW)+x] = ETW;
for (r=-ETW; r<=ETW; r++) bwdtbl[((r+ETW)<<ETW)+x] = ETW;
for (r=-ETW; r<=ETW; r++) degtbl2[((r+ETW)<<ETW)+x] = 0;
r = 0;
for (i=0; i<ETW; i++) {
- if (getbit(buf,i)==OP) {
+ if (bp_getbit(buf,i)==OP) {
r++;
} else {
r--;
r = 0;
for (i=ETW-1; i>=0; i--) {
- if (getbit(buf,i)==OP) {
+ if (bp_getbit(buf,i)==OP) {
r--;
} else {
r++;
r = 0;
for (i=0; i<ETW; i++) {
- if (getbit(buf,i)==OP) {
+ if (bp_getbit(buf,i)==OP) {
r++;
} else {
r--;
depthtbl[((r+ETW)<<ETW)+x] += (1<<(ETW-1));
}
- r = 0;
- for (i=0; i<ETW; i++) {
- if (getbit(buf,i)==OP) r++;
- }
- popCount[x] = r;
r = 0; m = 0; M = 0;
m = ETW+1; M = -ETW-1;
minmaxtbl_v[OPT_MIN | OPT_LEFT][x] = ETW+1;
deg = 0;
for (i=0; i<ETW; i++) {
- if (getbit(buf,i)==OP) {
+ if (bp_getbit(buf,i)==OP) {
r++;
if (r > M) {
M = r;
minmaxtbl_v[OPT_MAX | OPT_RIGHT][x] = -ETW-1;
minmaxtbl_v[OPT_MIN | OPT_RIGHT][x] = ETW+1;
for (i=0; i<ETW; i++) {
- if (getbit(buf,i)==OP) {
+ if (bp_getbit(buf,i)==OP) {
r++;
if (r >= M) {
M = r;
ith = 0;
r = 0;
for (i = 0; i < ETW; i++) {
- if (getbit(buf,i)==OP) {
+ if (bp_getbit(buf,i)==OP) {
r++;
} else {
r--;
}
}
-}
+};
-int bp_construct(bp *b,int n, pb *B, int opt)
+bp * bp_construct(int n, pb *B, int opt)
{
+ bp *b;
int i,j,d;
int m,M,ds;
int ns,nm;
int *md;
int m_ofs;
int r; // # of minimum values
-
-#if 0
- if (SB % D != 0) {
- printf("warning: SB=%d should be a multiple of D=%d\n",SB,D);
- // not necessarily?
- }
- if (SB % RRR != 0) {
- printf("warning: SB=%d should be a multiple of RRR=%d\n",SB,RRR);
- }
-#endif
+ bp_xmalloc(b, 1);
b->B = B;
b->n = n;
b->da_inorder = NULL;
b->da_postorder = NULL;
b->da_dfuds_leaf = NULL;
- mymalloc(b->da,1,0);
- darray_construct(b->da,n,B, opt & OPT_FAST_PREORDER_SELECT);
+ b->da = bp_darray_construct(n,B, opt & OPT_FAST_PREORDER_SELECT);
b->idx_size += b->da->idx_size;
- //Kim: comment this and the following, they polute the printing of the xpath library
- //printf("preorder rank/select table: %d bytes (%1.2f bpc)\n",b->da->idx_size,(double)b->da->idx_size*8/n);
+ //TODO check.
make_matchtbl();
ns = (n+SB-1)/SB;
- mymalloc(sm, ns, 0); b->idx_size += ns * sizeof(*sm);
- mymalloc(sM, ns, 0); b->idx_size += ns * sizeof(*sM);
+
+ bp_xmalloc(sm, ns);
+ b->idx_size += ns * sizeof(*sm);
+
+ bp_xmalloc(sM, ns);
+ b->idx_size += ns * sizeof(*sM);
+
b->sm = sm;
b->sM = sM;
+
if (opt & OPT_DEGREE) {
- mymalloc(sd, ns, 0); b->idx_size += ns * sizeof(*sd);
+ bp_xmalloc(sd, ns);
+ b->idx_size += ns * sizeof(*sd);
b->sd = sd;
- //printf("SB degree table: %d bytes (%1.2f bpc)\n",ns * sizeof(*sd), (double)ns * sizeof(*sd) * 8/n);
}
- //printf("SB table: %d bytes (%1.2f bpc)\n",ns * sizeof(*sm) * 2, (double)ns * sizeof(*sm)*2 * 8/n);
for (i=0; i<n; i++) {
if (i % SB == 0) {
- ds = depth(b,i);
+ ds = bp_depth(b,i);
m = M = ds;
r = 1;
} else {
- d = depth(b,i);
+ d = bp_depth(b,i);
if (d == m) r++;
if (d < m) {
m = d;
if (d > M) M = d;
}
if (i % SB == SB-1 || i==n-1) {
- ds = depth(b,(i/SB)*SB-1);
+ ds = bp_depth(b,(i/SB)*SB-1);
if (m - ds + SB < 0 || m - ds + SB > 255) {
printf("error m=%d ds=%d\n",m,ds);
}
}
}
-#if 0
- printf("sd: ");
- for (i=0;i<n/SB;i++) printf("%d ",sd[i]);
- printf("\n");
-#endif
-
nm = (n+MB-1)/MB;
m_ofs = 1 << blog(nm-1);
b->m_ofs = m_ofs;
- mymalloc(mm, nm + m_ofs, 0); b->idx_size += (nm+m_ofs) * sizeof(*mm);
- mymalloc(mM, nm + m_ofs, 0); b->idx_size += (nm+m_ofs) * sizeof(*mM);
+ bp_xmalloc(mm, nm + m_ofs);
+ b->idx_size += (nm+m_ofs) * sizeof(*mm);
+
+ bp_xmalloc(mM, nm + m_ofs);
+ b->idx_size += (nm+m_ofs) * sizeof(*mM);
b->mm = mm;
b->mM = mM;
+
if (opt & OPT_DEGREE) {
- mymalloc(md, nm + m_ofs, 0); b->idx_size += (nm+m_ofs) * sizeof(*md);
+ bp_xmalloc(md, nm + m_ofs);
+ b->idx_size += (nm+m_ofs) * sizeof(*md);
b->md = md;
- //printf("MB degree table: %d bytes (%1.2f bpc)\n",(nm+m_ofs) * sizeof(*md), (double)(nm+m_ofs) * sizeof(*md) * 8/n);
}
- //printf("MB table: %d bytes (%1.2f bpc)\n",(nm+m_ofs) * sizeof(*mm) * 2, (double)(nm+m_ofs) * sizeof(*mm)*2 * 8/n);
for (i=0; i<n; i++) {
- d = depth(b,i);
+ d = bp_depth(b,i);
if (i % MB == 0) {
m = M = d;
r = 1;
if (opt & OPT_DEGREE) md[m_ofs+ i/MB] = r;
}
}
-
+
for (j=m_ofs-1; j > 0; j--) {
m = 0;
if (j*2 < nm + m_ofs) m = mm[j*2];
}
-#if 0
- printf("md: ");
- for (i=0;i<m_ofs + n/MB;i++) printf("%d ",md[i]);
- printf("\n");
-#endif
-
if (opt & OPT_LEAF) {
- mymalloc(b->da_leaf,1,0);
- darray_pat_construct(b->da_leaf, n, B, 2, 0x2, opt & OPT_FAST_LEAF_SELECT);
- //printf("leaf rank/select table: %d bytes (%1.2f bpc)\n",b->da_leaf->idx_size,(double)b->da_leaf->idx_size*8/n);
- b->idx_size += b->da_leaf->idx_size;
+ b->da_leaf = bp_darray_pat_construct(n, B, 2, 0x2, opt & OPT_FAST_LEAF_SELECT);
+ b->idx_size += b->da_leaf->idx_size;
} else {
b->da_leaf = NULL;
}
if (opt & OPT_INORDER) {
- mymalloc(b->da_inorder,1,0);
- darray_pat_construct(b->da_inorder, n, B, 2, 0x1, opt & OPT_FAST_INORDER_SELECT);
- //printf("inorder rank/select table: %d bytes (%1.2f bpc)\n",b->da_inorder->idx_size,(double)b->da_inorder->idx_size*8/n);
+
+ b->da_inorder = bp_darray_pat_construct(n, B, 2, 0x1, opt & OPT_FAST_INORDER_SELECT);
b->idx_size += b->da_inorder->idx_size;
+
} else {
b->da_inorder = NULL;
}
if (opt & OPT_FAST_POSTORDER_SELECT) {
- mymalloc(b->da_postorder,1,0);
- darray_pat_construct(b->da_postorder, n, B, 1, 0x0, (opt & OPT_FAST_POSTORDER_SELECT) | OPT_NO_RANK);
- //printf("postorder rank/select table: %d bytes (%1.2f bpc)\n",b->da_postorder->idx_size,(double)b->da_postorder->idx_size*8/n);
+
+ b->da_postorder = bp_darray_pat_construct(n, B, 1, 0x0, (opt & OPT_FAST_POSTORDER_SELECT) | OPT_NO_RANK);
b->idx_size += b->da_postorder->idx_size;
+
} else {
b->da_postorder = NULL;
}
if (opt & OPT_DFUDS_LEAF) {
- mymalloc(b->da_dfuds_leaf,1,0);
- darray_pat_construct(b->da_dfuds_leaf, n, B, 2, 0x0, opt & OPT_FAST_DFUDS_LEAF_SELECT);
- //printf("dfuds leaf rank/select table: %d bytes (%1.2f bpc)\n",b->da_dfuds_leaf->idx_size,(double)b->da_dfuds_leaf->idx_size*8/n);
+ b->da_dfuds_leaf = bp_darray_pat_construct( n, B, 2, 0x0, opt & OPT_FAST_DFUDS_LEAF_SELECT);
b->idx_size += b->da_dfuds_leaf->idx_size;
+
} else {
b->da_dfuds_leaf = NULL;
}
- return 0;
+ return b;
}
-// destroyTree: frees the memory of tree.
-void destroyTree(bp *b) {
+void bp_delete(bp *b) {
if (!b) return; // nothing to free
- destroyDarray(b->da); // destroys da data structure
- if (b->da) free(b->da);
-
+ bp_darray_free(b->da); // destroys da data structure
if (b->sm) free(b->sm);
if (b->sM) free(b->sM);
if (b->sd) free(b->sd);
if (b->mm) free(b->mm);
if (b->mM) free(b->mM);
if (b->md) free(b->md);
-
- destroyDarray(b->da_leaf);
- if (b->da_leaf) free(b->da_leaf);
-
- destroyDarray(b->da_inorder);
- if (b->da_inorder) free(b->da_inorder);
-
- destroyDarray(b->da_postorder);
- if (b->da_postorder) free(b->da_postorder);
-
- destroyDarray(b->da_dfuds_leaf);
- if (b->da_dfuds_leaf) free(b->da_dfuds_leaf);
+
+ bp_darray_free(b->da_leaf);
+
+ bp_darray_free(b->da_inorder);
+
+ bp_darray_free(b->da_postorder);
+
+ bp_darray_free(b->da_dfuds_leaf);
+
+ bp_free(b);
}
// loadTree: load parentheses data structure from file
// By Diego Arroyuelo
-void loadTree(bp *b, FILE *fp) {
+bp * loadTree(FILE *fp) {
pb *B;
int n, opt;
exit(1);
}
- mymalloc(B,(n+D-1)/D,0);
-
+ bp_xmalloc(B, (n+D-1)/D);
+
if (fread(B, sizeof(pb), (n+D-1)/D, fp) != ((n+D-1)/D)) {
printf("Error: cannot read parentheses sequence from file\n");
exit(1);
printf("Error: cannot read opt in parentheses from file\n");
exit(1);
}
-
- bp_construct(b, n, B, opt);
-
-}
+ return bp_construct(n, B, opt);
+
+}
-int naive_fwd_excess(bp *b,int s, int rel)
+static int naive_fwd_excess(bp *b,int s, int rel)
{
int i,v,n;
pb *B;
n = b->n; B = b->B;
v = 0;
for (i=s+1; i<n; i++) {
- if (getbit(B,i)==OP) {
+ if (bp_getbit(B,i)==OP) {
v++;
} else {
v--;
return -1;
}
-int naive_bwd_excess(bp *b,int s, int rel)
+static int naive_bwd_excess(bp *b,int s, int rel)
{
int i,v;
pb *B;
B = b->B;
v = 0;
for (i=s; i>=0; i--) {
- if (getbit(B,i)==OP) {
+ if (bp_getbit(B,i)==OP) {
v--;
} else {
v++;
return -2;
}
-int naive_search_SB_l(bp *b, int i, int rel)
+static int naive_search_SB_l(bp *b, int i, int rel)
{
int il,v;
il = (i / SB) * SB;
for (; i>=il; i--) {
- if (getbit(b->B,i)==OP) {
+ if (bp_getbit(b->B,i)==OP) {
rel++;
} else {
rel--;
return -3;
}
-int naive_rmq(bp *b, int s, int t,int opt)
+int bp_naive_rmq(bp *b, int s, int t,int opt)
{
int d,i,dm,im;
if (opt & OPT_RIGHT) {
- d = dm = depth(b,t); im = t;
+ d = dm = bp_depth(b,t); im = t;
i = t-1;
while (i >= s) {
- if (getbit(b->B,i+1)==CP) {
+ if (bp_getbit(b->B,i+1)==CP) {
d++;
if (opt & OPT_MAX) {
if (d > dm) {
i--;
}
} else {
- d = dm = depth(b,s); im = s;
+ d = dm = bp_depth(b,s); im = s;
i = s+1;
while (i <= t) {
- if (getbit(b->B,i)==OP) {
+ if (bp_getbit(b->B,i)==OP) {
d++;
if (opt & OPT_MAX) {
if (d > dm) {
return im;
}
-int root_node(bp *b)
-{
- return 0;
-}
-int rank_open(bp *b, int s)
+int bp_rank_open(bp *b, int s)
{
- return darray_rank(b->da,s);
+ return bp_darray_rank(b->da, s);
}
-int rank_close(bp *b, int s)
+int bp_rank_close(bp *b, int s)
{
- return s+1 - darray_rank(b->da,s);
+ return s+1 - bp_darray_rank(b->da, s);
}
-int select_open(bp *b, int s)
+int bp_select_open(bp *b, int s)
{
if (b->opt & OPT_FAST_PREORDER_SELECT) {
- return darray_select(b->da,s,1);
+ return bp_darray_select(b->da,s,1);
} else {
- return darray_select_bsearch(b->da,s,getpat_preorder);
+ return bp_darray_select_bsearch(b->da, s, getpat_preorder);
}
}
-int select_close(bp *b, int s)
+int bp_select_close(bp *b, int s)
{
if (b->opt & OPT_FAST_POSTORDER_SELECT) {
- return darray_pat_select(b->da_postorder,s,getpat_postorder);
+ return bp_darray_pat_select(b->da_postorder, s, getpat_postorder);
} else {
return postorder_select_bsearch(b,s);
}
}
-///////////////////////////////////////////
-// find_close(bp *b,int s)
-// returns the matching close parenthesis of s
-///////////////////////////////////////////
-int find_close(bp *b,int s)
-{
- return fwd_excess(b,s,-1);
-}
-
-///////////////////////////////////////////
-// find_open(bp *b,int s)
-// returns the matching open parenthesis of s
-///////////////////////////////////////////
-int find_open(bp *b,int s)
-{
- int r;
- r = bwd_excess(b,s,0);
- if (r >= -1) return r+1;
- return -1;
-}
-
-///////////////////////////////////////////
-// parent(bp *b,int s)
-// returns the parent of s
-// -1 if s is the root
-///////////////////////////////////////////
-int parent(bp *b,int s)
-{
- int r;
- r = bwd_excess(b,s,-2);
- if (r >= -1) return r+1;
- return -1;
-}
-
-int enclose(bp *b,int s)
-{
- return parent(b,s);
-}
-
-///////////////////////////////////////////
-// level_ancestor(bp *b,int s,int d)
-// returns the ancestor of s with relative depth d (d < 0)
-// -1 if no such node
-///////////////////////////////////////////
-int level_ancestor(bp *b,int s,int d)
-{
- int r;
- r = bwd_excess(b,s,d-1);
- if (r >= -1) return r+1;
- return -1;
-}
///////////////////////////////////////////
-// lca(bp *b, int s, int t)
-// returns the lowest common ancestor of s and t
-///////////////////////////////////////////
-int lca(bp *b, int s, int t)
-{
- return parent(b,rmq(b,s,t,0)+1);
-}
-
-
-///////////////////////////////////////////
-// preorder_rank(bp *b,int s)
-// returns the preorder (>= 1) of node s (s >= 0)
-///////////////////////////////////////////
-int preorder_rank(bp *b,int s)
-{
- return darray_rank(b->da,s);
-}
-
-///////////////////////////////////////////
-// preorder_select(bp *b,int s)
-// returns the node with preorder s (s >= 1)
-// -1 if no such node
-///////////////////////////////////////////
-int preorder_select(bp *b,int s)
-{
- // no error handling
- if (b->opt & OPT_FAST_PREORDER_SELECT) {
- return darray_select(b->da,s,1);
- } else {
- return darray_select_bsearch(b->da,s,getpat_preorder);
- }
-}
-
-///////////////////////////////////////////
-// postorder_rank(bp *b,int s)
+// bp_postorder_rank(bp *b,int s)
// returns the postorder (>= 1) of node s (s >= 0)
// -1 if s-th bit is not OP
///////////////////////////////////////////
-int postorder_rank(bp *b,int s)
+int bp_postorder_rank(bp *b,int s)
{
int t;
- if (inspect(b,s) == CP) return -1;
- t = find_close(b,s);
+ if (bp_inspect(b,s) == CP) return -1;
+ t = bp_find_close(b,s);
// return t+1 - darray_rank(b->da,t);
- return rank_close(b,t);
+ return bp_rank_close(b,t);
}
-int postorder_select_bsearch(bp *b,int s)
+static int postorder_select_bsearch(bp *b,int s)
{
int l,r,m;
while (l < r) {
m = (l+r)/2;
//printf("m=%d rank=%d s=%d\n",m,m+1 - darray_rank(b->da,m),s);
- if (m+1 - darray_rank(b->da,m) >= s) {
+ if (m+1 - bp_darray_rank(b->da,m) >= s) {
r = m;
} else {
l = m+1;
}
///////////////////////////////////////////
-// postorder_select(bp *b,int s)
+// bp_postorder_select(bp *b,int s)
// returns the position of CP of the node with postorder s (>= 1)
///////////////////////////////////////////
-int postorder_select(bp *b,int s)
+int bp_postorder_select(bp *b,int s)
{
#if 0
if (b->opt & OPT_FAST_POSTORDER_SELECT) {
- return darray_pat_select(b->da_postorder,s,getpat_postorder);
+ return bp_darray_pat_select(b->da_postorder,s,getpat_postorder);
} else {
return postorder_select_bsearch(b->da,s);
}
#else
- return select_close(b,s);
+ return bp_select_close(b,s);
#endif
}
///////////////////////////////////////////
-// leaf_rank(bp *b,int s)
+// bp_leaf_rank(bp *b,int s)
// returns the number of leaves to the left of s
///////////////////////////////////////////
-int leaf_rank(bp *b,int s)
+int bp_leaf_rank(bp *b,int s)
{
if ((b->opt & OPT_LEAF) == 0) {
printf("leaf_rank: error!!! not supported\n");
if (s >= b->n-1) {
s = b->n-2;
}
- return darray_pat_rank(b->da_leaf,s,getpat_leaf);
+ return bp_darray_pat_rank(b->da_leaf,s,getpat_leaf);
}
///////////////////////////////////////////
// leaf_select(bp *b,int s)
// returns the position of s-th leaf
///////////////////////////////////////////
-int leaf_select(bp *b,int s)
+int bp_leaf_select(bp *b,int s)
{
if ((b->opt & OPT_LEAF) == 0) {
printf("leaf_select: error!!! not supported\n");
}
if (s > b->da_leaf->m) return -1;
if (b->opt & OPT_FAST_LEAF_SELECT) {
- return darray_pat_select(b->da_leaf,s,getpat_leaf);
+ return bp_darray_pat_select(b->da_leaf,s,getpat_leaf);
} else {
- return darray_select_bsearch(b->da_leaf,s,getpat_leaf);
+ return bp_darray_select_bsearch(b->da_leaf,s,getpat_leaf);
}
}
///////////////////////////////////////////
-// inorder_rank(bp *b,int s)
+// bp_inorder_rank(bp *b,int s)
// returns the number of ")(" (s >= 0)
///////////////////////////////////////////
-int inorder_rank(bp *b,int s)
+int bp_inorder_rank(bp *b,int s)
{
if ((b->opt & OPT_INORDER) == 0) {
printf("inorder_rank: error!!! not supported\n");
if (s >= b->n-1) {
s = b->n-2;
}
- return darray_pat_rank(b->da_inorder,s,getpat_inorder);
+ return bp_darray_pat_rank(b->da_inorder,s,getpat_inorder);
}
///////////////////////////////////////////
-// inorder_select(bp *b,int s)
+// bp_inorder_select(bp *b,int s)
// returns the s-th position of ")(" (s >= 1)
///////////////////////////////////////////
-int inorder_select(bp *b,int s)
+int bp_inorder_select(bp *b,int s)
{
if ((b->opt & OPT_INORDER) == 0) {
printf("inorder_select: error!!! not supported\n");
return -1;
}
if (b->opt & OPT_FAST_INORDER_SELECT) {
- return darray_pat_select(b->da_inorder,s,getpat_inorder);
+ return bp_darray_pat_select(b->da_inorder,s,getpat_inorder);
} else {
- return darray_select_bsearch(b->da_inorder,s,getpat_inorder);
+ return bp_darray_select_bsearch(b->da_inorder,s,getpat_inorder);
}
}
///////////////////////////////////////////
-// leftmost_leaf(bp *b, int s)
+// bp_leftmost_leaf(bp *b, int s)
///////////////////////////////////////////
-int leftmost_leaf(bp *b, int s)
+int bp_leftmost_leaf(bp *b, int s)
{
if ((b->opt & OPT_LEAF) == 0) {
printf("leftmost_leaf: error!!! not supported\n");
return -1;
}
- return leaf_select(b,leaf_rank(b,s)+1);
+ return bp_leaf_select(b, bp_leaf_rank(b,s)+1);
}
///////////////////////////////////////////
// rightmost_leaf(bp *b, int s)
///////////////////////////////////////////
-int rightmost_leaf(bp *b, int s)
+int bp_rightmost_leaf(bp *b, int s)
{
int t;
if ((b->opt & OPT_LEAF) == 0) {
printf("leftmost_leaf: error!!! not supported\n");
return -1;
}
- t = find_close(b,s);
- return leaf_select(b,leaf_rank(b,t));
-}
-
-
-
-///////////////////////////////////////////
-// inspect(bp *b, int s)
-// returns OP (==1) or CP (==0) at s-th bit (0 <= s < n)
-///////////////////////////////////////////
-int inspect(bp *b, int s)
-{
- if (s < 0 || s >= b->n) {
- printf("inspect: error s=%d is out of [%d,%d]\n",s,0,b->n-1);
- }
- return getbit(b->B,s);
+ t = bp_find_close(b,s);
+ return bp_leaf_select(b, bp_leaf_rank(b,t));
}
-int isleaf(bp *b, int s)
-{
- if (inspect(b,s) != OP) {
- printf("isleaf: error!!! B[%d] = OP\n",s);
- }
- if (inspect(b,s+1) == CP) return 1;
- else return 0;
-}
-
-///////////////////////////////////////////
-// subtree_size(bp *b, int s)
-// returns the number of nodes in the subtree of s
///////////////////////////////////////////
-int subtree_size(bp *b, int s)
-{
- return (find_close(b,s) - s + 1) / 2;
-}
-
-///////////////////////////////////////////
-// first_child(bp *b, int s)
-// returns the first child
-// -1 if s is a leaf
-///////////////////////////////////////////
-int first_child(bp *b, int s)
-{
- if (inspect(b,s+1) == CP) return -1;
- return s+1;
-}
-
-///////////////////////////////////////////
-// next_sibling(bp *b,int s)
-// returns the next sibling of parent(s)
-// -1 if s is the last child
-//////////////////////////////////////////
-int next_sibling(bp *b, int s)
-{
- int t;
- t = find_close(b,s)+1;
- if (t >= b->n) {
- printf("next_sibling: error s=%d t=%d\n",s,t);
- }
- if (inspect(b,t) == CP) return -1;
- return t;
-}
-
-///////////////////////////////////////////
-// prev_sibling(bp *b,int s)
+// bp_prev_sibling(bp *b,int s)
// returns the previous sibling of parent(s)
// -1 if s is the first child
//////////////////////////////////////////
-int prev_sibling(bp *b, int s)
+int bp_prev_sibling(bp *b, int s)
{
int t;
- if (s < 0) {
- printf("prev_sibling: error s=%d\n",s);
- }
if (s == 0) return -1;
- if (inspect(b,s-1) == OP) return -1;
- t = find_open(b,s-1);
+ if (bp_inspect(b,s-1) == OP) return -1;
+ t = bp_find_open(b,s-1);
return t;
}
///////////////////////////////////////////
-// deepest_node(bp *b,int s)
+// bp_deepest_node(bp *b,int s)
// returns the first node with the largest depth in the subtree of s
///////////////////////////////////////////
-int deepest_node(bp *b,int s)
+int bp_deepest_node(bp *b,int s)
{
int t,m;
- t = find_close(b,s);
- m = rmq(b,s,t, OPT_MAX);
+ t = bp_find_close(b,s);
+ m = bp_rmq(b,s,t, OPT_MAX);
return m;
}
///////////////////////////////////////////
-// subtree_height(bp *b,int s)
+// bp_subtree_height(bp *b,int s)
// returns the height of the subtree of s
// 0 if s is a leaf
///////////////////////////////////////////
-int subtree_height(bp *b,int s)
+int bp_subtree_height(bp *b,int s)
{
int t;
- t = deepest_node(b,s);
- return depth(b,t) - depth(b,s);
+ t = bp_deepest_node(b,s);
+ return bp_depth(b,t) - bp_depth(b,s);
}
-int naive_degree(bp *b, int s)
+int bp_naive_degree(bp *b, int s)
{
int t,d;
d = 0;
- t = first_child(b,s);
+ t = bp_first_child(b,s);
while (t >= 0) {
d++;
- t = next_sibling(b,t);
+ t = bp_next_sibling(b,t);
}
return d;
}
///////////////////////////////////////////
-// degree(bp *b, int s)
+// bp_degree(bp *b, int s)
// returns the number of children of s
// 0 if s is a leaf
///////////////////////////////////////////
-int degree(bp *b, int s)
+int bp_degree(bp *b, int s)
{
if (b->opt & OPT_DEGREE) {
- return fast_degree(b,s,b->n,0);
+ return bp_fast_degree(b,s,b->n,0);
} else {
- return naive_degree(b,s);
+ return bp_naive_degree(b,s);
}
}
-int naive_child(bp *b, int s, int d)
+int bp_naive_child(bp *b, int s, int d)
{
int t,i;
- t = first_child(b,s);
+ t = bp_first_child(b,s);
for (i = 1; i < d; i++) {
if (t == -1) break;
- t = next_sibling(b,t);
+ t = bp_next_sibling(b,t);
}
return t;
}
///////////////////////////////////////////
-// child(bp *b, int s, int d)
+// bp_child(bp *b, int s, int d)
// returns the d-th child of s (1 <= d <= degree(s))
// -1 if no such node
///////////////////////////////////////////
-int child(bp *b, int s, int d)
+int bp_child(bp *b, int s, int d)
{
int r;
if (b->opt & OPT_DEGREE) {
//return find_open(b,fast_degree(b,s,b->n,d));
- if (d==1) return first_child(b,s);
- r = fast_degree(b,s,b->n,d-1)+1;
- if (inspect(b,r) == CP) return -1;
+ if (d==1) return bp_first_child(b,s);
+ r = bp_fast_degree(b,s,b->n,d-1)+1;
+ if (bp_inspect(b,r) == CP) return -1;
return r;
} else {
- return naive_child(b,s,d);
+ return bp_naive_child(b,s,d);
}
}
-int naive_child_rank(bp *b, int t)
+int bp_naive_child_rank(bp *b, int t)
{
int v,d;
d = 0;
while (t != -1) {
d++;
- t = prev_sibling(b,t);
+ t = bp_prev_sibling(b,t);
}
return d;
}
///////////////////////////////////////////
-// child_rank(bp *b, int t)
+// bp_child_rank(bp *b, int t)
// returns d if t is the d-th child of the parent of t (d >= 1)
// 1 if t is the root
///////////////////////////////////////////
-int child_rank(bp *b, int t)
+int bp_child_rank(bp *b, int t)
{
int r;
- if (t == root_node(b)) return 1;
+ if (t == bp_root_node(b)) return 1;
if (b->opt & OPT_DEGREE) {
- r = parent(b,t);
- return fast_degree(b,r,t,0)+1;
+ r = bp_parent(b,t);
+ return bp_fast_degree(b,r,t,0)+1;
} else {
- return naive_child_rank(b,t);
+ return bp_naive_child_rank(b,t);
}
}
-
-///////////////////////////////////////////
-// is_ancestor(bp *b, int s, int t)
-// returns 1 if s is an ancestor of t
-// 0 otherwise
-///////////////////////////////////////////
-int is_ancestor(bp *b, int s, int t)
-{
- int v;
- v = find_close(b,s);
- if (s <= t && t <= v) return 1;
- return 0;
-}
-
///////////////////////////////////////////
// distance(bp *b, int s, int t)
// returns the length of the shortest path from s to t in the tree
///////////////////////////////////////////
-int distance(bp *b, int s, int t)
+int bp_distance(bp *b, int s, int t)
{
int v,d;
- v = lca(b,s,t);
- d = depth(b,v);
- return (depth(b,s) - d) + (depth(b,t) - d);
+ v = bp_lca(b,s,t);
+ d = bp_depth(b,v);
+ return (bp_depth(b,s) - d) + (bp_depth(b,t) - d);
}
///////////////////////////////////////////
// level_next(bp *b, int d)
///////////////////////////////////////////
-int level_next(bp *b,int s)
+int bp_level_next(bp *b,int s)
{
int t;
- t = fwd_excess(b,s,0);
+ t = bp_fwd_excess(b,s,0);
return t;
}
///////////////////////////////////////////
// level_prev(bp *b, int d)
///////////////////////////////////////////
-int level_prev(bp *b,int s)
+int bp_level_prev(bp *b,int s)
{
int t;
- t = bwd_excess(b,s,0);
+ t = bp_bwd_excess(b,s,0);
return t;
}
///////////////////////////////////////////
-// level_leftmost(bp *b, int d)
+// bp_level_leftmost(bp *b, int d)
///////////////////////////////////////////
-int level_leftmost(bp *b, int d)
+int bp_level_leftmost(bp *b, int d)
{
int t;
if (d < 1) return -1;
if (d == 1) return 0;
- t = fwd_excess(b,0,d);
+ t = bp_fwd_excess(b,0,d);
return t;
}
///////////////////////////////////////////
-// level_rigthmost(bp *b, int d)
+// bp_level_rigthmost(bp *b, int d)
///////////////////////////////////////////
int level_rigthmost(bp *b, int d)
{
int t;
if (d < 1) return -1;
if (d == 1) return 0;
- t = bwd_excess(b,0,d-1);
- return find_open(b,t);
+ t = bp_bwd_excess(b,0,d-1);
+ return bp_find_open(b,t);
}
///////////////////////////////////////////
// leaf_size(bp *b, int s)
///////////////////////////////////////////
-int leaf_size(bp *b, int s)
+int bp_leaf_size(bp *b, int s)
{
int t;
if ((b->opt & OPT_LEAF) == 0) {
printf("leaf_size: error!!! not supported\n");
return -1;
}
- t = find_close(b,s);
- return leaf_rank(b,t) - leaf_rank(b,s);
+ t = bp_find_close(b,s);
+ return bp_leaf_rank(b,t) - bp_leaf_rank(b,s);
}
#ifndef BP_H_
#define BP_H_
+
#ifdef __cplusplus
extern "C" {
#endif
#include <stdio.h>
#include <stdlib.h>
-#include "darray.h"
+#include "bp-darray.h"
+#include "bp-utils.h"
+
#define OP 1
#define CP 0
int opt;
} bp;
-int bp_construct(bp *b,int n, pb *B, int opt);
-void printbp(bp *b, int s, int t);
-
-
-int rank_open(bp *b, int s);
-int rank_close(bp *b, int s);
-int select_open(bp *b, int s);
-int select_close(bp *b, int s);
-
-
-int root_node(bp *b);
-int find_close(bp *b,int s);
-int find_open(bp *b,int s);
-int enclose(bp *b,int s);
-int parent(bp *b,int s);
-int level_ancestor(bp *b,int s,int d);
-int depth(bp *b, int s);
-int preorder_rank(bp *b,int s);
-int postorder_rank(bp *b,int s);
-int inspect(bp *b, int s);
-int isleaf(bp *b, int s);
-int rmq(bp *b, int s, int t, int opt);
-int subtree_size(bp *b, int s);
-int first_child(bp *b, int s);
-int next_sibling(bp *b, int s);
-int prev_sibling(bp *b, int s);
-int deepest_node(bp *b,int s);
-int subtree_height(bp *b,int s);
-int is_ancestor(bp *b, int s, int t);
-int distance(bp *b, int s, int t);
-int level_lefthmost(bp *b, int d);
-int level_rigthmost(bp *b, int d);
-int degree(bp *b,int s);
+pb getpat_preorder(pb *b);
+pb getpat_leaf(pb *b);
+pb getpat_inorder(pb *b);
+pb getpat_postorder(pb *b);
+
+bp * bp_construct(int n, pb *B, int opt);
+void bp_print(bp *b, int s, int t);
+
+int bp_fwd_excess(bp *b,int s, int rel);
+int bp_bwd_excess(bp *b,int s, int rel);
+int bp_rmq(bp *b, int s, int t, int opt);
+int bp_depth(bp *b, int s);
+
+int bp_rank_open(bp *b, int s);
+int bp_rank_close(bp *b, int s);
+int bp_select_open(bp *b, int s);
+int bp_select_close(bp *b, int s);
+
+
+static inline int bp_root_node(bp *b)
+{
+ return 0;
+}
+
+
+///////////////////////////////////////////
+// find_close(bp *b,int s)
+// returns the matching close parenthesis of s
+///////////////////////////////////////////
+static inline int bp_find_close(bp *b,int s)
+{
+ return bp_fwd_excess(b, s, -1);
+}
+
+///////////////////////////////////////////
+// bp_find_open(bp *b,int s)
+// returns the matching open parenthesis of s
+///////////////////////////////////////////
+static inline int bp_find_open(bp *b,int s)
+{
+ int r = bp_bwd_excess(b, s, 0);
+ return (r >= -1) ? r+1 : -1;
+}
+
+
+///////////////////////////////////////////
+// bp_parent(bp *b,int s)
+// returns the parent of s
+// -1 if s is the root
+///////////////////////////////////////////
+static inline int bp_parent(bp *b, int s)
+{
+ int r = bp_bwd_excess(b,s,-2);
+ return (r >= -1) ? r+1 : -1;
+}
+
+///////////////////////////////////////////
+// bp_parent_close(bp *b,int s)
+// returns the closing parenthesis of the parent
+// of s
+// -1 if s is the root
+///////////////////////////////////////////
+
+static inline int bp_parent_close(bp *b, int s)
+{
+ return bp_fwd_excess(b,s,-2);
+}
+
+
+static inline int bp_enclose(bp *b, int s)
+{
+ return bp_parent(b, s);;
+}
+
+
+///////////////////////////////////////////
+// bp_level_ancestor(bp *b,int s,int d)
+// returns the ancestor of s with relative depth d (d < 0)
+// -1 if no such node
+///////////////////////////////////////////
+static inline int bp_level_ancestor(bp *b,int s,int d)
+{
+ int r = bp_bwd_excess(b,s,d-1);
+ return (r >= -1) ? r+1 :-1;
+}
+
+///////////////////////////////////////////
+// bp_lca(bp *b, int s, int t)
+// returns the lowest common ancestor of s and t
+///////////////////////////////////////////
+static inline int bp_lca(bp *b, int s, int t)
+{
+ return bp_parent(b, bp_rmq(b,s,t,0)+1);
+}
+
+///////////////////////////////////////////
+// preorder_rank(bp *b,int s)
+// returns the preorder (>= 1) of node s (s >= 0)
+///////////////////////////////////////////
+
+static inline int bp_preorder_rank(bp *b,int s)
+{
+ return bp_darray_rank(b->da,s);
+}
+
+///////////////////////////////////////////
+// preorder_select(bp *b,int s)
+// returns the node with preorder s (s >= 1)
+// -1 if no such node
+///////////////////////////////////////////
+static inline int bp_preorder_select(bp *b,int s)
+{
+ // no error handling
+ if (b->opt & OPT_FAST_PREORDER_SELECT) {
+ return bp_darray_select(b->da,s,1);
+ } else {
+ return bp_darray_select_bsearch(b->da,s, getpat_preorder);
+ }
+}
+
+
+int bp_postorder_rank(bp *b,int s);
+
+///////////////////////////////////////////
+// inspect(bp *b, int s)
+// returns OP (==1) or CP (==0) at s-th bit (0 <= s < n)
+///////////////////////////////////////////
+static inline int bp_inspect(bp *b, int s)
+{
+ return bp_getbit(b->B,s);
+}
+
+static inline int bp_isleaf(bp *b, int s)
+{
+ return (bp_inspect(b, s+1) == CP);
+}
+
+///////////////////////////////////////////
+// bp_subtree_size(bp *b, int s)
+// returns the number of nodes in the subtree of s
+///////////////////////////////////////////
+static inline int bp_subtree_size(bp *b, int s)
+{
+ return (bp_find_close(b, s) - s + 1) / 2;
+}
+
+///////////////////////////////////////////
+// first_child(bp *b, int s)
+// returns the first child
+// -1 if s is a leaf
+///////////////////////////////////////////
+
+static inline int bp_first_child(bp *b, int s)
+{
+ return (bp_inspect(b, s+1) == CP) ? -1 : s+1;
+}
+
+
+///////////////////////////////////////////
+// next_sibling(bp *b,int s)
+// returns the next sibling of parent(s)
+// -1 if s is the last child
+//////////////////////////////////////////
+static inline int bp_next_sibling(bp *b, int s)
+{
+ int t;
+ t = bp_find_close(b,s) + 1;
+ return (bp_inspect(b, t) == CP) ? -1 : t;
+
+}
+
+int bp_prev_sibling(bp *b, int s);
+int bp_deepest_node(bp *b,int s);
+int bp_subtree_height(bp *b,int s);
+///////////////////////////////////////////
+// bp_is_ancestor(bp *b, int s, int t)
+// returns 1 if s is an ancestor of t
+// 0 otherwise
+///////////////////////////////////////////
+static inline int bp_is_ancestor(bp *b, int s, int t)
+{
+ int v;
+ if (s == 0) return 1;
+ v = bp_find_close(b, s);
+ return (s <= t && t <= v);
+}
+
+int bp_distance(bp *b, int s, int t);
+int bp_level_lefthmost(bp *b, int d);
+int bp_level_rigthmost(bp *b, int d);
+int bp_degree(bp *b,int s);
// not efficient
-int naive_degree(bp *b, int s);
-int naive_child(bp *b, int s, int d);
-int naive_child_rank(bp *b, int t);
-int naive_rmq(bp *b, int s, int t,int opt);
-int postorder_select(bp *b,int s);
+int bp_naive_degree(bp *b, int s);
+int bp_naive_child(bp *b, int s, int d);
+int bp_naive_child_rank(bp *b, int t);
+int bp_naive_rmq(bp *b, int s, int t,int opt);
+int bp_postorder_select(bp *b,int s);
// using preorder select index
int preorder_select(bp *b,int s);
// using leaf index
-int leaf_rank(bp *b,int s);
-int leaf_size(bp *b, int s);
-int leftmost_leaf(bp *b, int s);
-int rightmost_leaf(bp *b, int s);
+int bp_leaf_rank(bp *b,int s);
+int bp_leaf_size(bp *b, int s);
+int bp_leftmost_leaf(bp *b, int s);
+int bp_rightmost_leaf(bp *b, int s);
// using leaf select index
-int leaf_select(bp *b,int s);
+int bp_leaf_select(bp *b,int s);
// using inorder index
-int inorder_rank(bp *b,int s);
+int bp_inorder_rank(bp *b,int s);
// using inorder select index
-int inorder_select(bp *b,int s);
+int bp_inorder_select(bp *b,int s);
// using degree index
-int fast_degree(bp *b,int s, int t, int ith);
-int child_rank(bp *b, int t);
-int child(bp *b, int s, int d);
+int bp_fast_degree(bp *b,int s, int t, int ith);
+int bp_child_rank(bp *b, int t);
+int bp_child(bp *b, int s, int d);
// new functions for persistence purposes, added by Diego Arroyuelo
void saveTree(bp *b, FILE *fp);
-void loadTree(bp *b, FILE *fp);
-void destroyTree(bp *b);
+bp * loadTree(FILE *fp);
+void bp_delete(bp *b);
static inline int blog(int x)
return l;
}
-pb getpat_preorder(pb *b);
-pb getpat_leaf(pb *b);
-pb getpat_inorder(pb *b);
-pb getpat_postorder(pb *b);
-int fwd_excess(bp *b,int s, int rel);
-int bwd_excess(bp *b,int s, int rel);
extern int *matchtbl,*parenttbl;
void make_naivetbl(pb *B,int n);
-extern int popCount[1<<ETW];
+
extern int fwdtbl[(2*ETW+1)*(1<<ETW)];
extern int bwdtbl[(2*ETW+1)*(1<<ETW)];
extern int mintbl_li[1<<ETW], mintbl_lv[1<<ETW];
+++ /dev/null
-#include "bp.h"
-#include "utils.h"
-
-
-#ifndef min
-#define min(x,y) ((x)<(y)?(x):(y))
-#endif
-
-#ifndef max
-#define max(x,y) ((x)>(y)?(x):(y))
-#endif
-
-#define NOTFOUND -2
-#define CONTINUE -3
-#define END -4
-#define FOUND -5
-
-#define SBid(i) ((i)>>logSB)
-#define SBfirst(i) ((i) & (~(SB-1)))
-#define SBlast(i) ((i) | (SB-1))
-
-#define MBid(i) ((i)>>logMB)
-#define MBfirst(i) ((i) & (~(MB-1)))
-#define MBlast(i) ((i) | (MB-1))
-
-pb getpat_preorder(pb *b)
-{
- return *b;
-}
-
-pb getpat_postorder(pb *b)
-{
- return ~(*b);
-}
-
-pb getpat_leaf(pb *b)
-{
- pb w1,w2,w;
- w1 = b[0];
- w2 = (w1 << 1) + (b[1] >> (D-1));
- w = w1 & (~w2);
- return w;
-}
-
-pb getpat_inorder(pb *b)
-{
- pb w1,w2,w;
- w1 = b[0];
- w2 = (w1 << 1) + (b[1] >> (D-1));
- w = (~w1) & w2;
- return w;
-}
-
-pb getpat_dfuds_leaf(pb *b)
-{
- pb w1,w2,w;
- w1 = b[0];
- w2 = (w1 << 1) + (b[1] >> (D-1));
- w = (~w1) & (~w2);
- return w;
-}
-
-
-
-///////////////////////////////////////////
-// depth(bp *b, int s)
-// returns the depth of s
-// The root node has depth 1
-///////////////////////////////////////////
-int depth(bp *b, int s)
-{
- int d;
- if (s < 0) return 0;
- d = 2 * darray_rank(b->da,s) - (s+1);
-#if 0
- if (d != naive_depth(b,s)) {
- d = naive_depth(b,s);
- darray_rank(b->da,s);
- }
- //printf("depth(%d)=%d\n",s,d);
-#endif
- return d;
-}
-
-int fast_depth(bp *b, int s)
-{
- int d;
- darray *da;
- int r,j;
-
- s++;
- if ((s & (RRR-1)) != 0) {
- //printf("fast_depth:warning s=%d\n",s);
- return depth(b,s);
- }
- da = b->da;
- //d = 2 * (da->rl[s>>logR] + da->rm[s>>logRR] + da->rs[s>>logRRR]) - s;
-
- r = da->rl[s>>logR] + da->rm[s>>logRR];
- j = (s>>logRRR) & (RR/RRR-1);
- while (j > 0) {
- r += da->rs[((s>>logRR)<<(logRR-logRRR))+j-1];
- j--;
- }
- d = 2 * r - s;
-
- return d;
-}
-
-int search_SB_r(bp *b, int i, int rel)
-{
- int j,r,n,il;
- pb *p,x,w;
-
- n = b->n;
- il = min((SBid(i) + 1) << logSB,n);
- p = &b->B[i>>logD];
- while (i<il) {
- x = *p++;
- j = i & (D-1);
- x <<= j;
- j = D-j;
- while (j>0) {
- w = (x >> (D-ETW)) & ((1<<ETW)-1);
- if (rel >= -ETW && rel <= ETW) {
- r = fwdtbl[((rel+ETW)<<ETW)+w];
- if (r<ETW && r<j) {
- if (i+r >= n) return NOTFOUND;
- return i+r;
- }
- }
- r = min(j,ETW);
- rel -= 2*popcount(w)-r;
- x <<= r;
- i += r;
- j -= r;
- }
- }
- return CONTINUE;
-}
-
-int search_MB_r(bp *b, int i, int td)
-{
- int il,d;
- int m,M,n;
- pb *B;
-
- B = b->B;
- n = b->n;
-
- il = min((MBid(i) + 1) << logMB,n);
- for ( ; i < il; i+=SB) {
-#if (SB % RRR != 0)
- d = depth(b,i-1);
-#else
- d = fast_depth(b,i-1);
-#endif
- m = d + b->sm[SBid(i)] - SB;
- M = d + b->sM[SBid(i)] - 1;
- if (m <= td && td <= M) {
- return search_SB_r(b,i,td-d);
- }
- }
- if (i >= n) return NOTFOUND;
- return CONTINUE;
-}
-
-///////////////////////////////////////////
-// fwd_excess(bp *b,int s, int rel)
-// find the leftmost value depth(s)+rel to the right of s (exclusive)
-///////////////////////////////////////////
-int fwd_excess(bp *b,int s, int rel)
-{
- int i,n;
- int d,td;
- int m,M;
- int m_ofs;
- pb *B;
- n = b->n; B = b->B;
-
- i = s+1;
-
- d = search_SB_r(b,i,rel);
- if (d >= NOTFOUND) return d;
-
- i = min((SBid(i) + 1) << logSB, n);
- td = depth(b,s) + rel;
- d = search_MB_r(b,i,td);
- if (d >= NOTFOUND) return d;
-
- m_ofs = b->m_ofs;
-
- i = MBid(s) + m_ofs;
-
- while (i > 0) {
- if ((i&1) == 0) {
- i++;
- m = b->mm[i];
- M = b->mM[i];
-
- if (m <= td && td <= M) {
- while (i < m_ofs) {
- i <<= 1;
- m = b->mm[i];
- M = b->mM[i];
- if (!(m <= td && td <= M)) i++;
- }
- i -= m_ofs;
- i <<= logMB;
-
- return search_MB_r(b,i,td);
- };
-
- }
- i >>= 1;
- }
- return NOTFOUND;
-}
-
-
-int degree_SB_slow(bp *b, int i, int t, int rel, int *ans, int ith)
-{
- int j,r,n,il;
- pb *p,x,w,w2;
- int d, deg, v;
-
- n = t;
- il = min((SBid(i) + 1) << logSB,n);
- d = deg = 0;
-
- while (i < il) {
- if (getbit(b->B,i)==OP) {
- d++;
- } else {
- d--;
- }
- if (d < rel) { // reached the end
- if (ith > 0) {
- return NOTFOUND;
- } else {
- *ans = deg;
- return END;
- }
- }
- if (d == rel) { // found the same depth
- deg++;
- if (deg == ith) {
- *ans = i;
- return FOUND;
- }
- }
- i++;
- }
- *ans = deg;
- return CONTINUE;
-}
-
-int degree_SB(bp *b, int i, int t, int rel, int *ans, int ith)
-{
- int j,r,n,il;
- pb *p,x,w,w2;
- int d, deg, v;
- int degtmp,itmp;
- int ith2,d2;
-
- n = t;
- il = min((SBid(i) + 1) << logSB,n);
- d = deg = 0;
-
- p = &b->B[i>>logD];
- while (i < il) {
- x = *p++;
- j = i & (D-1);
- x <<= j;
- j = min(D-j,il-i);
- while (j>0) {
- w = (x >> (D-ETW)) & ((1<<ETW)-1);
- w2 = 0;
- if (j < ETW || il-i < ETW) {
- r = max(ETW-j,ETW-(il-i));
- w2 = (1<<r)-1;
- }
- v = minmaxtbl_v[0][w | w2];
- if (d + v < rel) {
- if (ith > 0) {
-#if 0
- for (r = 0; r < ETW; r++) {
- if (w & 0x80) {
- d++;
- } else {
- d--;
- if (d < rel) break;
- if (d == rel) {
- ith--;
- if (ith == 0) {
- *ans = i + r;
- return FOUND;
- }
- }
- }
- w <<= 1;
- }
- return NOTFOUND;
-#else
- r = childtbl2[rel-d+ETW][ith-1][w];
- if (r >= 0) {
- *ans = i + r;
- return FOUND;
- }
- return NOTFOUND;
-#endif
- }
- r = ETW-1-minmaxtbl_i[0][w | w2];
- w2 = (1<<r)-1;
- deg += degtbl2[((rel-d+ETW)<<ETW) + (w & (~w2))];
- *ans = deg;
- return END;
- }
- if (d + v == rel) {
- r = degtbl[w | w2];
- deg += r;
- if (ith > 0) {
- if (ith <= r) {
- *ans = i + childtbl[((ith-1)<<ETW) + (w | w2)];
- return FOUND;
- }
- ith -= r;
- }
- }
-
- r = min(j,ETW);
- d += 2*popcount(w)-r;
- x <<= r;
- i += r;
- j -= r;
- }
- }
-
- *ans = deg;
- return CONTINUE;
-}
-
-int degree_MB(bp *b, int i, int t, int td, int *ans,int ith)
-{
- int il,d;
- int m,M,n,r;
- pb *B;
- int deg,degtmp;
-
- d = 0;
- B = b->B;
- n = t;
-
- il = min((MBid(i) + 1) << logMB,n);
- deg = 0;
- for ( ; i+SB-1 < il; i+=SB) {
-#if (SB % RRR != 0)
- d = depth(b,i-1);
-#else
- d = fast_depth(b,i-1);
-#endif
- m = d + b->sm[SBid(i)] - SB;
- if (m < td) {
- r = degree_SB(b,i,n,td-d,°tmp,ith);
- if (ith > 0) {
- if (r == NOTFOUND) return NOTFOUND;
- *ans = degtmp;
- return FOUND;
- } else {
- *ans = deg + degtmp;
- return END;
- }
- }
- if (m == td) {
- if (ith > 0) {
- if (ith <= b->sd[SBid(i)]) break;
- ith -= b->sd[SBid(i)];
- }
- deg += b->sd[SBid(i)];
- }
- }
- if (i < il) {
-#if (SB % RRR != 0)
- d = depth(b,i-1);
-#else
- d = fast_depth(b,i-1);
-#endif
- r = degree_SB(b,i,n,td-d,°tmp,ith);
- if (ith > 0) {
- if (r == NOTFOUND) return NOTFOUND;
- if (r == FOUND) {
- *ans = degtmp;
- return FOUND;
- }
- } else {
- deg += degtmp;
- }
- }
- *ans = deg;
- if (i >= n) return END;
- return CONTINUE;
-}
-
-static int partition_range(int s,int t)
-{
- int i,j,h;
-
- printf("partition [%d,%d] => ",s,t);
- h = 1;
- while (s <= t) {
- if (s & h) {
- if (s+h-1 <= t) {
- printf("[%d,%d] ",s,s+h-1);
- s += h;
- }
- } else {
- if (s+h > t) break;
- }
- h <<= 1;
- }
- while (h > 0) {
- if (s+h-1 <= t) {
- printf("[%d,%d] ",s,s+h-1);
- s += h;
- }
- h >>= 1;
- }
- printf("\n");
-}
-
-
-
-
-///////////////////////////////////////////
-// fast_degree(bp *b,int s, int t, int ith)
-// returns the number of children of s, to the left of t
-// returns the position of (ith)-th child of s if (ith > 0)
-///////////////////////////////////////////
-int fast_degree(bp *b,int s, int t, int ith)
-{
- int i,j,n;
- int d,td;
- int m,M;
- int m_ofs;
- pb *B;
- int deg,degtmp;
- int sm,tm,ss,h;
-
- n = t;
- B = b->B;
-
- deg = 0;
-
- i = s+1;
- if (i != SBfirst(i)) {
- d = degree_SB(b,i,n,0,°tmp,ith);
- if (ith > 0) {
- if (d == NOTFOUND) return -1;
- if (d == FOUND) return degtmp;
- ith -= degtmp;
- }
- if (d == END) return degtmp;
- deg += degtmp;
- }
-
- td = depth(b,s);
-
- i = SBid(i+SB-1) << logSB;
-
- if (i != MBfirst(i)) {
- d = degree_MB(b,i,n,td,°tmp,ith);
- if (ith > 0) {
- if (d == NOTFOUND) return -1;
- if (d == FOUND) return degtmp;
- ith -= degtmp;
- deg += degtmp;
- } else {
- deg += degtmp;
- if (d == END) {
- return deg;
- }
- }
- }
-
-#if 0
- // sequential search
-
- sm = MBid(i+MB-1);
- tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
-
- m_ofs = b->m_ofs;
- sm += m_ofs; tm += m_ofs;
- for (i=sm; i<=tm; i++) {
- if (b->mm[i] < td) {
- break;
- }
- if (b->mm[i] == td) {
- if (ith > 0) {
- if (ith <= b->md[i]) break;
- ith -= b->md[i];
- }
- deg += b->md[i];
- }
- }
- ss = i - m_ofs;
-#else
- sm = MBid(i+MB-1);
- tm = MBid((n-1)+1)-1; // the rightmost MB fully contained in [0,n-1]
-
- m_ofs = b->m_ofs;
- sm += m_ofs; tm += m_ofs;
- ss = sm;
-
- //partition_range(sm,tm);
-
- //printf("partition [%d,%d] => ",sm,tm);
- h = 1;
- while (sm <= tm) {
- if (sm & h) {
- if (sm+h-1 <= tm) {
- //printf("[%d,%d] ",sm,sm+h-1);
- j = sm / h;
- if (b->mm[j] < td) {
- h >>= 1;
- break;
- }
- if (b->mm[j] == td) {
- if (ith > 0) {
- if (ith <= b->md[j]) {
- h >>= 1;
- break;
- }
- ith -= b->md[j];
- }
- deg += b->md[j];
- }
- sm += h;
- }
- } else {
- if (sm+h > tm) break;
- }
- h <<= 1;
- }
- while (h > 0) {
- if (sm+h-1 <= tm) {
- //printf("[%d,%d] ",sm,sm+h-1);
- j = sm / h;
- if (ith > 0) {
- if (b->mm[j] >= td) {
- if (b->mm[j] == td) {
- if (ith > b->md[j]) {
- ith -= b->md[j];
- sm += h;
- } else {
- deg += b->md[j];
- }
- } else {
- sm += h;
- }
- }
- } else {
- if (b->mm[j] >= td) {
- if (b->mm[j] == td) {
- deg += b->md[j];
- }
- sm += h;
- }
- }
- }
- h >>= 1;
- }
- //printf("\n");
- ss = sm;
-
- ss -= m_ofs;
-
-#endif
-
- ss <<= logMB;
-
- d = degree_MB(b,ss,n,td,°tmp,ith);
- if (ith > 0) {
- if (d == NOTFOUND) return -1;
- if (d == FOUND) return degtmp;
- }
- deg += degtmp;
- if (d == END) return deg;
- return deg;
-
- // unexpected (bug)
- printf("degree: ???\n");
- return -99;
-
-}
-
-int search_SB_l(bp *b, int i, int rel)
-{
- int j,r,il;
- pb *p,x,w;
-
- il = SBfirst(i);
-
- p = &b->B[i>>logD];
- while (i>=il) {
- x = *p--;
- j = (i & (D-1))+1;
- x >>= D-j;
- while (j>0) {
- w = x & ((1<<ETW)-1);
- if (rel >= -ETW && rel <= ETW) {
- r = bwdtbl[((rel+ETW)<<ETW)+w];
- if (r<ETW && r<j) {
- if (i-r < 0) return NOTFOUND;
- return i-r-1;
- }
- }
- r = min(j,ETW);
- rel += 2*popcount(w)-r;
- x >>= r;
- i -= r;
- j -= r;
- }
-
- }
- if (i < 0) return NOTFOUND;
- return CONTINUE;
-}
-
-int search_MB_l(bp *b, int i, int td)
-{
- int il,d;
- int m,M;
- pb *B;
-
-#if 0
- if (i % SB != SB-1) {
- printf("search_MB_l:error!!! i=%d SB=%d\n",i,i%SB);
- }
-#endif
- B = b->B;
-
- il = MBfirst(i);
- for ( ; i >= il; i-=SB) {
-#if (SB % RRR != 0)
- d = depth(b,i-SB);
-#else
- d = fast_depth(b,i-SB);
-#endif
- m = d + b->sm[SBid(i)] - SB;
- M = d + b->sM[SBid(i)] - 1;
- if (m <= td && td <= M) {
-#if (SB % RRR != 0)
- d = depth(b,i);
-#else
- d = fast_depth(b,i);
-#endif
- if (d == td) return i;
- return search_SB_l(b,i,td-d);
- }
- }
- return CONTINUE;
-}
-
-///////////////////////////////////////////
-// bwd_excess(bp *b,int s, int rel)
-// find the rightmost value depth(s)+rel to the left of s (exclusive)
-///////////////////////////////////////////
-int bwd_excess(bp *b,int s, int rel)
-{
- int i,n;
- int d,td;
- int m,M;
- int m_ofs;
- pb *B;
- n = b->n; B = b->B;
-
- i = s;
- d = search_SB_l(b,i,rel);
- if (d >= NOTFOUND) return d;
-
- i = SBfirst(i) -1;
-
- td = depth(b,s) + rel;
-
- d = search_MB_l(b,i,td);
- if (d >= NOTFOUND) return d;
-
- m_ofs = b->m_ofs;
- i = (s>>logMB) + m_ofs;
- while (i > 0) {
- if ((i&1) == 1) {
- i--;
- m = b->mm[i];
- M = b->mM[i];
- if (m <= td && td <= M) break;
- }
- i >>= 1;
- }
- if (i == 0) {
- if (td == 0) return -1;
- else return NOTFOUND;
- }
-
- while (i < m_ofs) {
- i = i*2 + 1;
- m = b->mm[i];
- M = b->mM[i];
- if (!(m <= td && td <= M)) i--;
- }
- i -= m_ofs;
- i = ((i+1)<<logMB)-1;
-
- d = search_MB_l(b,i,td);
- if (d >= NOTFOUND) return d;
-
- // unexpected (bug)
- printf("bwd_excess: ???\n");
- return -99;
-
-}
-
-int rmq_SB(bp *b, int s, int t, int opt, int *dm)
-{
- int i,d;
- int is,ds;
- pb *p,x,w,w2;
- int j,v,r;
- int lr;
- int op;
-
- lr = 0; if (opt & OPT_RIGHT) lr = 1;
- op = opt & (OPT_RIGHT | OPT_MAX);
-
- is = s; ds = d = 0;
- i = s+1;
-
-#if SB >= ETW
- p = &b->B[i>>logD];
- while (i <= t) {
- x = *p++;
- j = i & (D-1);
- x <<= j;
- j = min(D-j,t-i+1);
- while (j>0) {
- w = (x >> (D-ETW)) & ((1<<ETW)-1);
- w2 = 0;
- if (j < ETW || t-i < ETW-1) {
- r = max(ETW-j,ETW-1-(t-i));
- w2 = (1<<r)-1;
- }
-
- if (op & OPT_MAX) {
- v = minmaxtbl_v[op][w & (~w2)];
- if (d + v + lr > ds) {
- ds = d + v;
- is = i + minmaxtbl_i[op][w & (~w2)];
- }
- } else {
- v = minmaxtbl_v[op][w | w2];
- if (d + v < ds +lr) {
- ds = d + v;
- is = i + minmaxtbl_i[op][w | w2];
- }
- }
-
- r = min(j,ETW);
- d += 2*popcount(w)-r;
- x <<= r;
- i += r;
- j -= r;
- }
- }
-#else
- while (i <= t) {
- if (getbit(b->B,i)==OP) {
- d++;
- if (op & OPT_MAX) {
- if (d + lr > ds) {
- ds = d; is = i;
- }
- }
- } else {
- d--;
- if (!(op & OPT_MAX)) {
- if (d < ds + lr) {
- ds = d; is = i;
- }
- }
- }
- i++;
- }
-#endif
- *dm = ds;
- return is;
-}
-
-int rmq_MB(bp *b, int s, int t, int opt, int *dm)
-{
- int i,d,m;
- int mi,md;
- int lr;
-
- lr = 0; if (opt & OPT_RIGHT) lr = 1;
-
- md = *dm; mi = -1;
- for (i = s; i <= t; i++) {
-#if (SB % RRR != 0)
- d = depth(b,(i<<logSB)-1);
-#else
- d = fast_depth(b,(i<<logSB)-1);
-#endif
- if (opt & OPT_MAX) {
- m = d + b->sM[i] - 1;
- if (m + lr > md) {
- md = m; mi = i;
- }
- } else {
- m = d + b->sm[i] - SB;
- if (m < md + lr) {
- md = m; mi = i;
- }
- }
- }
- *dm = md;
- return mi;
-}
-
-
-
-
-///////////////////////////////////////////
-// rmq(bp *b, int s, int t, int opt)
-// returns the index of leftmost/rightmost minimum/maximum value
-// in the range [s,t] (inclusive)
-// returns the leftmost minimum if opt == 0
-// the rightmost one if (opt & OPT_RIGHT)
-// the maximum if (opt & OPT_MAX)
-///////////////////////////////////////////
-int rmq(bp *b, int s, int t, int opt)
-{
- int ss, ts; // SB index of s and t
- int sm, tm; // MB index of s and t
- int ds; // current min value
- int is; // current min index
- int ys; // type of current min index
- // 0: is is the index of min
- // 1: is is the SB index
- // 2: is is the MB index
- int m_ofs;
- int i,j,il,d,n;
- int dm;
- int lr;
-
- lr = 0; if (opt & OPT_RIGHT) lr = 1;
-
- n = b->n;
- if (s < 0 || t >= n || s > t) {
- printf("rmq: error s=%d t=%d n=%d\n",s,t,n);
- return -1;
- }
- if (s == t) return s;
-
-
- ////////////////////////////////////////////////////////////
- // search the SB of s
- ////////////////////////////////////////////////////////////
-
- il = min(SBlast(s),t);
- is = rmq_SB(b,s,il,opt,&dm);
- if (il == t) { // scan reached the end of the range
- return is;
- }
- ds = depth(b,s) + dm; ys = 0;
-
- ////////////////////////////////////////////////////////////
- // search the MB of s
- ////////////////////////////////////////////////////////////
-
- ss = SBid(s) + 1;
- il = min(SBid(MBlast(s)),SBid(t)-1);
- dm = ds;
- j = rmq_MB(b,ss,il,opt,&dm);
- //if (dm < ds + lr) {
- if (j >= 0) {
- ds = dm; is = j; ys = 1;
- }
-
- ////////////////////////////////////////////////////////////
- // search the tree
- ////////////////////////////////////////////////////////////
-
- sm = MBid(s) + 1;
- tm = MBid(t) - 1;
-
-#if 0
- // sequential search
- m_ofs = b->m_ofs;
- sm += m_ofs; tm += m_ofs;
- for (i=sm; i<=tm; i++) {
- if (opt & OPT_MAX) {
- if (b->mM[i] + lr > ds) {
- ds = b->mM[i]; is = i - m_ofs; ys = 2;
- }
- } else {
- if (b->mm[i] < ds + lr) {
- ds = b->mm[i]; is = i - m_ofs; ys = 2;
- }
- }
- }
-
-#else
- if (sm <= tm) {
- int h;
- h = blog(sm ^ tm);
-
- m_ofs = b->m_ofs;
- sm += m_ofs; tm += m_ofs;
-
- if (opt & OPT_MAX) {
- if (b->mM[sm] + lr > ds) {
- ds = b->mM[sm]; is = sm; ys = 2;
- }
- for (i=0; i<=h-2; i++) {
- j = sm>>i;
- if ((j&1) == 0) {
- if (b->mM[j+1] + lr > ds) {
- ds = b->mM[j+1]; is = j+1; ys = 2;
- }
- }
- }
- for (i=h-2; i>=0; i--) {
- j = tm>>i;
- if ((j&1)==1) {
- if (b->mM[j-1] + lr > ds) {
- ds = b->mM[j-1]; is = j-1; ys = 2;
- }
- }
- }
- if (b->mM[tm] + lr > ds) {
- ds = b->mM[tm]; is = tm; ys = 2;
- }
- if (ys == 2) {
- while (is < m_ofs) {
- is <<= 1;
- if (b->mM[is+1] + lr > b->mM[is]) is++;
- }
- is -= m_ofs;
- }
- } else { // MIN
- if (b->mm[sm] < ds + lr) {
- ds = b->mm[sm]; is = sm; ys = 2;
- }
- for (i=0; i<=h-2; i++) {
- j = sm>>i;
- if ((j&1) == 0) {
- if (b->mm[j+1] < ds + lr) {
- ds = b->mm[j+1]; is = j+1; ys = 2;
- }
- }
- }
- for (i=h-2; i>=0; i--) {
- j = tm>>i;
- if ((j&1)==1) {
- if (b->mm[j-1] < ds + lr) {
- ds = b->mm[j-1]; is = j-1; ys = 2;
- }
- }
- }
- if (b->mm[tm] < ds + lr) {
- ds = b->mm[tm]; is = tm; ys = 2;
- }
- if (ys == 2) {
- while (is < m_ofs) {
- is <<= 1;
- if (b->mm[is+1] < b->mm[is] + lr) is++;
- }
- is -= m_ofs;
- }
- }
- }
-
-#endif
-
- ////////////////////////////////////////////////////////////
- // search the MB of t
- ////////////////////////////////////////////////////////////
-
-
- ts = max(SBid(MBfirst(t)),SBid(s)+1);
- il = SBid(SBfirst(t)-1);
- dm = ds;
- j = rmq_MB(b,ts,il,opt,&dm);
- //if (dm < ds + lr) {
- if (j >= 0) {
- ds = dm; is = j; ys = 1;
- }
-
- ////////////////////////////////////////////////////////////
- // search the SB of t
- ////////////////////////////////////////////////////////////
-
- i = SBfirst(t);
- j = rmq_SB(b,i,t,opt,&dm);
- d = depth(b,i) + dm;
- if (opt & OPT_MAX) {
- if (d + lr > ds) {
- ds = d; is = j; ys = 0;
- }
- } else {
- if (d < ds + lr) {
- ds = d; is = j; ys = 0;
- }
- }
-
- ////////////////////////////////////////////////////////////
- // search the rest
- ////////////////////////////////////////////////////////////
-
- if (ys == 2) {
- ss = SBid(is << logMB);
- il = SBid(MBlast(is << logMB));
- if (opt & OPT_MAX) {
- dm = -n-1;
- } else {
- dm = n+1;
- }
- j = rmq_MB(b,ss,il,opt,&dm);
- ds = dm; is = j; ys = 1;
- }
-
- if (ys == 1) {
- ss = is << logSB;
- il = SBlast(is << logSB);
- is = rmq_SB(b,ss,il,opt,&dm);
- //ds = depth(b,ss) + dm; ys = 0;
- }
-
- return is;
-}
-
+++ /dev/null
-#include <stdio.h>\r
-#include <stdlib.h>\r
-#include "darray.h"\r
-#include "utils.h"\r
-\r
-//typedef unsigned char byte;\r
-//typedef unsigned short word;\r
-//typedef unsigned int dword;\r
-\r
-//typedef dword pb;\r
-//#define logD 5\r
-\r
-#define PBS (sizeof(pb)*8)\r
-#define D (1<<logD)\r
-#define logM 5\r
-#define M (1<<logM)\r
-#define logP 8\r
-#define P (1<<logP)\r
-#define logLL 16 // size of word\r
-#define LL (1<<logLL)\r
-#define logLLL 7\r
-//#define logLLL 5\r
-#define LLL (1<<logLLL)\r
-//#define logL 10\r
-//#define logL (logLL-3)\r
-#define logL (logLL-1-5)\r
-#define L (1<<logL)\r
-\r
-#ifndef min\r
- #define min(x,y) ((x)<(y)?(x):(y))\r
-#endif\r
-\r
-#define mymalloc(p,n,f) {p = (__typeof__(p)) malloc((n)*sizeof(*p)); if ((p)==NULL) {printf("not enough memory\n"); exit(1);}}\r
-\r
-int setbit(pb *B, int i,int x)\r
-{\r
- int j,l;\r
-\r
- j = i / D;\r
- l = i % D;\r
- if (x==0) B[j] &= (~(1<<(D-1-l)));\r
- else if (x==1) B[j] |= (1<<(D-1-l));\r
- else {\r
- printf("error setbit x=%d\n",x);\r
- exit(1);\r
- }\r
- return x;\r
-}\r
-\r
-int setbit_zero(pb *B, int i)\r
-{\r
- int j,l;\r
- j = i >> logD;\r
- l = i & (D-1);\r
- B[j] &= (~(1<<(D-1-l)));\r
- return 0;\r
-}\r
-\r
-int setbit_one(pb *B, int i)\r
-{\r
- int j,l;\r
- j = i >> logD;\r
- l = i & (D-1);\r
- B[j] |= (1<<(D-1-l));\r
- return 1;\r
-\r
-}\r
-\r
-\r
-\r
-int setbits(pb *B, int i, int d, int x)\r
-{\r
- int j;\r
-\r
- for (j=0; j<d; j++) {\r
- setbit(B,i+j,(x>>(d-j-1))&1);\r
- }\r
- return x;\r
-}\r
-\r
-int getbit(pb *B, int i)\r
-{\r
- int j,l;\r
-\r
- //j = i / D;\r
- //l = i % D;\r
- j = i >> logD;\r
- l = i & (D-1);\r
- return (B[j] >> (D-1-l)) & 1;\r
-}\r
-\r
-dword getbits(pb *B, int i, int d)\r
-{\r
- dword j,x;\r
-\r
- x = 0;\r
- for (j=0; j<d; j++) {\r
- x <<= 1;\r
- x += getbit(B,i+j);\r
- }\r
- return x;\r
-}\r
-\r
-int getpattern(pb *B, int i, int k, pb pat)\r
-{\r
- int j;\r
- int x;\r
- x = 1;\r
- for (j=0; j<k; j++) {\r
- x &= getbit(B,i+j) ^ (~(pat>>(k-1-j)));\r
- }\r
- //printf("getpattern(%d) = %d\n",i,x);\r
- return x;\r
-}\r
-\r
-\r
-static const unsigned int popCount[] = {\r
-0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,\r
-1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,\r
-1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,\r
-2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
-1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,\r
-2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
-2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
-3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,\r
-1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,\r
-2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
-2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
-3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,\r
-2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
-3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,\r
-3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,\r
-4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8\r
-};\r
-\r
-static int selecttbl[8*256];\r
-\r
-void make_selecttbl(void)\r
-{\r
- int i,x,r;\r
- pb buf[1];\r
-\r
- for (x = 0; x < 256; x++) {\r
- setbits(buf,0,8,x);\r
- for (r=0; r<8; r++) selecttbl[(r<<8)+x] = -1;\r
- r = 0;\r
- for (i=0; i<8; i++) {\r
- if (getbit(buf,i)) {\r
- selecttbl[(r<<8)+x] = i;\r
- r++;\r
- }\r
- }\r
- }\r
-}\r
-\r
-\r
-int darray_construct(darray *da, int n, pb *buf, int opt)\r
-{\r
- int i,j,k,m;\r
- int nl;\r
- int p,pp;\r
- int il,is,ml,ms;\r
- int r,m2;\r
- dword *s;\r
- int p1,p2,p3,p4,s1,s2,s3,s4;\r
-\r
- da->idx_size = 0;\r
-\r
- make_selecttbl();\r
-\r
-\r
- if (L/LLL == 0) {\r
- printf("ERROR: L=%d LLL=%d\n",L,LLL);\r
- exit(1);\r
- }\r
-\r
- m = 0;\r
- for (i=0; i<n; i++) m += getbit(buf,i);\r
- da->n = n;\r
- da->m = m;\r
- //printf("n=%d m=%d\n",n,m);\r
-\r
- da->buf = buf;\r
-\r
- if (opt & (~OPT_NO_RANK)) { // construct select table\r
-#if 0\r
- mymalloc(s,m,0);\r
- m = 0;\r
- for (i=0; i<n; i++) {\r
- if (getbit(buf,i)) {\r
- m++;\r
- s[m-1] = i;\r
- }\r
- }\r
-#endif \r
- nl = (m-1) / L + 1;\r
- mymalloc(da->lp,nl+1,1); da->idx_size += (nl+1)*sizeof(*da->lp);\r
- mymalloc(da->p,nl+1,1); da->idx_size += (nl+1)*sizeof(*da->p);\r
-#if 0\r
- printf("lp table: %d bytes (%1.2f bpc)\n",(nl+1)*sizeof(*da->lp), (double)(nl+1)*sizeof(*da->lp) * 8/n);\r
- printf("p table: %d bytes (%1.2f bpc)\n",(nl+1)*sizeof(*da->p), (double)(nl+1)*sizeof(*da->p) * 8/n);\r
-#endif \r
-\r
- for (r = 0; r < 2; r++) {\r
- s1 = s2 = s3 = s4 = 0;\r
- p1 = p2 = p3 = p4 = -1;\r
- \r
- ml = ms = 0;\r
- for (il = 0; il < nl; il++) {\r
- //pp = s[il*L];\r
- while (s1 <= il*L) {\r
- if (getbit(buf,p1+1)) s1++;\r
- p1++;\r
- }\r
- pp = p1;\r
- da->lp[il] = pp;\r
- i = min((il+1)*L-1,m-1);\r
- //p = s[i];\r
- while (s2 <= i) {\r
- if (getbit(buf,p2+1)) s2++;\r
- p2++;\r
- }\r
- p = p2;\r
- //printf("%d ",p-pp);\r
- if (p - pp >= LL) {\r
- if (r == 1) {\r
- for (is = 0; is < L; is++) {\r
- if (il*L+is >= m) break;\r
- //da->sl[ml*L+is] = s[il*L+is];\r
- while (s3 <= il*L+is) {\r
- if (getbit(buf,p3+1)) s3++;\r
- p3++;\r
- }\r
- da->sl[ml*L+is] = p3;\r
- }\r
- }\r
- da->p[il] = -(ml+1);\r
- ml++;\r
- } else {\r
- if (r == 1) {\r
- for (is = 0; is < L/LLL; is++) {\r
- if (il*L+is*LLL >= m) break;\r
- while (s4 <= il*L+is*LLL) {\r
- if (getbit(buf,p4+1)) s4++;\r
- p4++;\r
- }\r
- //da->ss[ms*(L/LLL)+is] = s[il*L+is*LLL] - pp;\r
- da->ss[ms*(L/LLL)+is] = p4 - pp;\r
- }\r
- }\r
- da->p[il] = ms;\r
- ms++;\r
- }\r
- }\r
- if (r == 0) {\r
- mymalloc(da->sl,ml*L+1,1); da->idx_size += (ml*L+1)*sizeof(*da->sl);\r
- mymalloc(da->ss,ms*(L/LLL)+1,1); da->idx_size += (ms*(L/LLL)+1)*sizeof(*da->ss);\r
-#if 0\r
- printf("sl table: %d bytes (%1.2f bpc)\n",(ml*L+1)*sizeof(*da->sl), (double)(ml*L+1)*sizeof(*da->sl) * 8/n);\r
- printf("ss table: %d bytes (%1.2f bpc)\n",(ms*(L/LLL)+1)*sizeof(*da->ss), (double)(ms*(L/LLL)+1)*sizeof(*da->ss) * 8/n);\r
-#endif\r
- }\r
- }\r
- //free(s);\r
- } else { // no select table\r
- da->lp = NULL;\r
- da->p = da->sl = NULL;\r
- da->ss = NULL;\r
- }\r
-\r
- // construct rank table\r
-\r
- if ((opt & OPT_NO_RANK) == 0) {\r
- mymalloc(da->rl,n/R1+2,1); da->idx_size += (n/R1+2)*sizeof(*da->rl);\r
- mymalloc(da->rm,n/RR+2,1); da->idx_size += (n/RR+2)*sizeof(*da->rm);\r
- mymalloc(da->rs,n/RRR+2,1); da->idx_size += (n/RRR+2)*sizeof(*da->rs);\r
-#if 0\r
- printf("rl table: %d bytes (%1.2f bpc)\n",(n/R1+2)*sizeof(*da->rl), (double)(n/R1+2)*sizeof(*da->rl) * 8/n);\r
- printf("rm table: %d bytes (%1.2f bpc)\n",(n/RR+2)*sizeof(*da->rm), (double)(n/RR+2)*sizeof(*da->rm) * 8/n);\r
- printf("rs table: %d bytes (%1.2f bpc)\n",(n/RRR+2)*sizeof(*da->rs), (double)(n/RRR+2)*sizeof(*da->rs) * 8/n);\r
-#endif\r
- r = 0;\r
- for (k=0; k<=n; k+=R1) {\r
- da->rl[k/R1] = r;\r
- m2 = 0;\r
- for (i=0; i<R1; i+=RR) {\r
- if (k+i <= n) da->rm[(k+i)/RR] = m2;\r
- m = 0;\r
- for (j=0; j<RR; j++) {\r
- if (k+i+j < n && getbit(buf,k+i+j)==1) m++;\r
- if (j % RRR == RRR-1) {\r
- if (k+i+j <= n) da->rs[(k+i+j)/RRR] = m;\r
- m2 += m;\r
- m = 0;\r
- }\r
- }\r
- if (m != 0) {\r
- printf("???\n");\r
- }\r
- //m2 += m;\r
- }\r
- r += m2;\r
- }\r
- }\r
-\r
- return 0;\r
-}\r
-\r
-// destroyDarray: frees the memory of darray\r
-// Added by Diego Arroyuelo\r
-void destroyDarray(darray *da) {\r
- if (!da) return; // nothing to free\r
- \r
- if (da->buf) free(da->buf);\r
- if (da->lp) free(da->lp);\r
- if (da->sl) free(da->sl);\r
- if (da->ss) free(da->ss);\r
- if (da->p) free(da->p);\r
- if (da->rl) free(da->rl);\r
- if (da->rm) free(da->rm);\r
- if (da->rs) free(da->rs);\r
-}\r
-\r
-int darray_rank0(darray *da, int i)\r
-{\r
- int r,j;\r
- pb *p;\r
-\r
-#if (RRR == D*2)\r
- r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];\r
- p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
- j = i & (RRR-1);\r
- if (j < D) r += popcount(*p >> (D-1-j));\r
- else r += popcount(*p) + popcount(p[1] >> (D-1-(j-D)));\r
-#else\r
-\r
- j = i & (RRR-1);\r
- if (j < RRR/2) {\r
- r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];\r
- p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
- while (j >= D) {\r
- r += popcount(*p++);\r
- j -= D;\r
- }\r
- r += popcount(*p >> (D-1-j));\r
- } else {\r
- j = RRR-1 - (i & (RRR-1));\r
- i += j+1;\r
- r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];\r
- p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
- while (j >= D) {\r
- r -= popcount(*--p);\r
- j -= D;\r
- }\r
- if (j > 0) r -= popcount(*--p << (D-j));\r
- }\r
-\r
-#endif\r
-\r
- return r;\r
-}\r
-\r
-int darray_rank(darray *da, int i)\r
-{\r
- int r,j,i_rr, i_rrr;\r
- pb *p;\r
- i_rr = i >> logRR;\r
- i_rrr = i >> logRRR;\r
- r = da->rl[i>>logR] + da->rm[i_rr];\r
-\r
- j = (i_rrr) & (RR/RRR-1);\r
- while (j > 0) {\r
- r += da->rs[((i_rr)<<(logRR-logRRR))+j-1];\r
- j--;\r
- }\r
-\r
- p = da->buf + ((i_rrr)<<(logRRR-logD));\r
- j = i & (RRR-1);\r
- while (j >= D) {\r
- r += popcount(*p++);\r
- j -= D;\r
- }\r
- r += popcount(*p >> (D-1-j));\r
-\r
- return r;\r
-}\r
-\r
-int darray_select_bsearch(darray *da, int i, pb (*getpat)(pb *))\r
-{\r
- int j;\r
- int l,r,m,n;\r
- int s;\r
- int t,x,rr;\r
- pb *p,w;\r
-\r
- // for debug\r
- //s = darray_select(da,i,1);\r
- //\r
- //printf("select(%d)=%d\n",i,s);\r
-\r
-\r
-\r
- if (i == 0) return -1;\r
-\r
- if (i > da->m) {\r
- return -1;\r
- }\r
-\r
- i--;\r
-\r
-\r
-\r
- n = da->m;\r
-\r
- t = i;\r
-\r
- l = 0; r = (n-1)>>logR;\r
- // find the smallest index x s.t. rl[x] >= t\r
- while (l < r) {\r
- m = (l+r)/2;\r
- //printf("m=%d rl[m+1]=%d t=%d\n",m,da->rl[m+1],t);\r
- if (da->rl[m+1] > t) { // m+1 is out of range\r
- r = m; // new r = m >= l\r
- } else {\r
- l = m+1; // new l = m+1 <= r\r
- }\r
- }\r
- x = l;\r
- t -= da->rl[x];\r
-\r
- x <<= logR;\r
-\r
- l = x >> logRR; r = (min(x+R1-1,n))>>logRR;\r
- while (l < r) {\r
- m = (l+r)/2;\r
- //printf("m=%d rm[m+1]=%d t=%d\n",m,da->rm[m+1],t);\r
- if (da->rm[m+1] > t) { // m+1 is out of range\r
- r = m;\r
- } else {\r
- l = m+1; // new l = m+1 <= r\r
- }\r
- }\r
- x = l;\r
- t -= da->rm[x];\r
-\r
- x <<= logRR;\r
-\r
-#if 0\r
- l = x >> logRRR; r = (min(x+RR-1,n))>>logRRR;\r
- while (l < r) {\r
- m = (l+r)/2;\r
- //printf("m=%d rs[m+1]=%d t=%d\n",m,da->rs[m+1],t);\r
- if (da->rs[m+1] > t) { // m+1 is out of range\r
- r = m;\r
- } else {\r
- l = m+1; // new l = m+1 <= r\r
- }\r
- }\r
- x = l;\r
- t -= da->rs[x];\r
-#else\r
- l = x >> logRRR;\r
- while (t > da->rs[l]) {\r
- t -= da->rs[l];\r
- l++;\r
- }\r
- x = l;\r
-#endif\r
-\r
- x <<= logRRR;\r
-\r
- p = &da->buf[x >> logD];\r
- while (1) {\r
- m = popcount(getpat(p));\r
- if (m > t) break;\r
- t -= m;\r
- x += D;\r
- p++;\r
- }\r
-\r
- w = getpat(p);\r
- while (1) {\r
- rr = popCount[w >> (D-8)];\r
- if (rr > t) break;\r
- t -= rr;\r
- x += 8;\r
- w <<= 8;\r
- }\r
- x += selecttbl[((t-0)<<8)+(w>>(D-8))];\r
-\r
-#if 0\r
- if (x != s) {\r
- printf("error x=%d s=%d\n",x,s);\r
- }\r
-#endif\r
- return x;\r
-}\r
-\r
-int darray_pat_rank(darray *da, int i, pb (*getpat)(pb *))\r
-{\r
- int r,j;\r
- pb *p;\r
-\r
- r = da->rl[i>>logR] + da->rm[i>>logRR];\r
- j = (i>>logRRR) & (RR/RRR-1);\r
- while (j > 0) {\r
- r += da->rs[((i>>logRR)<<(logRR-logRRR))+j-1];\r
- j--;\r
- }\r
-\r
- p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
- j = i & (RRR-1);\r
- while (j >= D) {\r
- r += popcount(getpat(p));\r
- p++;\r
- j -= D;\r
- }\r
- r += popcount(getpat(p) >> (D-1-j));\r
-\r
- return r;\r
-}\r
-\r
-\r
-int darray_select(darray *da, int i,int f)\r
-{\r
- int p,r;\r
- int il;\r
- int rr;\r
- pb x;\r
- pb *q;\r
-\r
- if (i == 0) return -1;\r
-\r
- if (i > da->m) {\r
- return -1;\r
- //printf("ERROR: m=%d i=%d\n",da->m,i);\r
- //exit(1);\r
- }\r
-\r
- i--;\r
-\r
- il = da->p[i>>logL];\r
- if (il < 0) {\r
- il = -il-1;\r
- p = da->sl[(il<<logL)+(i & (L-1))];\r
- } else {\r
- p = da->lp[i>>logL];\r
- p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];\r
- r = i - (i & (LLL-1));\r
-\r
- q = &(da->buf[p>>logD]);\r
-\r
- if (f == 1) {\r
- rr = p & (D-1);\r
- r -= popcount(*q >> (D-1-rr));\r
- p = p - rr;\r
- \r
- while (1) {\r
- rr = popcount(*q);\r
- if (r + rr >= i) break;\r
- r += rr;\r
- p += D;\r
- q++;\r
- }\r
- \r
- x = *q;\r
- while (1) {\r
- //rr = popcount(x >> (D-8));\r
- //rr = popcount(x >> (D-8));\r
- rr = popcount8(x >> (D-8));\r
- if (r + rr >= i) break;\r
- r += rr;\r
- p += 8;\r
- x <<= 8;\r
- }\r
- p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];\r
- } else {\r
- rr = p & (D-1);\r
- r -= popcount((~(*q)) >> (D-1-rr));\r
- p = p - rr;\r
- \r
- while (1) {\r
- rr = popcount(~(*q));\r
- if (r + rr >= i) break;\r
- r += rr;\r
- p += D;\r
- q++;\r
- }\r
- \r
- x = ~(*q);\r
-\r
- while (1) {\r
- //rr = popcount(x >> (D-8));\r
- //rr = popCount[x >> (D-8)];\r
- rr = popcount8(x >> (D-8));\r
- if (r + rr >= i) break;\r
- r += rr;\r
- p += 8;\r
- x <<= 8;\r
- }\r
- p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];\r
- }\r
- }\r
- return p;\r
-}\r
-\r
-int darray_pat_select(darray *da, int i, pb (*getpat)(pb *))\r
-{\r
- int p,r;\r
- int il;\r
- int rr;\r
- pb x;\r
- pb *q;\r
-\r
- if (i == 0) return -1;\r
-\r
- if (i > da->m) {\r
- return -1;\r
- //printf("ERROR: m=%d i=%d\n",da->m,i);\r
- //exit(1);\r
- }\r
-\r
- i--;\r
-\r
- il = da->p[i>>logL];\r
- if (il < 0) {\r
- il = -il-1;\r
- p = da->sl[(il<<logL)+(i & (L-1))];\r
- } else {\r
- p = da->lp[i>>logL];\r
- p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];\r
- r = i - (i & (LLL-1));\r
-\r
- q = &(da->buf[p>>logD]);\r
-\r
- rr = p & (D-1);\r
- r -= popcount(getpat(q) >> (D-1-rr));\r
- p = p - rr;\r
- \r
- while (1) {\r
- rr = popcount(getpat(q));\r
- if (r + rr >= i) break;\r
- r += rr;\r
- p += D;\r
- q++;\r
- }\r
- \r
- x = getpat(q);\r
- while (1) {\r
- //rr = popcount(x >> (D-8));\r
- //rr = popCount[x >> (D-8)];\r
- rr = popcount8(x >> (D-8));\r
- if (r + rr >= i) break;\r
- r += rr;\r
- p += 8;\r
- x <<= 8;\r
- }\r
- p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];\r
- }\r
- return p;\r
-}\r
-\r
-int darray_pat_construct(darray *da, int n, pb *buf, int k, pb pat, int opt)\r
-{\r
- int i;\r
- pb *b;\r
- mymalloc(b,(n+D-1)/D,0);\r
-\r
- for (i=0; i<n-k+1; i++) {\r
- setbit(b,i,getpattern(buf,i,k,pat));\r
- }\r
- for (i=n-k+1; i<n; i++) {\r
- setbit(b,i,0);\r
- }\r
-\r
- darray_construct(da,n,b,opt);\r
- da->buf = buf;\r
-\r
- free(b);\r
- \r
- return 0;\r
-}\r
+++ /dev/null
-#ifndef DARRAY_H_
-#define DARRAY_H_
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-typedef unsigned char byte;
-typedef unsigned short word;
-typedef unsigned int dword;
-
-#define OPT_NO_RANK (1<<30)
-
-
-typedef dword pb;
-
-#define logD 5
-#define D (1<<logD)
-
-#define logR 16
-#define R1 (1<<logR)
-#define logRR 10
-//#define logRR 8
-#define RR (1<<logRR)
-#define logRRR 7
-#define RRR (1<<logRRR)
-
-typedef struct {
- int n,m;
- int size;
- pb *buf;
- dword *lp;
- dword *sl;
- word *ss;
- dword *p;
-
- dword *rl;
- word *rm;
- byte *rs;
-
- int idx_size;
-} darray;
-
-int setbit(pb *B, int i,int x);
-int setbits(pb *B, int i, int d, int x);
-int getbit(pb *B, int i);
-dword getbits(pb *B, int i, int d);
-//unsigned int popcount(pb x);
-
-int darray_construct(darray *da, int n, pb *buf,int opt);
-int darray_select(darray *da, int i,int f);
-int darray_rank(darray *da, int i);
-int darray_pat_construct(darray *da, int n, pb *buf, int k, pb pat, int opt);
-int darray_pat_select(darray *da, int i, pb (*getpat)(pb *));
-int darray_pat_rank(darray *da, int i, pb (*getpat)(pb *));
-
-int darray_select_bsearch(darray *da, int i, pb (*getpat)(pb *));
-
-// Added by Diego Arroyuelo
-void destroyDarray(darray *da);
-
-#ifdef __cplusplus
-}
-#endif
-
-
-#endif
+++ /dev/null
-#ifndef UTILS_H_
-#define UTILS_H_
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-
-#ifdef HAS_NATIVE_POPCOUNT
-static inline unsigned int popcount(pb n){
- asm ("popcnt %1, %0" : "=r" (n) : "0" (n));
- return n;
-}
-
-static inline unsigned int popcount8(pb n) {
- return popcount(n & 0xff);
-}
-
-#else
-
-static unsigned int popcount8(pb x)
-{
- dword r;
- r = x;
- r = ((r & 0xaa)>>1) + (r & 0x55);
- r = ((r & 0xcc)>>2) + (r & 0x33);
- r = ((r>>4) + r) & 0x0f;
- return r;
-}
-
-static inline unsigned int
-popcount(pb x)
-{
- uint m1 = 0x55555555;
- uint m2 = 0xc30c30c3;
- x -= (x >> 1) & m1;
- x = (x & m2) + ((x >> 2) & m2) + ((x >> 4) & m2);
- x += x >> 6;
- return (x + (x >> 12) + (x >> 24)) & 0x3f;
-}
-
-#endif
-
-#ifdef __cplusplus
- }
-#endif
-
-
-#endif
projects / SXSI / libbp.git / commitdiff