5 typedef unsigned int qword;
8 typedef unsigned long long qword;
11 #define PBS (sizeof(uint)*8)
17 #define logLL 16 // size of word
23 #define LLL (1<<logLLL)
25 //#define logL (logLL-3)
26 #define logL (logLL-1-5)
40 int __setbit(uint *B, int i,int x) {
45 if (x==0) B[j] &= (~(1<<(D-1-l)));
46 else if (x==1) B[j] |= (1<<(D-1-l));
48 printf("error __setbit x=%d\n",x);
55 int __setbit2(uchar *B, int i,int x) {
60 if (x==0) B[j] &= (~(1<<(8-1-l)));
61 else if (x==1) B[j] |= (1<<(8-1-l));
63 printf("error __setbit2 x=%d\n",x);
70 int __setbits(uint *B, int i, int d, int x) {
74 __setbit(B,i+j,(x>>(d-j-1))&1);
80 int __getbit(uint *B, int i) {
87 return (B[j] >> (D-1-l)) & 1;
91 int __getbit2(uchar *B, int i) {
98 return (B[j] >> (8-1-l)) & 1;
103 uint __getbits(uint *B, int i, int d) {
109 x = (((qword)B[0]) << D) + B[1];
115 x = (((qword)B[0])<<D)+B[1];
118 x &= (((qword)1L<<D)-1)<<D;
130 uint __getbits(uint *B, int i, int d) {
134 for (j=0; j<d; j++) {
136 x += __getbit(B,i+j);
142 static const unsigned int _popCount[] = {
143 0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
144 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
145 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
146 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
147 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
148 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
149 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
150 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
151 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
152 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
153 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
154 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
155 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
156 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
157 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
158 4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8
161 static unsigned int __selecttbl[8*256];
163 void make___selecttbl(void) {
167 for (x = 0; x < 256; x++) {
168 __setbits(buf,0,8,x);
169 for (r=0; r<8; r++) __selecttbl[(r<<8)+x] = -1;
171 for (i=0; i<8; i++) {
172 if (__getbit(buf,i)) {
173 __selecttbl[(r<<8)+x] = i;
180 unsigned int __popCount(uint x) {
184 r = r - ((r>>1) & 0x77777777) - ((r>>2) & 0x33333333) - ((r>>3) & 0x11111111);
185 r = ((r + (r>>4)) & 0x0f0f0f0f) % 0xff;
188 r = ((r & 0xaaaaaaaa)>>1) + (r & 0x55555555);
189 r = ((r & 0xcccccccc)>>2) + (r & 0x33333333);
190 //r = ((r & 0xf0f0f0f0)>>4) + (r & 0x0f0f0f0f);
191 r = ((r>>4) + r) & 0x0f0f0f0f;
192 //r = ((r & 0xff00ff00)>>8) + (r & 0x00ff00ff);
194 //r = ((r & 0xffff0000)>>16) + (r & 0x0000ffff);
195 r = ((r>>16) + r) & 63;
197 r = _popCount[x & 0xff];
199 r += _popCount[x & 0xff];
201 r += _popCount[x & 0xff];
203 r += _popCount[x & 0xff];
209 unsigned int __popCount8(uint x) {
213 r = ((r & 0xaa)>>1) + (r & 0x55);
214 r = ((r & 0xcc)>>2) + (r & 0x33);
215 r = ((r>>4) + r) & 0x0f;
217 r = _popCount[x & 0xff];
222 int selectd2_save(selectd2 * s, FILE * fp) {
224 wr += fwrite(&s->n,sizeof(uint),1,fp);
225 wr += fwrite(&s->m,sizeof(uint),1,fp);
226 wr += fwrite(&s->size,sizeof(uint),1,fp);
227 wr += fwrite(&s->ss_len,sizeof(uint),1,fp);
228 wr += fwrite(&s->sl_len,sizeof(uint),1,fp);
229 wr += fwrite(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
230 uint nl = (s->m-1) / L + 1;
231 wr += fwrite(s->lp,sizeof(uint),nl+1,fp);
232 wr += fwrite(s->p,sizeof(uint),nl+1,fp);
233 wr += fwrite(s->ss,sizeof(ushort),s->ss_len,fp);
234 wr += fwrite(s->sl,sizeof(uint),s->sl_len,fp);
235 if(wr!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
240 int selectd2_load(selectd2 * s, FILE * fp) {
242 rd += fread(&s->n,sizeof(uint),1,fp);
243 rd += fread(&s->m,sizeof(uint),1,fp);
244 rd += fread(&s->size,sizeof(uint),1,fp);
245 rd += fread(&s->ss_len,sizeof(uint),1,fp);
246 rd += fread(&s->sl_len,sizeof(uint),1,fp);
247 s->buf = new uchar[(s->n+7)/8+1];
248 rd += fread(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
249 uint nl = (s->m-1) / L + 1;
250 s->lp = new uint[nl+1];
251 rd += fread(s->lp,sizeof(uint),nl+1,fp);
252 s->p = new uint[nl+1];
253 rd += fread(s->p,sizeof(uint),nl+1,fp);
254 s->ss = new ushort[s->ss_len];
255 rd += fread(s->ss,sizeof(ushort),s->ss_len,fp);
256 s->sl = new uint[s->sl_len];
257 rd += fread(s->sl,sizeof(uint),s->sl_len,fp);
258 if(rd!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
263 void selectd2_free(selectd2 * s) {
271 int selectd2_construct(selectd2 *select, int n, uchar *buf) {
282 printf("ERROR: L=%d LLL=%d\n",L,LLL);
287 for (i=0; i<n; i++) m += __getbit2(buf,i);
290 //printf("n=%d m=%d\n",n,m);
296 for (i=0; i<n; i++) {
297 if (__getbit2(buf,i)) {
304 select->size = 0; //ignoring buf, shared with selects3
305 select->lp = new uint[nl+1];
306 for(int k=0;k<nl+1;k++) select->lp[k]=0;
307 select->size += (nl+1)*sizeof(uint);
308 select->p = new uint[nl+1];
309 for(int k=0;k<nl+1;k++) select->p[k]=0;
310 select->size += (nl+1)*sizeof(uint);
312 for (r = 0; r < 2; r++) {
314 for (il = 0; il < nl; il++) {
317 i = min((il+1)*L-1,m-1);
319 //printf("%d ",p-pp);
322 for (is = 0; is < L; is++) {
323 if (il*L+is >= m) break;
324 select->sl[ml*L+is] = s[il*L+is];
327 select->p[il] = -((ml<<logL)+1);
332 for (is = 0; is < L/LLL; is++) {
333 if (il*L+is*LLL >= m) break;
334 select->ss[ms*(L/LLL)+is] = s[il*L+is*LLL] - pp;
337 select->p[il] = ms << (logL-logLLL);
342 select->sl = new uint[ml*L+1];
343 for(int k=0;k<ml*L+1;k++) select->sl[k]=0;
344 select->size += sizeof(uint)*(ml*L+1);
345 select->sl_len = ml*L+1;
346 select->ss = new ushort[ms*(L/LLL)+1];
347 for(int k=0;k<ms*(L/LLL)+1;k++) select->ss[k]=0;
348 select->ss_len = ms*(L/LLL)+1;
349 select->size += sizeof(ushort)*(ms*(L/LLL)+1);
357 int selectd2_select(selectd2 *select, int i,int f) {
363 if (i == 0) return -1;
367 printf("ERROR: m=%d i=%d\n",select->m,i);
374 il = select->p[i>>logL];
377 //p = select->sl[(il<<logL)+(i & (L-1))];
378 p = select->sl[il+(i & (L-1))];
381 p = select->lp[i>>logL];
382 //p += select->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
383 p += select->ss[il+((i & (L-1))>>logLLL)];
384 r = i - (i & (LLL-1));
386 q = &(select->buf[p>>3]);
390 r -= _popCount[*q >> (8-1-rr)];
395 if (r + rr >= i) break;
400 p = (q - select->buf) << 3;
401 p += __selecttbl[((i-r-1)<<8)+(*q)];
405 r -= _popCount[(*q ^ 0xff) >> (8-1-rr)];
409 rr = _popCount[*q ^ 0xff];
410 if (r + rr >= i) break;
415 p = (q - select->buf) << 3;
416 p += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
423 int selectd2_select2(selectd2 *select, int i,int f, int *st, int *en) {
436 printf("ERROR: m=%d i=%d\n",select->m,i);
443 il = select->p[i>>logL];
446 //p = select->sl[(il<<logL)+(i & (L-1))];
447 p = select->sl[il+(i & (L-1))];
449 if ((i>>logL) == ((i+1)>>logL)) {
450 p2 = select->sl[il+((i+1) & (L-1))];
453 p2 = selectd2_select(select,i+2,f);
457 p = select->lp[i>>logL];
458 //p += select->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
459 p += select->ss[il+((i & (L-1))>>logLLL)];
460 r = i - (i & (LLL-1));
462 q = &(select->buf[p>>3]);
466 r -= _popCount[*q >> (8-1-rr)];
471 if (r + rr >= i) break;
476 p = (q - select->buf) << 3;
477 p += __selecttbl[((i-r-1)<<8)+(*q)];
479 if ((i>>logL) == ((i+1)>>logL)) {
483 if (r + rr >= i) break;
487 p2 = (q - select->buf) << 3;
488 p2 += __selecttbl[((i-r-1)<<8)+(*q)];
491 p2 = selectd2_select(select,i+2,f);
497 r -= _popCount[(*q ^ 0xff) >> (8-1-rr)];
501 rr = _popCount[*q ^ 0xff];
502 if (r + rr >= i) break;
507 p = (q - select->buf) << 3;
508 p += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
510 if ((i>>logL) == ((i+1)>>logL)) {
513 rr = _popCount[*q ^ 0xff];
514 if (r + rr >= i) break;
518 p2 = (q - select->buf) << 3;
519 p2 += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
522 p2 = selectd2_select(select,i+2,f);
532 int selects3_save(selects3 * s, FILE * fp) {
534 wr += fwrite(&s->n,sizeof(uint),1,fp);
535 wr += fwrite(&s->m,sizeof(uint),1,fp);
536 wr += fwrite(&s->size,sizeof(uint),1,fp);
537 wr += fwrite(&s->d,sizeof(uint),1,fp);
538 wr += fwrite(&s->hi_len,sizeof(uint),1,fp);
539 wr += fwrite(&s->low_len,sizeof(uint),1,fp);
540 wr += fwrite(s->hi,sizeof(uchar),s->hi_len,fp);
541 wr += fwrite(s->low,sizeof(uint),s->low_len,fp);
542 if(wr!=(6+s->hi_len+s->low_len))
544 if(selectd2_save(s->sd0,fp)) return 2;
545 if(selectd2_save(s->sd1,fp)) return 3;
549 int selects3_load(selects3 * s, FILE * fp) {
551 rd += fread(&s->n,sizeof(uint),1,fp);
552 rd += fread(&s->m,sizeof(uint),1,fp);
553 rd += fread(&s->size,sizeof(uint),1,fp);
554 rd += fread(&s->d,sizeof(uint),1,fp);
555 rd += fread(&s->hi_len,sizeof(uint),1,fp);
556 rd += fread(&s->low_len,sizeof(uint),1,fp);
557 s->hi = new uchar[s->hi_len];
558 rd += fread(s->hi,sizeof(uchar),s->hi_len,fp);
559 s->low = new uint[s->low_len];
560 rd += fread(s->low,sizeof(uint),s->low_len,fp);
561 if(rd!=(6+s->hi_len+s->low_len))
563 s->sd0 = new selectd2;
564 if(selectd2_load(s->sd0,fp)) return 2;
565 s->sd1 = new selectd2;
566 if(selectd2_load(s->sd1,fp)) return 3;
567 delete [] s->sd0->buf;
568 delete [] s->sd1->buf;
574 void selects3_free(selects3 * s) {
577 //delete [] s->sd0->buf;
578 selectd2_free(s->sd0);
580 selectd2_free(s->sd1);
584 int selects3_construct(selects3 *select, int n, uint *buf) {
592 for (i=0; i<n; i++) m += __getbit(buf,i);
596 if (m == 0) return 0;
607 buf2 = new uchar[(2*m+8-1)/8+1];
608 for(int k=0;k<(2*m+8-1)/8+1;k++) buf2[k]=0;
609 select->hi_len = (2*m+8-1)/8+1;
610 low = new uint[(d*m+PBS-1)/PBS+1];
611 for(uint k=0;k<(d*m+PBS-1)/PBS+1;k++) low[k]=0;
612 select->low_len = (d*m+PBS-1)/PBS+1;
616 select->size = sizeof(uchar)*((2*m+8-1)/8+1) + sizeof(uint)*((d*m+PBS-1)/PBS+1);
618 for (i=0; i<m*2; i++) __setbit2(buf2,i,0);
621 for (i=0; i<n; i++) {
622 if (__getbit(buf,i)) {
623 __setbit2(buf2,(i>>d)+m,1);
624 __setbits(low,m*d,d,i & ((1<<d)-1));
631 select->size += 2*sizeof(selectd2);
633 selectd2_construct(sd1,m*2,buf2);
636 for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
637 selectd2_construct(sd0,m*2,buf2);
640 for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
645 int selects3_select(selects3 *select, int i) {
650 printf("ERROR: m=%d i=%d\n",select->m,i);
655 if (i == 0) return -1;
659 x = selectd2_select(select->sd1,i,1) - (i-1);
661 x += __getbits(select->low,(i-1)*d,d);
666 int selects3_rank(selects3 *select, int i) {
676 y = selectd2_select(select->sd0,ii,0)+1;
677 // selectd2_select2(select->sd0,ii,0,&y1,&y2);
679 //printf("y %d y1 %d %d\n",y,y1,y2-y1);
689 if (((z << r) & 0x80) == 0) break;
690 w = __getbits(q,x*d,d);