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 setbits(pb *B, int i, int d, int x)
\r
53 for (j=0; j<d; j++) {
\r
54 setbit(B,i+j,(x>>(d-j-1))&1);
\r
59 int getbit(pb *B, int i)
\r
67 return (B[j] >> (D-1-l)) & 1;
\r
70 dword getbits(pb *B, int i, int d)
\r
75 for (j=0; j<d; j++) {
\r
82 int getpattern(pb *B, int i, int k, pb pat)
\r
87 for (j=0; j<k; j++) {
\r
88 x &= getbit(B,i+j) ^ (~(pat>>(k-1-j)));
\r
90 //printf("getpattern(%d) = %d\n",i,x);
\r
95 static const unsigned int popCount[] = {
\r
96 0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
\r
97 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
\r
98 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
\r
99 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
100 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
\r
101 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
102 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
103 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
\r
104 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
\r
105 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
106 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
107 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
\r
108 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
\r
109 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
\r
110 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
\r
111 4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8
\r
114 static int selecttbl[8*256];
\r
116 unsigned int popcount(pb x)
\r
121 r = r - ((r>>1) & 0x77777777) - ((r>>2) & 0x33333333) - ((r>>3) & 0x11111111);
\r
122 r = ((r + (r>>4)) & 0x0f0f0f0f) % 0xff;
\r
125 r = ((r & 0xaaaaaaaa)>>1) + (r & 0x55555555);
\r
126 r = ((r & 0xcccccccc)>>2) + (r & 0x33333333);
\r
127 //r = ((r & 0xf0f0f0f0)>>4) + (r & 0x0f0f0f0f);
\r
128 r = ((r>>4) + r) & 0x0f0f0f0f;
\r
129 //r = ((r & 0xff00ff00)>>8) + (r & 0x00ff00ff);
\r
131 //r = ((r & 0xffff0000)>>16) + (r & 0x0000ffff);
\r
132 r = ((r>>16) + r) & 63;
\r
134 r = popCount[x & 0xff];
\r
136 r += popCount[x & 0xff];
\r
138 r += popCount[x & 0xff];
\r
140 r += popCount[x & 0xff];
\r
145 unsigned int popcount8(pb x)
\r
150 r = ((r & 0xaa)>>1) + (r & 0x55);
\r
151 r = ((r & 0xcc)>>2) + (r & 0x33);
\r
152 r = ((r>>4) + r) & 0x0f;
\r
154 r = popCount[x & 0xff];
\r
159 void make_selecttbl(void)
\r
164 for (x = 0; x < 256; x++) {
\r
165 setbits(buf,0,8,x);
\r
166 for (r=0; r<8; r++) selecttbl[(r<<8)+x] = -1;
\r
168 for (i=0; i<8; i++) {
\r
169 if (getbit(buf,i)) {
\r
170 selecttbl[(r<<8)+x] = i;
\r
178 int darray_construct(darray *da, int n, pb *buf, int opt)
\r
186 int p1,p2,p3,p4,s1,s2,s3,s4;
\r
194 printf("ERROR: L=%d LLL=%d\n",L,LLL);
\r
199 for (i=0; i<n; i++) m += getbit(buf,i);
\r
202 //printf("n=%d m=%d\n",n,m);
\r
206 if (opt & (~OPT_NO_RANK)) { // construct select table
\r
210 for (i=0; i<n; i++) {
\r
211 if (getbit(buf,i)) {
\r
217 nl = (m-1) / L + 1;
\r
218 mymalloc(da->lp,nl+1,1); da->idx_size += (nl+1)*sizeof(*da->lp);
\r
219 mymalloc(da->p,nl+1,1); da->idx_size += (nl+1)*sizeof(*da->p);
\r
221 printf("lp table: %d bytes (%1.2f bpc)\n",(nl+1)*sizeof(*da->lp), (double)(nl+1)*sizeof(*da->lp) * 8/n);
\r
222 printf("p table: %d bytes (%1.2f bpc)\n",(nl+1)*sizeof(*da->p), (double)(nl+1)*sizeof(*da->p) * 8/n);
\r
225 for (r = 0; r < 2; r++) {
\r
226 s1 = s2 = s3 = s4 = 0;
\r
227 p1 = p2 = p3 = p4 = -1;
\r
230 for (il = 0; il < nl; il++) {
\r
232 while (s1 <= il*L) {
\r
233 if (getbit(buf,p1+1)) s1++;
\r
238 i = min((il+1)*L-1,m-1);
\r
241 if (getbit(buf,p2+1)) s2++;
\r
245 //printf("%d ",p-pp);
\r
246 if (p - pp >= LL) {
\r
248 for (is = 0; is < L; is++) {
\r
249 if (il*L+is >= m) break;
\r
250 //da->sl[ml*L+is] = s[il*L+is];
\r
251 while (s3 <= il*L+is) {
\r
252 if (getbit(buf,p3+1)) s3++;
\r
255 da->sl[ml*L+is] = p3;
\r
258 da->p[il] = -(ml+1);
\r
262 for (is = 0; is < L/LLL; is++) {
\r
263 if (il*L+is*LLL >= m) break;
\r
264 while (s4 <= il*L+is*LLL) {
\r
265 if (getbit(buf,p4+1)) s4++;
\r
268 //da->ss[ms*(L/LLL)+is] = s[il*L+is*LLL] - pp;
\r
269 da->ss[ms*(L/LLL)+is] = p4 - pp;
\r
277 mymalloc(da->sl,ml*L+1,1); da->idx_size += (ml*L+1)*sizeof(*da->sl);
\r
278 mymalloc(da->ss,ms*(L/LLL)+1,1); da->idx_size += (ms*(L/LLL)+1)*sizeof(*da->ss);
\r
280 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
281 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
286 } else { // no select table
\r
288 da->p = da->sl = NULL;
\r
292 // construct rank table
\r
294 if ((opt & OPT_NO_RANK) == 0) {
\r
295 mymalloc(da->rl,n/R1+2,1); da->idx_size += (n/R1+2)*sizeof(*da->rl);
\r
296 mymalloc(da->rm,n/RR+2,1); da->idx_size += (n/RR+2)*sizeof(*da->rm);
\r
297 mymalloc(da->rs,n/RRR+2,1); da->idx_size += (n/RRR+2)*sizeof(*da->rs);
\r
299 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
300 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
301 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
304 for (k=0; k<=n; k+=R1) {
\r
307 for (i=0; i<R1; i+=RR) {
\r
308 if (k+i <= n) da->rm[(k+i)/RR] = m2;
\r
310 for (j=0; j<RR; j++) {
\r
311 if (k+i+j < n && getbit(buf,k+i+j)==1) m++;
\r
312 if (j % RRR == RRR-1) {
\r
313 if (k+i+j <= n) da->rs[(k+i+j)/RRR] = m;
\r
330 // destroyDarray: frees the memory of darray
\r
331 // Added by Diego Arroyuelo
\r
332 void destroyDarray(darray *da) {
\r
333 if (!da) return; // nothing to free
\r
335 if (da->buf) free(da->buf);
\r
336 if (da->lp) free(da->lp);
\r
337 if (da->sl) free(da->sl);
\r
338 if (da->ss) free(da->ss);
\r
339 if (da->p) free(da->p);
\r
340 if (da->rl) free(da->rl);
\r
341 if (da->rm) free(da->rm);
\r
342 if (da->rs) free(da->rs);
\r
345 int darray_rank0(darray *da, int i)
\r
351 r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];
\r
352 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
354 if (j < D) r += popcount(*p >> (D-1-j));
\r
355 else r += popcount(*p) + popcount(p[1] >> (D-1-(j-D)));
\r
360 r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];
\r
361 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
363 r += popcount(*p++);
\r
366 r += popcount(*p >> (D-1-j));
\r
368 j = RRR-1 - (i & (RRR-1));
\r
370 r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];
\r
371 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
373 r -= popcount(*--p);
\r
376 if (j > 0) r -= popcount(*--p << (D-j));
\r
384 int darray_rank(darray *da, int i)
\r
389 r = da->rl[i>>logR] + da->rm[i>>logRR];
\r
390 j = (i>>logRRR) & (RR/RRR-1);
\r
392 r += da->rs[((i>>logRR)<<(logRR-logRRR))+j-1];
\r
396 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
399 r += popcount(*p++);
\r
402 r += popcount(*p >> (D-1-j));
\r
407 int darray_select_bsearch(darray *da, int i, pb (*getpat)(pb *))
\r
416 //s = darray_select(da,i,1);
\r
418 //printf("select(%d)=%d\n",i,s);
\r
422 if (i == 0) return -1;
\r
436 l = 0; r = (n-1)>>logR;
\r
437 // find the smallest index x s.t. rl[x] >= t
\r
440 //printf("m=%d rl[m+1]=%d t=%d\n",m,da->rl[m+1],t);
\r
441 if (da->rl[m+1] > t) { // m+1 is out of range
\r
442 r = m; // new r = m >= l
\r
444 l = m+1; // new l = m+1 <= r
\r
452 l = x >> logRR; r = (min(x+R1-1,n))>>logRR;
\r
455 //printf("m=%d rm[m+1]=%d t=%d\n",m,da->rm[m+1],t);
\r
456 if (da->rm[m+1] > t) { // m+1 is out of range
\r
459 l = m+1; // new l = m+1 <= r
\r
468 l = x >> logRRR; r = (min(x+RR-1,n))>>logRRR;
\r
471 //printf("m=%d rs[m+1]=%d t=%d\n",m,da->rs[m+1],t);
\r
472 if (da->rs[m+1] > t) { // m+1 is out of range
\r
475 l = m+1; // new l = m+1 <= r
\r
482 while (t > da->rs[l]) {
\r
491 p = &da->buf[x >> logD];
\r
493 m = popcount(getpat(p));
\r
502 rr = popCount[w >> (D-8)];
\r
508 x += selecttbl[((t-0)<<8)+(w>>(D-8))];
\r
512 printf("error x=%d s=%d\n",x,s);
\r
518 int darray_pat_rank(darray *da, int i, pb (*getpat)(pb *))
\r
523 r = da->rl[i>>logR] + da->rm[i>>logRR];
\r
524 j = (i>>logRRR) & (RR/RRR-1);
\r
526 r += da->rs[((i>>logRR)<<(logRR-logRRR))+j-1];
\r
530 p = da->buf + ((i>>logRRR)<<(logRRR-logD));
\r
533 r += popcount(getpat(p));
\r
537 r += popcount(getpat(p) >> (D-1-j));
\r
543 int darray_select(darray *da, int i,int f)
\r
551 if (i == 0) return -1;
\r
555 //printf("ERROR: m=%d i=%d\n",da->m,i);
\r
561 il = da->p[i>>logL];
\r
564 p = da->sl[(il<<logL)+(i & (L-1))];
\r
566 p = da->lp[i>>logL];
\r
567 p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
\r
568 r = i - (i & (LLL-1));
\r
570 q = &(da->buf[p>>logD]);
\r
574 r -= popcount(*q >> (D-1-rr));
\r
579 if (r + rr >= i) break;
\r
587 //rr = popcount(x >> (D-8));
\r
588 rr = popCount[x >> (D-8)];
\r
589 //rr = popcount8(x >> (D-8));
\r
590 if (r + rr >= i) break;
\r
595 p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];
\r
598 r -= popcount((~(*q)) >> (D-1-rr));
\r
602 rr = popcount(~(*q));
\r
603 if (r + rr >= i) break;
\r
612 //rr = popcount(x >> (D-8));
\r
613 rr = popCount[x >> (D-8)];
\r
614 //rr = popcount8(x >> (D-8));
\r
615 if (r + rr >= i) break;
\r
620 p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];
\r
626 int darray_pat_select(darray *da, int i, pb (*getpat)(pb *))
\r
634 if (i == 0) return -1;
\r
638 //printf("ERROR: m=%d i=%d\n",da->m,i);
\r
644 il = da->p[i>>logL];
\r
647 p = da->sl[(il<<logL)+(i & (L-1))];
\r
649 p = da->lp[i>>logL];
\r
650 p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
\r
651 r = i - (i & (LLL-1));
\r
653 q = &(da->buf[p>>logD]);
\r
656 r -= popcount(getpat(q) >> (D-1-rr));
\r
660 rr = popcount(getpat(q));
\r
661 if (r + rr >= i) break;
\r
669 //rr = popcount(x >> (D-8));
\r
670 rr = popCount[x >> (D-8)];
\r
671 //rr = popcount8(x >> (D-8));
\r
672 if (r + rr >= i) break;
\r
677 p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];
\r
682 int darray_pat_construct(darray *da, int n, pb *buf, int k, pb pat, int opt)
\r
686 mymalloc(b,(n+D-1)/D,0);
\r
688 for (i=0; i<n-k+1; i++) {
\r
689 setbit(b,i,getpattern(buf,i,k,pat));
\r
691 for (i=n-k+1; i<n; i++) {
\r
695 darray_construct(da,n,b,opt);
\r