--- /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
projects / SXSI / libbp.git / blobdiff