--- /dev/null
+POPCOUNT=$(shell grep -q popcnt /proc/cpuinfo && echo 1)
+
+ifeq ($(POPCOUNT), 1)
+ POPCOUNT_FLAG=-DHAS_NATIVE_POPCOUNT
+else
+ POPCOUNT_FLAG=
+endif
+
+ifeq ($(VERBOSE), true)
+ HIDE=
+else
+ HIDE=@
+endif
+
+ifeq ($(DEBUG), true)
+ OPT_FLAGS=-O0 -g $(POPCOUNT_FLAG) -fno-PIC -static
+else
+ OPT_FLAGS=-O4 $(POPCOUNT_FLAG) -fno-PIC -static
+endif
+
+
+
+INC_FLAGS=-I.
+CFLAGS= $(INC_FLAGS) $(OPT_FLAGS)
+CXXFLAGS= $(INC_FLAGS) $(OPT_FLAGS)
+CC=gcc
+
+
+OBJECTS_BP=bp.o darray.o bpcore.o
+LIB_BP=libbpa.a
+
+all: depend $(LIB_BP)
+
+$(LIB_BP): $(OBJECTS_BP)
+ @echo [BP]
+ $(HIDE) ar rcs libbp.a $(OBJECTS_BP)
+
+%o: %c
+ @echo [C] $@
+ $(HIDE) $(CC) -c $(CFLAGS) $< -o $@
+
+%o: %cpp
+ @echo [C++] $@
+ $(HIDE) $(CC) -c $(CXXFLAGS) $< -o $@
+
+depend:
+ @echo [DEPEND]
+ $(HIDE) (gcc -MM *.c) > $@
+
+clean:
+ @echo [CLEAN]
+ $(HIDE) rm -f *.[ao] depend
+
+-include depend
+
--- /dev/null
+#include "bp.h"
+
+#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; })
+
+//#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);
+
+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) {
+ d++;
+ } else {
+ d--;
+ }
+ }
+ return d;
+}
+
+void printbp(bp *b, int s, int t)
+{
+ int i,c,d;
+ d = 0;
+ for (i=s; i<=t; i++) {
+ if (getbit(b->B,i)==OP) {
+ c = '(';
+ d++;
+ } else {
+ c = ')';
+ d--;
+ }
+ putchar(c);
+ }
+}
+
+int *matchtbl,*parenttbl;
+void make_naivetbl(pb *B,int n)
+{
+ int i,v;
+ int *stack,s;
+
+ mymalloc(matchtbl,n,0);
+ mymalloc(parenttbl,n,0);
+ mymalloc(stack,n,0);
+
+ for (i=0; i<n; i++) matchtbl[i] = parenttbl[i] = -1;
+
+ s = 0;
+ v = 0;
+ stack[0] = -1;
+ for (i=0; i<n; i++) {
+ if (getbit(B,i)==OP) {
+ v++;
+ if (v > 0) {
+ parenttbl[i] = stack[v-1];
+ stack[v] = i;
+ }
+ } else {
+ if (v > 0) {
+ matchtbl[stack[v]] = i; // close
+ matchtbl[i] = stack[v]; // open
+ }
+ v--;
+ }
+ }
+
+ free(stack);
+}
+
+int popCount[1<<ETW];
+int fwdtbl[(2*ETW+1)*(1<<ETW)];
+int bwdtbl[(2*ETW+1)*(1<<ETW)];
+#if 0
+int mintbl_li[1<<ETW], mintbl_lv[1<<ETW];
+int mintbl_ri[1<<ETW], mintbl_rv[1<<ETW];
+int maxtbl_li[1<<ETW], maxtbl_lv[1<<ETW];
+int maxtbl_ri[1<<ETW], maxtbl_rv[1<<ETW];
+#endif
+int minmaxtbl_i[4][1<<ETW], minmaxtbl_v[4][1<<ETW];
+int degtbl[1<<ETW];
+int degtbl2[(2*ETW+1)*(1<<ETW)];
+int childtbl[(ETW)*(1<<ETW)];
+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;
+ int m,M;
+ pb buf[1];
+ int deg;
+ if (inited)
+ return;
+ inited = 1;
+ for (x = 0; x < (1<<ETW); x++) {
+ 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;
+ for (r=-ETW; r<=ETW; r++) depthtbl[((r+ETW)<<ETW)+x] = 0;
+
+ r = 0;
+ for (i=0; i<ETW; i++) {
+ if (getbit(buf,i)==OP) {
+ r++;
+ } else {
+ r--;
+ }
+ if (fwdtbl[((r+ETW)<<ETW)+x] == ETW) fwdtbl[((r+ETW)<<ETW)+x] = i;
+ }
+
+ r = 0;
+ for (i=ETW-1; i>=0; i--) {
+ if (getbit(buf,i)==OP) {
+ r--;
+ } else {
+ r++;
+ }
+ if (bwdtbl[((r+ETW)<<ETW)+x] == ETW) bwdtbl[((r+ETW)<<ETW)+x] = ETW-1-i;
+ }
+
+ r = 0;
+ for (i=0; i<ETW; i++) {
+ if (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;
+ //maxtbl_lv[x] = -ETW-1;
+ //mintbl_lv[x] = ETW+1;
+ minmaxtbl_v[OPT_MAX | OPT_LEFT][x] = -ETW-1;
+ minmaxtbl_v[OPT_MIN | OPT_LEFT][x] = ETW+1;
+ deg = 0;
+ for (i=0; i<ETW; i++) {
+ if (getbit(buf,i)==OP) {
+ r++;
+ if (r > M) {
+ M = r;
+ //maxtbl_li[x] = i; maxtbl_lv[x] = r;
+ minmaxtbl_i[OPT_MAX | OPT_LEFT][x] = i;
+ minmaxtbl_v[OPT_MAX | OPT_LEFT][x] = r;
+ }
+ } else {
+ r--;
+ if (r == m) {
+ deg++;
+ childtbl[((deg-1)<<ETW) + x] = i;
+ }
+ if (r < m) {
+ m = r;
+ //mintbl_li[x] = i; mintbl_lv[x] = r;
+ minmaxtbl_i[OPT_MIN | OPT_LEFT][x] = i;
+ minmaxtbl_v[OPT_MIN | OPT_LEFT][x] = r;
+ deg = 1;
+ childtbl[((deg-1)<<ETW) + x] = i;
+ }
+ }
+ if (r <= m) degtbl2[((r+ETW)<<ETW)+x]++;
+ }
+ degtbl[x] = deg;
+
+ r = 0; m = 0; M = 0;
+ //maxtbl_rv[x] = -ETW-1;
+ //mintbl_rv[x] = ETW+1;
+ 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) {
+ r++;
+ if (r >= M) {
+ M = r;
+ //maxtbl_ri[x] = i; maxtbl_rv[x] = r;
+ minmaxtbl_i[OPT_MAX | OPT_RIGHT][x] = i;
+ minmaxtbl_v[OPT_MAX | OPT_RIGHT][x] = r;
+ }
+ } else {
+ r--;
+ if (r <= m) {
+ m = r;
+ //mintbl_ri[x] = i; mintbl_rv[x] = r;
+ minmaxtbl_i[OPT_MIN | OPT_RIGHT][x] = i;
+ minmaxtbl_v[OPT_MIN | OPT_RIGHT][x] = r;
+ }
+ }
+ }
+
+ for (i = 0; i < ETW; i++) {
+ for (j = -ETW; j <= ETW; j++) {
+ childtbl2[j+ETW][i][x] = -1;
+ }
+ }
+
+ for (j=-ETW; j<=ETW; j++) {
+ int ith;
+ ith = 0;
+ r = 0;
+ for (i = 0; i < ETW; i++) {
+ if (getbit(buf,i)==OP) {
+ r++;
+ } else {
+ r--;
+ if (r < j) break;
+ if (r == j) {
+ ith++;
+ childtbl2[j+ETW][ith-1][x] = i;
+ }
+ }
+ }
+ }
+ }
+
+}
+
+
+int bp_construct(bp *b,int n, pb *B, int opt)
+{
+ int i,j,d;
+ int m,M,ds;
+ int ns,nm;
+ byte *sm, *sM;
+ byte *sd;
+ int *mm, *mM;
+ 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
+
+ b->B = B;
+ b->n = n;
+ b->opt = opt;
+ b->idx_size = 0;
+ b->sm = NULL;
+ b->sM = NULL;
+ b->sd = NULL;
+ b->mm = NULL;
+ b->mM = NULL;
+ b->md = NULL;
+ b->da_leaf = NULL;
+ 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->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);
+
+ 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);
+ b->sm = sm;
+ b->sM = sM;
+ if (opt & OPT_DEGREE) {
+ mymalloc(sd, ns, 0); 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);
+ m = M = ds;
+ r = 1;
+ } else {
+ d = depth(b,i);
+ if (d == m) r++;
+ if (d < m) {
+ m = d;
+ r = 1;
+ }
+ if (d > M) M = d;
+ }
+ if (i % SB == SB-1 || i==n-1) {
+ ds = 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 (M - ds + 1 < 0 || M - ds + 1 > 255) {
+ printf("error M=%d ds=%d\n",M,ds);
+ }
+ sm[i/SB] = m - ds + SB;
+ sM[i/SB] = M - ds + 1;
+ if (opt & OPT_DEGREE) sd[i/SB] = r;
+ }
+ }
+
+#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);
+ b->mm = mm;
+ b->mM = mM;
+ if (opt & OPT_DEGREE) {
+ mymalloc(md, nm + m_ofs, 0); 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);
+ if (i % MB == 0) {
+ m = M = d;
+ r = 1;
+ } else {
+ if (d == m) r++;
+ if (d < m) {
+ m = d;
+ r = 1;
+ }
+ if (d > M) M = d;
+ }
+ if (i % MB == MB-1 || i==n-1) {
+ mm[m_ofs+ i/MB] = m;
+ mM[m_ofs+ i/MB] = M;
+ 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 (j*2+1 < nm + m_ofs) m = min(m,mm[j*2+1]);
+ M = 0;
+ if (j*2 < nm + m_ofs) M = mM[j*2];
+ if (j*2+1 < nm + m_ofs) M = max(M,mM[j*2+1]);
+ mm[j] = m; mM[j] = M;
+ if (opt & OPT_DEGREE) {
+ d = 0;
+ if (j*2 < nm + m_ofs) d = md[j*2];
+ if (j*2+1 < nm + m_ofs) {
+ if (mm[j*2] == mm[j*2+1]) d += md[j*2+1];
+ if (mm[j*2] > mm[j*2+1]) d = md[j*2+1];
+ }
+ md[j] = d;
+ }
+ }
+ mm[0] = -1;
+ mM[0] = mM[1];
+ if (opt & OPT_DEGREE) {
+ md[0] = -1;
+ }
+
+
+#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;
+ } 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->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->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->idx_size += b->da_dfuds_leaf->idx_size;
+ } else {
+ b->da_dfuds_leaf = NULL;
+ }
+
+ return 0;
+}
+
+// destroyTree: frees the memory of tree.
+void destroyTree(bp *b) {
+ if (!b) return; // nothing to free
+
+ destroyDarray(b->da); // destroys da data structure
+ if (b->da) free(b->da);
+
+ 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);
+}
+
+
+// saveTree: saves parentheses data structure to file
+// By Diego Arroyuelo
+void saveTree(bp *b, FILE *fp) {
+
+ if (fwrite(&(b->n), sizeof(int), 1, fp) != 1) {
+ printf("Error: cannot save number of parentheses to file\n");
+ exit(1);
+ }
+
+ if (fwrite(b->B, sizeof(pb), (b->n+D-1)/D, fp) != ((b->n+D-1)/D)) {
+ printf("Error: cannot save parentheses sequence to file\n");
+ exit(1);
+ }
+
+ if (fwrite(&(b->opt), sizeof(int), 1, fp) != 1) {
+ printf("Error: cannot save opt in parentheses to file\n");
+ exit(1);
+ }
+}
+
+// loadTree: load parentheses data structure from file
+// By Diego Arroyuelo
+void loadTree(bp *b, FILE *fp) {
+
+ pb *B;
+ int n, opt;
+
+ if (fread(&n, sizeof(int), 1, fp) != 1) {
+ printf("Error: cannot read number of parentheses from file\n");
+ exit(1);
+ }
+
+ mymalloc(B,(n+D-1)/D,0);
+
+ 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);
+ }
+
+ if (fread(&opt, sizeof(int), 1, fp) != 1) {
+ printf("Error: cannot read opt in parentheses from file\n");
+ exit(1);
+ }
+
+ bp_construct(b, n, B, opt);
+
+}
+
+
+
+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) {
+ v++;
+ } else {
+ v--;
+ }
+ if (v == rel) return i;
+ }
+ return -1;
+}
+
+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) {
+ v--;
+ } else {
+ v++;
+ }
+ if (v == rel) return i-1;
+ }
+ return -2;
+}
+
+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) {
+ rel++;
+ } else {
+ rel--;
+ }
+ if (rel == 0) return i-1;
+ }
+ if (i < 0) return -2;
+ return -3;
+}
+
+int 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;
+ i = t-1;
+ while (i >= s) {
+ if (getbit(b->B,i+1)==CP) {
+ d++;
+ if (opt & OPT_MAX) {
+ if (d > dm) {
+ dm = d; im = i;
+ }
+ }
+ } else {
+ d--;
+ if (!(opt & OPT_MAX)) {
+ if (d < dm) {
+ dm = d; im = i;
+ }
+ }
+ }
+ i--;
+ }
+ } else {
+ d = dm = depth(b,s); im = s;
+ i = s+1;
+ while (i <= t) {
+ if (getbit(b->B,i)==OP) {
+ d++;
+ if (opt & OPT_MAX) {
+ if (d > dm) {
+ dm = d; im = i;
+ }
+ }
+ } else {
+ d--;
+ if (!(opt & OPT_MAX)) {
+ if (d < dm) {
+ dm = d; im = i;
+ }
+ }
+ }
+ i++;
+ }
+ }
+ return im;
+}
+
+int root_node(bp *b)
+{
+ return 0;
+}
+
+
+int rank_open(bp *b, int s)
+{
+ return darray_rank(b->da,s);
+}
+
+int rank_close(bp *b, int s)
+{
+ return s+1 - darray_rank(b->da,s);
+}
+
+int select_open(bp *b, int s)
+{
+ if (b->opt & OPT_FAST_PREORDER_SELECT) {
+ return darray_select(b->da,s,1);
+ } else {
+ return darray_select_bsearch(b->da,s,getpat_preorder);
+ }
+}
+
+int select_close(bp *b, int s)
+{
+ if (b->opt & OPT_FAST_POSTORDER_SELECT) {
+ return 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)
+// 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 t;
+ if (inspect(b,s) == CP) return -1;
+ t = find_close(b,s);
+ // return t+1 - darray_rank(b->da,t);
+ return rank_close(b,t);
+}
+
+int postorder_select_bsearch(bp *b,int s)
+{
+ int l,r,m;
+
+ if (s == 0) return -1;
+
+ if (s > b->da->n - b->da->m) {
+ return -1;
+ }
+ l = 0; r = b->da->n - 1;
+
+ 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) {
+ r = m;
+ } else {
+ l = m+1;
+ }
+ }
+ return l;
+}
+
+///////////////////////////////////////////
+// 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)
+{
+#if 0
+ if (b->opt & OPT_FAST_POSTORDER_SELECT) {
+ return darray_pat_select(b->da_postorder,s,getpat_postorder);
+ } else {
+ return postorder_select_bsearch(b->da,s);
+ }
+#else
+ return select_close(b,s);
+#endif
+}
+
+///////////////////////////////////////////
+// leaf_rank(bp *b,int s)
+// returns the number of leaves to the left of s
+///////////////////////////////////////////
+int leaf_rank(bp *b,int s)
+{
+ if ((b->opt & OPT_LEAF) == 0) {
+ printf("leaf_rank: error!!! not supported\n");
+ return -1;
+ }
+ if (s >= b->n-1) {
+ s = b->n-2;
+ }
+ return 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)
+{
+ if ((b->opt & OPT_LEAF) == 0) {
+ printf("leaf_select: error!!! not supported\n");
+ return -1;
+ }
+ 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);
+ } else {
+ return darray_select_bsearch(b->da_leaf,s,getpat_leaf);
+ }
+}
+
+
+///////////////////////////////////////////
+// inorder_rank(bp *b,int s)
+// returns the number of ")(" (s >= 0)
+///////////////////////////////////////////
+int inorder_rank(bp *b,int s)
+{
+ if ((b->opt & OPT_INORDER) == 0) {
+ printf("inorder_rank: error!!! not supported\n");
+ return -1;
+ }
+ if (s >= b->n-1) {
+ s = b->n-2;
+ }
+ return darray_pat_rank(b->da_inorder,s,getpat_inorder);
+}
+
+///////////////////////////////////////////
+// inorder_select(bp *b,int s)
+// returns the s-th position of ")(" (s >= 1)
+///////////////////////////////////////////
+int 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);
+ } else {
+ return darray_select_bsearch(b->da_inorder,s,getpat_inorder);
+ }
+}
+
+///////////////////////////////////////////
+// leftmost_leaf(bp *b, int s)
+///////////////////////////////////////////
+int 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);
+}
+
+///////////////////////////////////////////
+// rightmost_leaf(bp *b, int s)
+///////////////////////////////////////////
+int 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);
+}
+
+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)
+// returns the previous sibling of parent(s)
+// -1 if s is the first child
+//////////////////////////////////////////
+int 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);
+ return t;
+}
+
+///////////////////////////////////////////
+// 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 t,m;
+ t = find_close(b,s);
+ m = rmq(b,s,t, OPT_MAX);
+ return m;
+}
+
+///////////////////////////////////////////
+// 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 t;
+ t = deepest_node(b,s);
+ return depth(b,t) - depth(b,s);
+}
+
+int naive_degree(bp *b, int s)
+{
+ int t,d;
+ d = 0;
+ t = first_child(b,s);
+ while (t >= 0) {
+ d++;
+ t = next_sibling(b,t);
+ }
+ return d;
+}
+
+///////////////////////////////////////////
+// degree(bp *b, int s)
+// returns the number of children of s
+// 0 if s is a leaf
+///////////////////////////////////////////
+int degree(bp *b, int s)
+{
+ if (b->opt & OPT_DEGREE) {
+ return fast_degree(b,s,b->n,0);
+ } else {
+ return naive_degree(b,s);
+ }
+}
+
+int naive_child(bp *b, int s, int d)
+{
+ int t,i;
+ t = first_child(b,s);
+ for (i = 1; i < d; i++) {
+ if (t == -1) break;
+ t = next_sibling(b,t);
+ }
+ return t;
+}
+
+///////////////////////////////////////////
+// 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 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;
+ return r;
+ } else {
+ return naive_child(b,s,d);
+ }
+
+}
+
+
+int naive_child_rank(bp *b, int t)
+{
+ int v,d;
+ d = 0;
+ while (t != -1) {
+ d++;
+ t = prev_sibling(b,t);
+ }
+ return d;
+}
+
+///////////////////////////////////////////
+// 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 r;
+ if (t == root_node(b)) return 1;
+ if (b->opt & OPT_DEGREE) {
+ r = parent(b,t);
+ return fast_degree(b,r,t,0)+1;
+ } else {
+ return 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 v,d;
+ v = lca(b,s,t);
+ d = depth(b,v);
+ return (depth(b,s) - d) + (depth(b,t) - d);
+}
+
+///////////////////////////////////////////
+// level_next(bp *b, int d)
+///////////////////////////////////////////
+int level_next(bp *b,int s)
+{
+ int t;
+ t = fwd_excess(b,s,0);
+ return t;
+}
+
+///////////////////////////////////////////
+// level_prev(bp *b, int d)
+///////////////////////////////////////////
+int level_prev(bp *b,int s)
+{
+ int t;
+ t = bwd_excess(b,s,0);
+ return t;
+}
+
+///////////////////////////////////////////
+// level_leftmost(bp *b, int d)
+///////////////////////////////////////////
+int level_leftmost(bp *b, int d)
+{
+ int t;
+ if (d < 1) return -1;
+ if (d == 1) return 0;
+ t = fwd_excess(b,0,d);
+ return t;
+}
+
+///////////////////////////////////////////
+// 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);
+}
+
+///////////////////////////////////////////
+// leaf_size(bp *b, int s)
+///////////////////////////////////////////
+int 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);
+}
--- /dev/null
+#ifndef BP_H_
+#define BP_H_
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+
+#include <stdio.h>
+#include <stdlib.h>
+#include "darray.h"
+
+#define OP 1
+#define CP 0
+
+#define OPT_MIN 0
+#define OPT_MAX 1
+#define OPT_LEFT 0
+#define OPT_RIGHT 2
+
+#define OPT_LEAF (1<<0)
+#define OPT_INORDER (1<<1)
+#define OPT_DEGREE (1<<2)
+#define OPT_FAST_PREORDER_SELECT (1<<3)
+#define OPT_FAST_LEAF_SELECT (1<<4)
+#define OPT_FAST_INORDER_SELECT (1<<5)
+#define OPT_FAST_POSTORDER_SELECT (1<<6)
+#define OPT_DFUDS_LEAF (1<<7)
+#define OPT_FAST_DFUDS_LEAF_SELECT (1<<8)
+
+//#define logSB 9
+#define logSB 7
+//#define logSB 2
+#define SB (1<<logSB)
+//#define logMB 16
+//#define logMB 12
+#define logMB 15
+//#define logMB 3
+#define MB (1<<logMB)
+
+#define ETW 8 // width of excess lookup table
+#define W1 2 // branching factor
+
+
+
+typedef struct {
+ int n;
+ pb *B;
+ darray *da;
+ byte *sm, *sM;
+ byte *sd;
+ int *mm, *mM;
+ int *md;
+ int m_ofs;
+
+ darray *da_leaf;
+ darray *da_inorder;
+ darray *da_postorder;
+ darray *da_dfuds_leaf;
+ int idx_size;
+ 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);
+
+// 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);
+
+// 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);
+
+// using leaf select index
+int leaf_select(bp *b,int s);
+
+// using inorder index
+int inorder_rank(bp *b,int s);
+
+// using inorder select index
+int 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);
+
+
+// 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);
+
+
+inline int blog(int x)
+{
+ int l;
+ l = 0;
+ while (x>0) {
+ x>>=1;
+ l++;
+ }
+ 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];
+extern int mintbl_ri[1<<ETW], mintbl_rv[1<<ETW];
+extern int maxtbl_li[1<<ETW], maxtbl_lv[1<<ETW];
+extern int maxtbl_ri[1<<ETW], maxtbl_rv[1<<ETW];
+
+extern int minmaxtbl_i[4][1<<ETW], minmaxtbl_v[4][1<<ETW];
+extern int degtbl[1<<ETW];
+extern int degtbl2[(2*ETW+1)*(1<<ETW)];
+extern int childtbl[(ETW)*(1<<ETW)];
+extern int depthtbl[(2*ETW+1)*(1<<ETW)];
+extern int childtbl2[2*ETW+1][ETW][(1<<ETW)];
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif
--- /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