5 //typedef unsigned char byte;
\r
6 //typedef unsigned short word;
\r
7 //typedef unsigned int dword;
\r
12 #define PBS (sizeof(pb)*8)
\r
18 #define logLL 16 // size of word
\r
19 #define LL (1<<logLL)
\r
22 #define LLL (1<<logLLL)
\r
24 //#define logL (logLL-3)
\r
25 #define logL (logLL-1-5)
\r
29 #define min(x,y) ((x)<(y)?(x):(y))
\r
32 #define mymalloc(p,n,f) {p = (__typeof__(p)) malloc((n)*sizeof(*p)); if ((p)==NULL) {printf("not enough memory\n"); exit(1);}}
\r
34 int setbit(pb *B, int i,int x)
\r
40 if (x==0) B[j] &= (~(1<<(D-1-l)));
\r
41 else if (x==1) B[j] |= (1<<(D-1-l));
\r
43 printf("error setbit x=%d\n",x);
\r
49 int setbit_zero(pb *B, int i)
\r
54 B[j] &= (~(1<<(D-1-l)));
\r
58 int setbit_one(pb *B, int i)
\r
63 B[j] |= (1<<(D-1-l));
\r
70 int setbits(pb *B, int i, int d, int x)
\r
74 for (j=0; j<d; j++) {
\r
75 setbit(B,i+j,(x>>(d-j-1))&1);
\r
80 int getbit(pb *B, int i)
\r
88 return (B[j] >> (D-1-l)) & 1;
\r
91 dword getbits(pb *B, int i, int d)
\r
96 for (j=0; j<d; j++) {
\r
103 int getpattern(pb *B, int i, int k, pb pat)
\r
108 for (j=0; j<k; j++) {
\r
109 x &= getbit(B,i+j) ^ (~(pat>>(k-1-j)));
\r
111 //printf("getpattern(%d) = %d\n",i,x);
\r
116 static const unsigned int popCount[] = {
\r
117 0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
\r
118 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
\r
119 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
\r
120 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
121 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
\r
122 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
123 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
124 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
\r
125 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
\r
126 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
127 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
128 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
\r
129 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
130 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
\r
131 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
\r
132 4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8
\r
135 static int selecttbl[8*256];
\r
137 unsigned int popcount_old(pb x)
\r
142 r = r - ((r>>1) & 0x77777777) - ((r>>2) & 0x33333333) - ((r>>3) & 0x11111111);
\r
143 r = ((r + (r>>4)) & 0x0f0f0f0f) % 0xff;
\r
146 r = ((r & 0xaaaaaaaa)>>1) + (r & 0x55555555);
\r
147 r = ((r & 0xcccccccc)>>2) + (r & 0x33333333);
\r
148 //r = ((r & 0xf0f0f0f0)>>4) + (r & 0x0f0f0f0f);
\r
149 r = ((r>>4) + r) & 0x0f0f0f0f;
\r
150 //r = ((r & 0xff00ff00)>>8) + (r & 0x00ff00ff);
\r
152 //r = ((r & 0xffff0000)>>16) + (r & 0x0000ffff);
\r
153 r = ((r>>16) + r) & 63;
\r
155 r = popCount[x & 0xff];
\r
157 r += popCount[x & 0xff];
\r
159 r += popCount[x & 0xff];
\r
161 r += popCount[x & 0xff];
\r
166 inline unsigned int
\r
169 uint m1 = 0x55555555;
\r
170 uint m2 = 0xc30c30c3;
\r
171 x -= (x >> 1) & m1;
\r
172 x = (x & m2) + ((x >> 2) & m2) + ((x >> 4) & m2);
\r
174 return (x + (x >> 12) + (x >> 24)) & 0x3f;
\r
178 unsigned int popcount8(pb x)
\r
183 r = ((r & 0xaa)>>1) + (r & 0x55);
\r
184 r = ((r & 0xcc)>>2) + (r & 0x33);
\r
185 r = ((r>>4) + r) & 0x0f;
\r
187 r = popCount[x & 0xff];
\r
192 void make_selecttbl(void)
\r
197 for (x = 0; x < 256; x++) {
\r
198 setbits(buf,0,8,x);
\r
199 for (r=0; r<8; r++) selecttbl[(r<<8)+x] = -1;
\r
201 for (i=0; i<8; i++) {
\r
202 if (getbit(buf,i)) {
\r
203 selecttbl[(r<<8)+x] = i;
\r
211 int darray_construct(darray *da, int n, pb *buf, int opt)
\r
219 int p1,p2,p3,p4,s1,s2,s3,s4;
\r
227 printf("ERROR: L=%d LLL=%d\n",L,LLL);
\r
232 for (i=0; i<n; i++) m += getbit(buf,i);
\r
235 //printf("n=%d m=%d\n",n,m);
\r
239 if (opt & (~OPT_NO_RANK)) { // construct select table
\r
243 for (i=0; i<n; i++) {
\r
244 if (getbit(buf,i)) {
\r
250 nl = (m-1) / L + 1;
\r
251 mymalloc(da->lp,nl+1,1); da->idx_size += (nl+1)*sizeof(*da->lp);
\r
252 mymalloc(da->p,nl+1,1); da->idx_size += (nl+1)*sizeof(*da->p);
\r
254 printf("lp table: %d bytes (%1.2f bpc)\n",(nl+1)*sizeof(*da->lp), (double)(nl+1)*sizeof(*da->lp) * 8/n);
\r
255 printf("p table: %d bytes (%1.2f bpc)\n",(nl+1)*sizeof(*da->p), (double)(nl+1)*sizeof(*da->p) * 8/n);
\r
258 for (r = 0; r < 2; r++) {
\r
259 s1 = s2 = s3 = s4 = 0;
\r
260 p1 = p2 = p3 = p4 = -1;
\r
263 for (il = 0; il < nl; il++) {
\r
265 while (s1 <= il*L) {
\r
266 if (getbit(buf,p1+1)) s1++;
\r
271 i = min((il+1)*L-1,m-1);
\r
274 if (getbit(buf,p2+1)) s2++;
\r
278 //printf("%d ",p-pp);
\r
279 if (p - pp >= LL) {
\r
281 for (is = 0; is < L; is++) {
\r
282 if (il*L+is >= m) break;
\r
283 //da->sl[ml*L+is] = s[il*L+is];
\r
284 while (s3 <= il*L+is) {
\r
285 if (getbit(buf,p3+1)) s3++;
\r
288 da->sl[ml*L+is] = p3;
\r
291 da->p[il] = -(ml+1);
\r
295 for (is = 0; is < L/LLL; is++) {
\r
296 if (il*L+is*LLL >= m) break;
\r
297 while (s4 <= il*L+is*LLL) {
\r
298 if (getbit(buf,p4+1)) s4++;
\r
301 //da->ss[ms*(L/LLL)+is] = s[il*L+is*LLL] - pp;
\r
302 da->ss[ms*(L/LLL)+is] = p4 - pp;
\r
310 mymalloc(da->sl,ml*L+1,1); da->idx_size += (ml*L+1)*sizeof(*da->sl);
\r
311 mymalloc(da->ss,ms*(L/LLL)+1,1); da->idx_size += (ms*(L/LLL)+1)*sizeof(*da->ss);
\r
313 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
314 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
319 } else { // no select table
\r
321 da->p = da->sl = NULL;
\r
325 // construct rank table
\r
327 if ((opt & OPT_NO_RANK) == 0) {
\r
328 mymalloc(da->rl,n/R1+2,1); da->idx_size += (n/R1+2)*sizeof(*da->rl);
\r
329 mymalloc(da->rm,n/RR+2,1); da->idx_size += (n/RR+2)*sizeof(*da->rm);
\r
330 mymalloc(da->rs,n/RRR+2,1); da->idx_size += (n/RRR+2)*sizeof(*da->rs);
\r
332 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
333 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
334 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
337 for (k=0; k<=n; k+=R1) {
\r
340 for (i=0; i<R1; i+=RR) {
\r
341 if (k+i <= n) da->rm[(k+i)/RR] = m2;
\r
343 for (j=0; j<RR; j++) {
\r
344 if (k+i+j < n && getbit(buf,k+i+j)==1) m++;
\r
345 if (j % RRR == RRR-1) {
\r
346 if (k+i+j <= n) da->rs[(k+i+j)/RRR] = m;
\r
363 // destroyDarray: frees the memory of darray
\r
364 // Added by Diego Arroyuelo
\r
365 void destroyDarray(darray *da) {
\r
366 if (!da) return; // nothing to free
\r
368 if (da->buf) free(da->buf);
\r
369 if (da->lp) free(da->lp);
\r
370 if (da->sl) free(da->sl);
\r
371 if (da->ss) free(da->ss);
\r
372 if (da->p) free(da->p);
\r
373 if (da->rl) free(da->rl);
\r
374 if (da->rm) free(da->rm);
\r
375 if (da->rs) free(da->rs);
\r
378 int darray_rank0(darray *da, int i)
\r
384 r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];
\r
385 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
387 if (j < D) r += popcount(*p >> (D-1-j));
\r
388 else r += popcount(*p) + popcount(p[1] >> (D-1-(j-D)));
\r
393 r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];
\r
394 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
396 r += popcount(*p++);
\r
399 r += popcount(*p >> (D-1-j));
\r
401 j = RRR-1 - (i & (RRR-1));
\r
403 r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];
\r
404 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
406 r -= popcount(*--p);
\r
409 if (j > 0) r -= popcount(*--p << (D-j));
\r
417 int darray_rank(darray *da, int i)
\r
422 r = da->rl[i>>logR] + da->rm[i>>logRR];
\r
423 j = (i>>logRRR) & (RR/RRR-1);
\r
425 r += da->rs[((i>>logRR)<<(logRR-logRRR))+j-1];
\r
429 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
432 r += popcount(*p++);
\r
435 r += popcount(*p >> (D-1-j));
\r
440 int darray_select_bsearch(darray *da, int i, pb (*getpat)(pb *))
\r
449 //s = darray_select(da,i,1);
\r
451 //printf("select(%d)=%d\n",i,s);
\r
455 if (i == 0) return -1;
\r
469 l = 0; r = (n-1)>>logR;
\r
470 // find the smallest index x s.t. rl[x] >= t
\r
473 //printf("m=%d rl[m+1]=%d t=%d\n",m,da->rl[m+1],t);
\r
474 if (da->rl[m+1] > t) { // m+1 is out of range
\r
475 r = m; // new r = m >= l
\r
477 l = m+1; // new l = m+1 <= r
\r
485 l = x >> logRR; r = (min(x+R1-1,n))>>logRR;
\r
488 //printf("m=%d rm[m+1]=%d t=%d\n",m,da->rm[m+1],t);
\r
489 if (da->rm[m+1] > t) { // m+1 is out of range
\r
492 l = m+1; // new l = m+1 <= r
\r
501 l = x >> logRRR; r = (min(x+RR-1,n))>>logRRR;
\r
504 //printf("m=%d rs[m+1]=%d t=%d\n",m,da->rs[m+1],t);
\r
505 if (da->rs[m+1] > t) { // m+1 is out of range
\r
508 l = m+1; // new l = m+1 <= r
\r
515 while (t > da->rs[l]) {
\r
524 p = &da->buf[x >> logD];
\r
526 m = popcount(getpat(p));
\r
535 rr = popCount[w >> (D-8)];
\r
541 x += selecttbl[((t-0)<<8)+(w>>(D-8))];
\r
545 printf("error x=%d s=%d\n",x,s);
\r
551 int darray_pat_rank(darray *da, int i, pb (*getpat)(pb *))
\r
556 r = da->rl[i>>logR] + da->rm[i>>logRR];
\r
557 j = (i>>logRRR) & (RR/RRR-1);
\r
559 r += da->rs[((i>>logRR)<<(logRR-logRRR))+j-1];
\r
563 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
566 r += popcount(getpat(p));
\r
570 r += popcount(getpat(p) >> (D-1-j));
\r
576 int darray_select(darray *da, int i,int f)
\r
584 if (i == 0) return -1;
\r
588 //printf("ERROR: m=%d i=%d\n",da->m,i);
\r
594 il = da->p[i>>logL];
\r
597 p = da->sl[(il<<logL)+(i & (L-1))];
\r
599 p = da->lp[i>>logL];
\r
600 p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
\r
601 r = i - (i & (LLL-1));
\r
603 q = &(da->buf[p>>logD]);
\r
607 r -= popcount(*q >> (D-1-rr));
\r
612 if (r + rr >= i) break;
\r
620 //rr = popcount(x >> (D-8));
\r
621 rr = popCount[x >> (D-8)];
\r
622 //rr = popcount8(x >> (D-8));
\r
623 if (r + rr >= i) break;
\r
628 p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];
\r
631 r -= popcount((~(*q)) >> (D-1-rr));
\r
635 rr = popcount(~(*q));
\r
636 if (r + rr >= i) break;
\r
645 //rr = popcount(x >> (D-8));
\r
646 rr = popCount[x >> (D-8)];
\r
647 //rr = popcount8(x >> (D-8));
\r
648 if (r + rr >= i) break;
\r
653 p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];
\r
659 int darray_pat_select(darray *da, int i, pb (*getpat)(pb *))
\r
667 if (i == 0) return -1;
\r
671 //printf("ERROR: m=%d i=%d\n",da->m,i);
\r
677 il = da->p[i>>logL];
\r
680 p = da->sl[(il<<logL)+(i & (L-1))];
\r
682 p = da->lp[i>>logL];
\r
683 p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
\r
684 r = i - (i & (LLL-1));
\r
686 q = &(da->buf[p>>logD]);
\r
689 r -= popcount(getpat(q) >> (D-1-rr));
\r
693 rr = popcount(getpat(q));
\r
694 if (r + rr >= i) break;
\r
702 //rr = popcount(x >> (D-8));
\r
703 rr = popCount[x >> (D-8)];
\r
704 //rr = popcount8(x >> (D-8));
\r
705 if (r + rr >= i) break;
\r
710 p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];
\r
715 int darray_pat_construct(darray *da, int n, pb *buf, int k, pb pat, int opt)
\r
719 mymalloc(b,(n+D-1)/D,0);
\r
721 for (i=0; i<n-k+1; i++) {
\r
722 setbit(b,i,getpattern(buf,i,k,pat));
\r
724 for (i=n-k+1; i<n; i++) {
\r
728 darray_construct(da,n,b,opt);
\r