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];
162 static int built = 0;
164 void make___selecttbl(void) {
170 for (x = 0; x < 256; x++) {
171 __setbits(buf,0,8,x);
172 for (r=0; r<8; r++) __selecttbl[(r<<8)+x] = -1;
174 for (i=0; i<8; i++) {
175 if (__getbit(buf,i)) {
176 __selecttbl[(r<<8)+x] = i;
184 unsigned int __popCount(uint x) {
188 r = r - ((r>>1) & 0x77777777) - ((r>>2) & 0x33333333) - ((r>>3) & 0x11111111);
189 r = ((r + (r>>4)) & 0x0f0f0f0f) % 0xff;
192 r = ((r & 0xaaaaaaaa)>>1) + (r & 0x55555555);
193 r = ((r & 0xcccccccc)>>2) + (r & 0x33333333);
194 //r = ((r & 0xf0f0f0f0)>>4) + (r & 0x0f0f0f0f);
195 r = ((r>>4) + r) & 0x0f0f0f0f;
196 //r = ((r & 0xff00ff00)>>8) + (r & 0x00ff00ff);
198 //r = ((r & 0xffff0000)>>16) + (r & 0x0000ffff);
199 r = ((r>>16) + r) & 63;
201 r = _popCount[x & 0xff];
203 r += _popCount[x & 0xff];
205 r += _popCount[x & 0xff];
207 r += _popCount[x & 0xff];
213 unsigned int __popCount8(uint x) {
217 r = ((r & 0xaa)>>1) + (r & 0x55);
218 r = ((r & 0xcc)>>2) + (r & 0x33);
219 r = ((r>>4) + r) & 0x0f;
221 r = _popCount[x & 0xff];
227 int selectd2_save(selectd2 * s, FILE * fp) {
229 wr += fwrite(&s->n,sizeof(uint),1,fp);
230 wr += fwrite(&s->m,sizeof(uint),1,fp);
231 wr += fwrite(&s->size,sizeof(uint),1,fp);
232 wr += fwrite(&s->ss_len,sizeof(uint),1,fp);
233 wr += fwrite(&s->sl_len,sizeof(uint),1,fp);
234 wr += fwrite(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
235 uint nl = (s->m-1) / L + 1;
236 wr += fwrite(s->lp,sizeof(uint),nl+1,fp);
237 wr += fwrite(s->p,sizeof(uint),nl+1,fp);
238 wr += fwrite(s->ss,sizeof(ushort),s->ss_len,fp);
239 wr += fwrite(s->sl,sizeof(uint),s->sl_len,fp);
240 if(wr!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
246 int selectd2_load(selectd2 * s, FILE * fp) {
248 rd += fread(&s->n,sizeof(uint),1,fp);
249 rd += fread(&s->m,sizeof(uint),1,fp);
250 rd += fread(&s->size,sizeof(uint),1,fp);
251 rd += fread(&s->ss_len,sizeof(uint),1,fp);
252 rd += fread(&s->sl_len,sizeof(uint),1,fp);
253 s->buf = new uchar[(s->n+7)/8+1];
254 rd += fread(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
255 uint nl = (s->m-1) / L + 1;
256 s->lp = new uint[nl+1];
257 rd += fread(s->lp,sizeof(uint),nl+1,fp);
258 s->p = new uint[nl+1];
259 rd += fread(s->p,sizeof(uint),nl+1,fp);
260 s->ss = new ushort[s->ss_len];
261 rd += fread(s->ss,sizeof(ushort),s->ss_len,fp);
262 s->sl = new uint[s->sl_len];
263 rd += fread(s->sl,sizeof(uint),s->sl_len,fp);
264 if(rd!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
270 void selectd2_free(selectd2 * s) {
279 int selectd2_construct(selectd2 *select, int n, uchar *buf) {
290 printf("ERROR: L=%d LLL=%d\n",L,LLL);
295 for (i=0; i<n; i++) m += __getbit2(buf,i);
298 //printf("n=%d m=%d\n",n,m);
304 for (i=0; i<n; i++) {
305 if (__getbit2(buf,i)) {
312 select->size = 0; //ignoring buf, shared with selects3
313 select->lp = new uint[nl+1];
314 for(int k=0;k<nl+1;k++) select->lp[k]=0;
315 select->size += (nl+1)*sizeof(uint);
316 select->p = new uint[nl+1];
317 for(int k=0;k<nl+1;k++) select->p[k]=0;
318 select->size += (nl+1)*sizeof(uint);
320 for (r = 0; r < 2; r++) {
322 for (il = 0; il < nl; il++) {
325 i = min((il+1)*L-1,m-1);
327 //printf("%d ",p-pp);
330 for (is = 0; is < L; is++) {
331 if (il*L+is >= m) break;
332 select->sl[ml*L+is] = s[il*L+is];
335 select->p[il] = -((ml<<logL)+1);
340 for (is = 0; is < L/LLL; is++) {
341 if (il*L+is*LLL >= m) break;
342 select->ss[ms*(L/LLL)+is] = s[il*L+is*LLL] - pp;
345 select->p[il] = ms << (logL-logLLL);
350 select->sl = new uint[ml*L+1];
351 for(int k=0;k<ml*L+1;k++) select->sl[k]=0;
352 select->size += sizeof(uint)*(ml*L+1);
353 select->sl_len = ml*L+1;
354 select->ss = new ushort[ms*(L/LLL)+1];
355 for(int k=0;k<ms*(L/LLL)+1;k++) select->ss[k]=0;
356 select->ss_len = ms*(L/LLL)+1;
357 select->size += sizeof(ushort)*(ms*(L/LLL)+1);
365 int selectd2_select(selectd2 *select, int i,int f) {
371 if (i == 0) return -1;
375 printf("ERROR: m=%d i=%d\n",select->m,i);
382 il = select->p[i>>logL];
385 //p = select->sl[(il<<logL)+(i & (L-1))];
386 p = select->sl[il+(i & (L-1))];
389 p = select->lp[i>>logL];
390 //p += select->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
391 p += select->ss[il+((i & (L-1))>>logLLL)];
392 r = i - (i & (LLL-1));
394 q = &(select->buf[p>>3]);
398 r -= _popCount[*q >> (8-1-rr)];
403 if (r + rr >= i) break;
408 p = (q - select->buf) << 3;
409 p += __selecttbl[((i-r-1)<<8)+(*q)];
413 r -= _popCount[(*q ^ 0xff) >> (8-1-rr)];
417 rr = _popCount[*q ^ 0xff];
418 if (r + rr >= i) break;
423 p = (q - select->buf) << 3;
424 p += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
431 int selectd2_select2(selectd2 *select, int i,int f, int *st, int *en) {
444 printf("ERROR: m=%d i=%d\n",select->m,i);
451 il = select->p[i>>logL];
454 //p = select->sl[(il<<logL)+(i & (L-1))];
455 p = select->sl[il+(i & (L-1))];
457 if ((i>>logL) == ((i+1)>>logL)) {
458 p2 = select->sl[il+((i+1) & (L-1))];
461 p2 = selectd2_select(select,i+2,f);
465 p = select->lp[i>>logL];
466 //p += select->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
467 p += select->ss[il+((i & (L-1))>>logLLL)];
468 r = i - (i & (LLL-1));
470 q = &(select->buf[p>>3]);
474 r -= _popCount[*q >> (8-1-rr)];
479 if (r + rr >= i) break;
484 p = (q - select->buf) << 3;
485 p += __selecttbl[((i-r-1)<<8)+(*q)];
487 if ((i>>logL) == ((i+1)>>logL)) {
491 if (r + rr >= i) break;
495 p2 = (q - select->buf) << 3;
496 p2 += __selecttbl[((i-r-1)<<8)+(*q)];
499 p2 = selectd2_select(select,i+2,f);
505 r -= _popCount[(*q ^ 0xff) >> (8-1-rr)];
509 rr = _popCount[*q ^ 0xff];
510 if (r + rr >= i) break;
515 p = (q - select->buf) << 3;
516 p += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
518 if ((i>>logL) == ((i+1)>>logL)) {
521 rr = _popCount[*q ^ 0xff];
522 if (r + rr >= i) break;
526 p2 = (q - select->buf) << 3;
527 p2 += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
530 p2 = selectd2_select(select,i+2,f);
540 int selects3_save(selects3 * s, FILE * fp) {
542 wr += fwrite(&s->n,sizeof(uint),1,fp);
543 wr += fwrite(&s->m,sizeof(uint),1,fp);
544 wr += fwrite(&s->size,sizeof(uint),1,fp);
545 wr += fwrite(&s->d,sizeof(uint),1,fp);
546 wr += fwrite(&s->hi_len,sizeof(uint),1,fp);
547 wr += fwrite(&s->low_len,sizeof(uint),1,fp);
548 wr += fwrite(s->hi,sizeof(uchar),s->hi_len,fp);
549 wr += fwrite(s->low,sizeof(uint),s->low_len,fp);
550 if(wr!=(6+s->hi_len+s->low_len))
552 if(selectd2_save(s->sd0,fp)) return 2;
553 if(selectd2_save(s->sd1,fp)) return 3;
558 int selects3_load(selects3 * s, FILE * fp) {
560 rd += fread(&s->n,sizeof(uint),1,fp);
561 rd += fread(&s->m,sizeof(uint),1,fp);
562 rd += fread(&s->size,sizeof(uint),1,fp);
563 rd += fread(&s->d,sizeof(uint),1,fp);
564 rd += fread(&s->hi_len,sizeof(uint),1,fp);
565 rd += fread(&s->low_len,sizeof(uint),1,fp);
566 s->hi = new uchar[s->hi_len];
567 rd += fread(s->hi,sizeof(uchar),s->hi_len,fp);
568 s->low = new uint[s->low_len];
569 rd += fread(s->low,sizeof(uint),s->low_len,fp);
570 if(rd!=(6+s->hi_len+s->low_len))
572 s->sd0 = new selectd2;
573 if(selectd2_load(s->sd0,fp)) return 2;
574 s->sd1 = new selectd2;
575 if(selectd2_load(s->sd1,fp)) return 3;
576 delete [] s->sd0->buf;
577 delete [] s->sd1->buf;
584 void selects3_free(selects3 * s) {
587 //delete [] s->sd0->buf;
588 selectd2_free(s->sd0);
590 selectd2_free(s->sd1);
595 int selects3_construct(selects3 *select, int n, uint *buf) {
603 for (i=0; i<n; i++) m += __getbit(buf,i);
607 if (m == 0) return 0;
618 buf2 = new uchar[(2*m+8-1)/8+1];
619 for(int k=0;k<(2*m+8-1)/8+1;k++) buf2[k]=0;
620 select->hi_len = (2*m+8-1)/8+1;
621 low = new uint[(d*m+PBS-1)/PBS+1];
622 for(uint k=0;k<(d*m+PBS-1)/PBS+1;k++) low[k]=0;
623 select->low_len = (d*m+PBS-1)/PBS+1;
627 select->size = sizeof(uchar)*((2*m+8-1)/8+1) + sizeof(uint)*((d*m+PBS-1)/PBS+1);
629 for (i=0; i<m*2; i++) __setbit2(buf2,i,0);
632 for (i=0; i<n; i++) {
633 if (__getbit(buf,i)) {
634 __setbit2(buf2,(i>>d)+m,1);
635 __setbits(low,m*d,d,i & ((1<<d)-1));
642 select->size += 2*sizeof(selectd2);
644 selectd2_construct(sd1,m*2,buf2);
647 for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
648 selectd2_construct(sd0,m*2,buf2);
651 for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
655 //selects3 * lasts3=NULL;
659 int selects3_select(selects3 *select, int i) {
664 printf("ERROR: m=%d i=%d\n",select->m,i);
669 if (i == 0) return -1;
672 /*if(select->lasti==(uint)i-1) {
673 while(!__getbit2(select->sd1->buf,++select->lasts));
676 select->lasts = selectd2_select(select->sd1,i,1);
679 //lasts3 = select; */
680 x = selectd2_select(select->sd1,i,1) - (i-1);
681 //x = (select->lasts-(i-1)) << d;
683 x += __getbits(select->low,(i-1)*d,d);
688 int selects3_selectnext(selects3 *select, int i) {
689 //return selects3_select(select,selects3_rank(select,i)+1);
697 y = selectd2_select(select->sd0,ii,0)+1;
706 if (((z << r) & 0x80) == 0) {
710 w = __getbits(q,x*d,d);
713 if(__getbit2(select->hi,(8*y+r))) k2++;
721 if(__getbit2(select->hi,(8*y+r))) k2++;
732 for(int kk=0;kk<8-r;kk++) {
733 if(__getbit2(select->hi,c)) {
741 while(select->hi[pp]==0) {
745 while(!__getbit2(select->hi,c)) c++;
748 return __getbits(q,x*d,d)+((c)<<d);
751 int selects3_rank(selects3 *select, int i) {
762 y = selectd2_select(select->sd0,ii,0)+1;
763 // selectd2_select2(select->sd0,ii,0,&y1,&y2);
765 //printf("y %d y1 %d %d\n",y,y1,y2-y1);
775 if (((z << r) & 0x80) == 0) break;
776 w = __getbits(q,x*d,d);