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++) {
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++) {
169 for (r=0; r<8; r++) selecttbl[(r<<8)+x] = -1;
171 for (i=0; i<8; 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++) {
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);