+++ /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