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;
181 unsigned int __popCount(uint x) {
185 r = r - ((r>>1) & 0x77777777) - ((r>>2) & 0x33333333) - ((r>>3) & 0x11111111);
186 r = ((r + (r>>4)) & 0x0f0f0f0f) % 0xff;
189 r = ((r & 0xaaaaaaaa)>>1) + (r & 0x55555555);
190 r = ((r & 0xcccccccc)>>2) + (r & 0x33333333);
191 //r = ((r & 0xf0f0f0f0)>>4) + (r & 0x0f0f0f0f);
192 r = ((r>>4) + r) & 0x0f0f0f0f;
193 //r = ((r & 0xff00ff00)>>8) + (r & 0x00ff00ff);
195 //r = ((r & 0xffff0000)>>16) + (r & 0x0000ffff);
196 r = ((r>>16) + r) & 63;
198 r = _popCount[x & 0xff];
200 r += _popCount[x & 0xff];
202 r += _popCount[x & 0xff];
204 r += _popCount[x & 0xff];
210 unsigned int __popCount8(uint x) {
214 r = ((r & 0xaa)>>1) + (r & 0x55);
215 r = ((r & 0xcc)>>2) + (r & 0x33);
216 r = ((r>>4) + r) & 0x0f;
218 r = _popCount[x & 0xff];
224 int selectd2_save(selectd2 * s, FILE * fp) {
226 wr += fwrite(&s->n,sizeof(uint),1,fp);
227 wr += fwrite(&s->m,sizeof(uint),1,fp);
228 wr += fwrite(&s->size,sizeof(uint),1,fp);
229 wr += fwrite(&s->ss_len,sizeof(uint),1,fp);
230 wr += fwrite(&s->sl_len,sizeof(uint),1,fp);
231 wr += fwrite(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
232 uint nl = (s->m-1) / L + 1;
233 wr += fwrite(s->lp,sizeof(uint),nl+1,fp);
234 wr += fwrite(s->p,sizeof(uint),nl+1,fp);
235 wr += fwrite(s->ss,sizeof(ushort),s->ss_len,fp);
236 wr += fwrite(s->sl,sizeof(uint),s->sl_len,fp);
237 if(wr!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
243 int selectd2_load(selectd2 * s, FILE * fp) {
245 rd += fread(&s->n,sizeof(uint),1,fp);
246 rd += fread(&s->m,sizeof(uint),1,fp);
247 rd += fread(&s->size,sizeof(uint),1,fp);
248 rd += fread(&s->ss_len,sizeof(uint),1,fp);
249 rd += fread(&s->sl_len,sizeof(uint),1,fp);
250 s->buf = new uchar[(s->n+7)/8+1];
251 rd += fread(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
252 uint nl = (s->m-1) / L + 1;
253 s->lp = new uint[nl+1];
254 rd += fread(s->lp,sizeof(uint),nl+1,fp);
255 s->p = new uint[nl+1];
256 rd += fread(s->p,sizeof(uint),nl+1,fp);
257 s->ss = new ushort[s->ss_len];
258 rd += fread(s->ss,sizeof(ushort),s->ss_len,fp);
259 s->sl = new uint[s->sl_len];
260 rd += fread(s->sl,sizeof(uint),s->sl_len,fp);
261 if(rd!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
267 void selectd2_free(selectd2 * s) {
276 int selectd2_construct(selectd2 *select, int n, uchar *buf) {
287 printf("ERROR: L=%d LLL=%d\n",L,LLL);
292 for (i=0; i<n; i++) m += __getbit2(buf,i);
295 //printf("n=%d m=%d\n",n,m);
301 for (i=0; i<n; i++) {
302 if (__getbit2(buf,i)) {
309 select->size = 0; //ignoring buf, shared with selects3
310 select->lp = new uint[nl+1];
311 for(int k=0;k<nl+1;k++) select->lp[k]=0;
312 select->size += (nl+1)*sizeof(uint);
313 select->p = new uint[nl+1];
314 for(int k=0;k<nl+1;k++) select->p[k]=0;
315 select->size += (nl+1)*sizeof(uint);
317 for (r = 0; r < 2; r++) {
319 for (il = 0; il < nl; il++) {
322 i = min((il+1)*L-1,m-1);
324 //printf("%d ",p-pp);
327 for (is = 0; is < L; is++) {
328 if (il*L+is >= m) break;
329 select->sl[ml*L+is] = s[il*L+is];
332 select->p[il] = -((ml<<logL)+1);
337 for (is = 0; is < L/LLL; is++) {
338 if (il*L+is*LLL >= m) break;
339 select->ss[ms*(L/LLL)+is] = s[il*L+is*LLL] - pp;
342 select->p[il] = ms << (logL-logLLL);
347 select->sl = new uint[ml*L+1];
348 for(int k=0;k<ml*L+1;k++) select->sl[k]=0;
349 select->size += sizeof(uint)*(ml*L+1);
350 select->sl_len = ml*L+1;
351 select->ss = new ushort[ms*(L/LLL)+1];
352 for(int k=0;k<ms*(L/LLL)+1;k++) select->ss[k]=0;
353 select->ss_len = ms*(L/LLL)+1;
354 select->size += sizeof(ushort)*(ms*(L/LLL)+1);
362 int selectd2_select(selectd2 *select, int i,int f) {
368 if (i == 0) return -1;
372 printf("ERROR: m=%d i=%d\n",select->m,i);
379 il = select->p[i>>logL];
382 //p = select->sl[(il<<logL)+(i & (L-1))];
383 p = select->sl[il+(i & (L-1))];
386 p = select->lp[i>>logL];
387 //p += select->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
388 p += select->ss[il+((i & (L-1))>>logLLL)];
389 r = i - (i & (LLL-1));
391 q = &(select->buf[p>>3]);
395 r -= _popCount[*q >> (8-1-rr)];
400 if (r + rr >= i) break;
405 p = (q - select->buf) << 3;
406 p += __selecttbl[((i-r-1)<<8)+(*q)];
410 r -= _popCount[(*q ^ 0xff) >> (8-1-rr)];
414 rr = _popCount[*q ^ 0xff];
415 if (r + rr >= i) break;
420 p = (q - select->buf) << 3;
421 p += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
428 int selectd2_select2(selectd2 *select, int i,int f, int *st, int *en) {
441 printf("ERROR: m=%d i=%d\n",select->m,i);
448 il = select->p[i>>logL];
451 //p = select->sl[(il<<logL)+(i & (L-1))];
452 p = select->sl[il+(i & (L-1))];
454 if ((i>>logL) == ((i+1)>>logL)) {
455 p2 = select->sl[il+((i+1) & (L-1))];
458 p2 = selectd2_select(select,i+2,f);
462 p = select->lp[i>>logL];
463 //p += select->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
464 p += select->ss[il+((i & (L-1))>>logLLL)];
465 r = i - (i & (LLL-1));
467 q = &(select->buf[p>>3]);
471 r -= _popCount[*q >> (8-1-rr)];
476 if (r + rr >= i) break;
481 p = (q - select->buf) << 3;
482 p += __selecttbl[((i-r-1)<<8)+(*q)];
484 if ((i>>logL) == ((i+1)>>logL)) {
488 if (r + rr >= i) break;
492 p2 = (q - select->buf) << 3;
493 p2 += __selecttbl[((i-r-1)<<8)+(*q)];
496 p2 = selectd2_select(select,i+2,f);
502 r -= _popCount[(*q ^ 0xff) >> (8-1-rr)];
506 rr = _popCount[*q ^ 0xff];
507 if (r + rr >= i) break;
512 p = (q - select->buf) << 3;
513 p += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
515 if ((i>>logL) == ((i+1)>>logL)) {
518 rr = _popCount[*q ^ 0xff];
519 if (r + rr >= i) break;
523 p2 = (q - select->buf) << 3;
524 p2 += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
527 p2 = selectd2_select(select,i+2,f);
537 int selects3_save(selects3 * s, FILE * fp) {
539 wr += fwrite(&s->n,sizeof(uint),1,fp);
540 wr += fwrite(&s->m,sizeof(uint),1,fp);
541 wr += fwrite(&s->size,sizeof(uint),1,fp);
542 wr += fwrite(&s->d,sizeof(uint),1,fp);
543 wr += fwrite(&s->hi_len,sizeof(uint),1,fp);
544 wr += fwrite(&s->low_len,sizeof(uint),1,fp);
545 wr += fwrite(s->hi,sizeof(uchar),s->hi_len,fp);
546 wr += fwrite(s->low,sizeof(uint),s->low_len,fp);
547 if(wr!=(6+s->hi_len+s->low_len))
549 if(selectd2_save(s->sd0,fp)) return 2;
550 if(selectd2_save(s->sd1,fp)) return 3;
555 int selects3_load(selects3 * s, FILE * fp) {
557 rd += fread(&s->n,sizeof(uint),1,fp);
558 rd += fread(&s->m,sizeof(uint),1,fp);
559 rd += fread(&s->size,sizeof(uint),1,fp);
560 rd += fread(&s->d,sizeof(uint),1,fp);
561 rd += fread(&s->hi_len,sizeof(uint),1,fp);
562 rd += fread(&s->low_len,sizeof(uint),1,fp);
563 s->hi = new uchar[s->hi_len];
564 rd += fread(s->hi,sizeof(uchar),s->hi_len,fp);
565 s->low = new uint[s->low_len];
566 rd += fread(s->low,sizeof(uint),s->low_len,fp);
567 if(rd!=(6+s->hi_len+s->low_len))
569 s->sd0 = new selectd2;
570 if(selectd2_load(s->sd0,fp)) return 2;
571 s->sd1 = new selectd2;
572 if(selectd2_load(s->sd1,fp)) return 3;
573 delete [] s->sd0->buf;
574 delete [] s->sd1->buf;
581 void selects3_free(selects3 * s) {
584 //delete [] s->sd0->buf;
585 selectd2_free(s->sd0);
587 selectd2_free(s->sd1);
592 int selects3_construct(selects3 *select, int n, uint *buf) {
600 for (i=0; i<n; i++) m += __getbit(buf,i);
604 if (m == 0) return 0;
615 buf2 = new uchar[(2*m+8-1)/8+1];
616 for(int k=0;k<(2*m+8-1)/8+1;k++) buf2[k]=0;
617 select->hi_len = (2*m+8-1)/8+1;
618 low = new uint[(d*m+PBS-1)/PBS+1];
619 for(uint k=0;k<(d*m+PBS-1)/PBS+1;k++) low[k]=0;
620 select->low_len = (d*m+PBS-1)/PBS+1;
624 select->size = sizeof(uchar)*((2*m+8-1)/8+1) + sizeof(uint)*((d*m+PBS-1)/PBS+1);
626 for (i=0; i<m*2; i++) __setbit2(buf2,i,0);
629 for (i=0; i<n; i++) {
630 if (__getbit(buf,i)) {
631 __setbit2(buf2,(i>>d)+m,1);
632 __setbits(low,m*d,d,i & ((1<<d)-1));
639 select->size += 2*sizeof(selectd2);
641 selectd2_construct(sd1,m*2,buf2);
644 for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
645 selectd2_construct(sd0,m*2,buf2);
648 for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
652 //selects3 * lasts3=NULL;
656 int selects3_select(selects3 *select, int i) {
661 printf("ERROR: m=%d i=%d\n",select->m,i);
666 if (i == 0) return -1;
669 if(select->lasti==(uint)i-1) {
670 while(!__getbit2(select->sd1->buf,++select->lasti));
673 select->lasts = selectd2_select(select->sd1,i,1);
677 x = (select->lasts-(i-1)) << d;
678 x += __getbits(select->low,(i-1)*d,d);
683 int selects3_selectnext(selects3 *select, int i) {
684 return selects3_select(select,selects3_rank(select,i)+1);
692 y = selectd2_select(select->sd0,ii,0)+1;
701 if (((z << r) & 0x80) == 0) {
705 w = __getbits(q,x*d,d);
708 if(__getbit2(select->hi,(8*y+r))) k2++;
716 if(__getbit2(select->hi,(8*y+r))) k2++;
727 for(int kk=0;kk<8-r;kk++) {
728 if(__getbit2(select->hi,c)) {
736 while(select->hi[pp]==0) {
740 while(!__getbit2(select->hi,c)) c++;
743 return __getbits(q,x*d,d)+((c)<<d);
746 int selects3_rank(selects3 *select, int i) {
756 y = selectd2_select(select->sd0,ii,0)+1;
757 // selectd2_select2(select->sd0,ii,0,&y1,&y2);
759 //printf("y %d y1 %d %d\n",y,y1,y2-y1);
769 if (((z << r) & 0x80) == 0) break;
770 w = __getbits(q,x*d,d);