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