src/static_bitsequence/sdarray.cpp

   1 #include <sdarray.h>
   2 using std::min;
   3 using std::max;
   4 #if 0
   5 typedef unsigned int qword;
   6 #define logD 4
   7 #else
   8 typedef unsigned long long qword;
   9 #define logD 5
  10 #endif
  11 #define PBS (sizeof(uint)*8)
  12 #define D (1<<logD)
  13 #define logM 5
  14 #define M (1<<logM)
  15 #define logP 8
  16 #define P (1<<logP)
  17 #define logLL 16                 // size of word
  18 #define LL (1<<logLL)
  19 //#define logLLL 7
  20 #define logLLL 5
  21 //#define LLL 128
  22 //#define LLL 32
  23 #define LLL (1<<logLLL)
  24 //#define logL 10
  25 //#define logL (logLL-3)
  26 #define logL (logLL-1-5)
  27 #define L (1<<logL)
  28
  29 int __blog(int x) {
  30   int l;
  31   l = 0;
  32   while (x>0) {
  33     x>>=1;
  34     l++;
  35   }
  36   return l;
  37 }
  38
  39
  40 int __setbit(uint *B, int i,int x) {
  41   int j,l;
  42   //printf("%u\n",D);
  43   j = i / D;
  44   l = i % D;
  45   if (x==0) B[j] &= (~(1<<(D-1-l)));
  46   else if (x==1) B[j] |= (1<<(D-1-l));
  47   else {
  48     printf("error __setbit x=%d\n",x);
  49     exit(1);
  50   }
  51   return x;
  52 }
  53
  54
  55 int __setbit2(uchar *B, int i,int x) {
  56   int j,l;
  57
  58   j = i / 8;
  59   l = i % 8;
  60   if (x==0) B[j] &= (~(1<<(8-1-l)));
  61   else if (x==1) B[j] |= (1<<(8-1-l));
  62   else {
  63     printf("error __setbit2 x=%d\n",x);
  64     exit(1);
  65   }
  66   return x;
  67 }
  68
  69
  70 int __setbits(uint *B, int i, int d, int x) {
  71   int j;
  72
  73   for (j=0; j<d; j++) {
  74     __setbit(B,i+j,(x>>(d-j-1))&1);
  75   }
  76   return x;
  77 }
  78
  79
  80 int __getbit(uint *B, int i) {
  81   int j,l;
  82
  83   //j = i / D;
  84   //l = i % D;
  85   j = i >> logD;
  86   l = i & (D-1);
  87   return (B[j] >> (D-1-l)) & 1;
  88 }
  89
  90
  91 static int __getbit2(uchar *B, int i) {
  92   int j,l;
  93
  94   //j = i / D;
  95   //l = i % D;
  96   j = i >> 3;
  97   l = i & (8-1);
  98   return (B[j] >> (8-1-l)) & 1;
  99 }
 100
 101
 102 #if 1
 103 uint __getbits_aux(uint *B, int i, int d) {
 104   qword x,z;
 105   int j = i;
 106   B += (i >> logD);
 107   i &= (D-1);
 108   if (d==8 && (j & 7) == 0) {
 109     i = (24 - i) >> 3;
 110     x = (uint) (((unsigned char*) B)[i]);
 111   } else if (i+d <= 2*D) {
 112     x = (((qword)B[0]) << D) + B[1];
 113     x <<= i;
 114     x >>= (D*2-1-d);
 115     x >>= 1;
 116   }
 117   else {
 118     fprintf(stderr, "Warning: %d, %d\n", D, d);
 119     x = (((qword)B[0])<<D)+B[1];
 120     z = (x<<D)+B[2];
 121     x <<= i;
 122     x &= (((qword)1L<<D)-1)<<D;
 123     z <<= i;
 124     z >>= D;
 125     x += z;
 126     x >>= (2*D-d);
 127   }
 128
 129   return x;
 130 }
 131
 132 uint __getbits(uint *B, int i, int d)
 133 {
 134   ulong x;
 135 //  uint y;
 136   // y = __getbits_aux(B, i,d);
 137   B += (i >> logD);
 138   i &= (D-1);
 139   x = ((ulong *) B)[0];
 140   x = (x << 32)|(x >> 32);
 141   x <<= i;
 142   x >>= 2*D - d;
 143 //  fprintf(stderr, "slow: %i, fast: %i\n",
 144 //        y, (uint) x);
 145   return x;
 146 }
 147
 148 #endif
 149
 150 #if 0
 151 uint __getbits(uint *B, int i, int d) {
 152   uint j,x;
 153
 154   x = 0;
 155   for (j=0; j<d; j++) {
 156     x <<= 1;
 157     x += __getbit(B,i+j);
 158   }
 159   return x;
 160 }
 161 #endif
 162
 163
 164
 165 #ifdef HAS_NATIVE_POPCOUNT
 166 static inline unsigned int popcount(unsigned int n){
 167   asm ("popcnt %1, %0" : "=r" (n) : "0" (n));
 168   return n;
 169 }
 170
 171 static inline unsigned int popcount8(unsigned int n) {
 172   return popcount(n & 0xff);
 173 }
 174
 175 static inline unsigned int _fast_popcount(int x)
 176 {
 177   return popcount8(x);
 178 }
 179
 180
 181 #else
 182
 183 static const unsigned int _popCount[] = {
 184   0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
 185   1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
 186   1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
 187   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 188   1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
 189   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 190   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 191   3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
 192   1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
 193   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 194   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 195   3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
 196   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 197   3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
 198   3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
 199   4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8
 200 };
 201
 202 static inline unsigned int
 203 _fast_popcount(int x) {
 204   return _popCount[x];
 205 }
 206 #endif
 207
 208
 209
 210 static unsigned int __selecttbl[8*256];
 211 static int built = 0;
 212
 213 void make___selecttbl(void) {
 214   if(built) return;
 215   built = 1;
 216   int i,x,r;
 217   uint buf[1] = { 0 };
 218
 219   for (x = 0; x < 256; x++) {
 220     __setbits(buf,0,8,x);
 221     for (r=0; r<8; r++) __selecttbl[(r<<8)+x] = -1;
 222     r = 0;
 223     for (i=0; i<8; i++) {
 224       if (__getbit(buf,i)) {
 225         __selecttbl[(r<<8)+x] = i;
 226         r++;
 227       }
 228     }
 229   }
 230 }
 231
 232
 233 unsigned int __popCount(uint x) {
 234   uint r;
 235   #if 0
 236   r = x;
 237   r = r - ((r>>1) & 0x77777777) - ((r>>2) & 0x33333333) - ((r>>3) & 0x11111111);
 238   r = ((r + (r>>4)) & 0x0f0f0f0f) % 0xff;
 239   #elif 1
 240   r = x;
 241   r = ((r & 0xaaaaaaaa)>>1) + (r & 0x55555555);
 242   r = ((r & 0xcccccccc)>>2) + (r & 0x33333333);
 243   //r = ((r & 0xf0f0f0f0)>>4) + (r & 0x0f0f0f0f);
 244   r = ((r>>4) + r) & 0x0f0f0f0f;
 245   //r = ((r & 0xff00ff00)>>8) + (r & 0x00ff00ff);
 246   r = (r>>8) + r;
 247   //r = ((r & 0xffff0000)>>16) + (r & 0x0000ffff);
 248   r = ((r>>16) + r) & 63;
 249   #else
 250   r = _popCount[x & 0xff];
 251   x >>= 8;
 252   r += _popCount[x & 0xff];
 253   x >>= 8;
 254   r += _popCount[x & 0xff];
 255   x >>= 8;
 256   r += _popCount[x & 0xff];
 257   #endif
 258   return r;
 259 }
 260
 261
 262 unsigned int __popCount8(uint x) {
 263   uint r;
 264   #if 1
 265   r = x;
 266   r = ((r & 0xaa)>>1) + (r & 0x55);
 267   r = ((r & 0xcc)>>2) + (r & 0x33);
 268   r = ((r>>4) + r) & 0x0f;
 269   #else
 270   r = _popCount[x & 0xff];
 271   #endif
 272   return r;
 273 }
 274
 275
 276 int selectd2_save(selectd2 * s, FILE * fp) {
 277   uint wr = 0;
 278   wr += fwrite(&s->n,sizeof(uint),1,fp);
 279   wr += fwrite(&s->m,sizeof(uint),1,fp);
 280   wr += fwrite(&s->size,sizeof(uint),1,fp);
 281   wr += fwrite(&s->ss_len,sizeof(uint),1,fp);
 282   wr += fwrite(&s->sl_len,sizeof(uint),1,fp);
 283   wr += fwrite(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
 284   uint nl = (s->m-1) / L + 1;
 285   wr += fwrite(s->lp,sizeof(uint),nl+1,fp);
 286   wr += fwrite(s->p,sizeof(uint),nl+1,fp);
 287   wr += fwrite(s->ss,sizeof(ushort),s->ss_len,fp);
 288   wr += fwrite(s->sl,sizeof(uint),s->sl_len,fp);
 289   if(wr!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
 290     return 1;
 291   return 0;
 292 }
 293
 294
 295 int selectd2_load(selectd2 * s, FILE * fp) {
 296   uint rd = 0;
 297   rd += fread(&s->n,sizeof(uint),1,fp);
 298   rd += fread(&s->m,sizeof(uint),1,fp);
 299   rd += fread(&s->size,sizeof(uint),1,fp);
 300   rd += fread(&s->ss_len,sizeof(uint),1,fp);
 301   rd += fread(&s->sl_len,sizeof(uint),1,fp);
 302   s->buf = new uchar[(s->n+7)/8+1];
 303   rd += fread(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
 304   uint nl = (s->m-1) / L + 1;
 305   s->lp = new uint[nl+1];
 306   rd += fread(s->lp,sizeof(uint),nl+1,fp);
 307   s->p = new uint[nl+1];
 308   rd += fread(s->p,sizeof(uint),nl+1,fp);
 309   s->ss = new ushort[s->ss_len];
 310   rd += fread(s->ss,sizeof(ushort),s->ss_len,fp);
 311   s->sl = new uint[s->sl_len];
 312   rd += fread(s->sl,sizeof(uint),s->sl_len,fp);
 313   if(rd!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
 314     return 1;
 315   return 0;
 316 }
 317
 318
 319 void selectd2_free(selectd2 * s) {
 320   //delete [] s->buf;
 321   delete [] s->lp;
 322   delete [] s->p;
 323   delete [] s->ss;
 324   delete [] s->sl;
 325 }
 326
 327
 328 int selectd2_construct(selectd2 *select, int n, uchar *buf) {
 329   int i,m;
 330   int nl;
 331   int p,pp;
 332   int il,is,ml,ms;
 333   int r;
 334   uint *s;
 335
 336   make___selecttbl();
 337
 338   if (L/LLL == 0) {
 339     printf("ERROR: L=%d LLL=%d\n",L,LLL);
 340     exit(1);
 341   }
 342
 343   m = 0;
 344   for (i=0; i<n; i++) m += __getbit2(buf,i);
 345   select->n = n;
 346   select->m = m;
 347   //printf("n=%d m=%d\n",n,m);
 348
 349   select->buf = buf;
 350
 351   s = new uint[m];
 352   m = 0;
 353   for (i=0; i<n; i++) {
 354     if (__getbit2(buf,i)) {
 355       m++;
 356       s[m-1] = i;
 357     }
 358   }
 359
 360   nl = (m-1) / L + 1;
 361   select->size = 0;              //ignoring buf, shared with selects3
 362   select->lp = new uint[nl+1];
 363   for(int k=0;k<nl+1;k++) select->lp[k]=0;
 364   select->size += (nl+1)*sizeof(uint);
 365   select->p = new uint[nl+1];
 366   for(int k=0;k<nl+1;k++) select->p[k]=0;
 367   select->size += (nl+1)*sizeof(uint);
 368
 369   for (r = 0; r < 2; r++) {
 370     ml = ms = 0;
 371     for (il = 0; il < nl; il++) {
 372       pp = s[il*L];
 373       select->lp[il] = pp;
 374       i = min((il+1)*L-1,m-1);
 375       p = s[i];
 376       //printf("%d ",p-pp);
 377       if (p - pp >= LL) {
 378         if (r == 1) {
 379           for (is = 0; is < L; is++) {
 380             if (il*L+is >= m) break;
 381             select->sl[ml*L+is] = s[il*L+is];
 382           }
 383         }
 384         select->p[il] = -((ml<<logL)+1);
 385         ml++;
 386       }
 387       else {
 388         if (r == 1) {
 389           for (is = 0; is < L/LLL; is++) {
 390             if (il*L+is*LLL >= m) break;
 391             select->ss[ms*(L/LLL)+is] = s[il*L+is*LLL] - pp;
 392           }
 393         }
 394         select->p[il] = ms << (logL-logLLL);
 395         ms++;
 396       }
 397     }
 398     if (r == 0) {
 399       select->sl = new uint[ml*L+1];
 400       for(int k=0;k<ml*L+1;k++) select->sl[k]=0;
 401       select->size += sizeof(uint)*(ml*L+1);
 402       select->sl_len = ml*L+1;
 403       select->ss = new ushort[ms*(L/LLL)+1];
 404       for(int k=0;k<ms*(L/LLL)+1;k++) select->ss[k]=0;
 405       select->ss_len = ms*(L/LLL)+1;
 406       select->size += sizeof(ushort)*(ms*(L/LLL)+1);
 407     }
 408   }
 409   delete [] s;
 410   return 0;
 411 }
 412
 413
 414 int selectd2_select1(selectd2 *select, int i) {
 415   int p,r;
 416   int il;
 417   int rr;
 418   uchar *q;
 419   if (i <= 0) return -1;
 420   i--;
 421
 422   il = select->p[i>>logL];
 423   if (il < 0) {
 424     il = -il-1;
 425     //p = select->sl[(il<<logL)+(i & (L-1))];
 426     p = select->sl[il+(i & (L-1))];
 427   }
 428   else {
 429     p = select->lp[i>>logL];
 430     p += select->ss[il+((i & (L-1))>>logLLL)];
 431     r = i - (i & (LLL-1));
 432
 433     q = &(select->buf[p>>3]);
 434
 435     rr = p & (8-1);
 436     r -= _fast_popcount(*q >> (8-1-rr));
 437
 438     while (1) {
 439         //rr = _popCount[*q];
 440       rr = _fast_popcount(*q);
 441       if (r + rr >= i) break;
 442       r += rr;
 443       //p += 8;
 444       q++;
 445     }
 446       p = (q - select->buf) << 3;
 447       p += __selecttbl[((i-r-1)<<8)+(*q)];
 448   }
 449   return p;
 450 }
 451
 452 int selectd2_select0(selectd2 *select, int i) {
 453   int p,r;
 454   int il;
 455   int rr;
 456   uchar *q;
 457   if (i <= 0) return -1;
 458   i--;
 459
 460   il = select->p[i>>logL];
 461   if (il < 0) {
 462     il = -il-1;
 463     //p = select->sl[(il<<logL)+(i & (L-1))];
 464     p = select->sl[il+(i & (L-1))];
 465   }
 466   else {
 467     p = select->lp[i>>logL];
 468     //p += select->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
 469     p += select->ss[il+((i & (L-1))>>logLLL)];
 470     r = i - (i & (LLL-1));
 471
 472     q = &(select->buf[p>>3]);
 473
 474     rr = p & (8-1);
 475
 476     r -= _fast_popcount((*q ^ 0xff) >> (8-1-rr));
 477
 478     while (1) {
 479       //rr = _popCount[*q ^ 0xff];
 480       rr = _fast_popcount(*q ^ 0xff);
 481       if (r + rr >= i) break;
 482       r += rr;
 483       //p += 8;
 484       q++;
 485     }
 486     p = (q - select->buf) << 3;
 487     p += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
 488   }
 489   return p;
 490 }
 491
 492 int selectd2_select(selectd2 *select, int i,int f) {
 493   return f ? selectd2_select1(select, i) :
 494     selectd2_select0(select, i);
 495 }
 496
 497
 498 int selectd2_select2(selectd2 *select, int i,int f, int *st, int *en) {
 499   int p,r,p2;
 500   int il;
 501   int rr;
 502   uchar *q;
 503
 504   if (i == 0) {
 505     *st = -1;
 506     return -1;
 507   }
 508
 509   #if 0
 510   if (i > select->m) {
 511     printf("ERROR: m=%d i=%d\n",select->m,i);
 512     exit(1);
 513   }
 514   #endif
 515
 516   i--;
 517
 518   il = select->p[i>>logL];
 519   if (il < 0) {
 520     il = -il-1;
 521     //p = select->sl[(il<<logL)+(i & (L-1))];
 522     p = select->sl[il+(i & (L-1))];
 523
 524     if ((i>>logL) == ((i+1)>>logL)) {
 525       p2 = select->sl[il+((i+1) & (L-1))];
 526     }
 527     else {
 528       p2 = selectd2_select(select,i+2,f);
 529     }
 530   }
 531   else {
 532     p = select->lp[i>>logL];
 533     //p += select->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
 534     p += select->ss[il+((i & (L-1))>>logLLL)];
 535     r = i - (i & (LLL-1));
 536
 537     q = &(select->buf[p>>3]);
 538
 539     if (f == 1) {
 540       rr = p & (8-1);
 541       //r -= _popCount[*q >> (8-1-rr)];
 542       r -= _fast_popcount(*q >> (8-1-rr));
 543       //p = p - rr;
 544
 545       while (1) {
 546         //rr = _popCount[*q];
 547         rr = _fast_popcount(*q);
 548         if (r + rr >= i) break;
 549         r += rr;
 550         //p += 8;
 551         q++;
 552       }
 553       p = (q - select->buf) << 3;
 554       p += __selecttbl[((i-r-1)<<8)+(*q)];
 555
 556       if ((i>>logL) == ((i+1)>>logL)) {
 557         i++;
 558         while (1) {
 559           //rr = _popCount[*q];
 560           r = _fast_popcount(*q);
 561           if (r + rr >= i) break;
 562           r += rr;
 563           q++;
 564         }
 565         p2 = (q - select->buf) << 3;
 566         p2 += __selecttbl[((i-r-1)<<8)+(*q)];
 567       }
 568       else {
 569         p2 = selectd2_select(select,i+2,f);
 570       }
 571
 572     }
 573     else {
 574       rr = p & (8-1);
 575       //r -= _popCount[(*q ^ 0xff) >> (8-1-rr)];
 576       r -= _fast_popcount((*q ^ 0xff) >> (8-1-rr));
 577       //p = p - rr;
 578
 579       while (1) {
 580         //rr = _popCount[*q ^ 0xff];
 581         rr = _fast_popcount(*q ^ 0xff);
 582         if (r + rr >= i) break;
 583         r += rr;
 584         //p += 8;
 585         q++;
 586       }
 587       p = (q - select->buf) << 3;
 588       p += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
 589
 590       if ((i>>logL) == ((i+1)>>logL)) {
 591         i++;
 592         while (1) {
 593           //rr = _popCount[*q ^ 0xff];
 594           rr = _fast_popcount(*q ^ 0xff);
 595           if (r + rr >= i) break;
 596           r += rr;
 597           q++;
 598         }
 599         p2 = (q - select->buf) << 3;
 600         p2 += __selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
 601       }
 602       else {
 603         p2 = selectd2_select(select,i+2,f);
 604       }
 605     }
 606   }
 607   *st = p;
 608   *en = p2;
 609   return p;
 610 }
 611
 612
 613 int selects3_save(selects3 * s, FILE * fp) {
 614   uint wr = 0;
 615   wr += fwrite(&s->n,sizeof(uint),1,fp);
 616   wr += fwrite(&s->m,sizeof(uint),1,fp);
 617   wr += fwrite(&s->size,sizeof(uint),1,fp);
 618   wr += fwrite(&s->d,sizeof(uint),1,fp);
 619   wr += fwrite(&s->hi_len,sizeof(uint),1,fp);
 620   wr += fwrite(&s->low_len,sizeof(uint),1,fp);
 621   wr += fwrite(s->hi,sizeof(uchar),s->hi_len,fp);
 622   wr += fwrite(s->low,sizeof(uint),s->low_len,fp);
 623   if(wr!=(6+s->hi_len+s->low_len))
 624     return 1;
 625   if(selectd2_save(s->sd0,fp)) return 2;
 626   if(selectd2_save(s->sd1,fp)) return 3;
 627   return 0;
 628 }
 629
 630
 631 int selects3_load(selects3 * s, FILE * fp) {
 632   uint rd = 0;
 633   rd += fread(&s->n,sizeof(uint),1,fp);
 634   rd += fread(&s->m,sizeof(uint),1,fp);
 635   rd += fread(&s->size,sizeof(uint),1,fp);
 636   rd += fread(&s->d,sizeof(uint),1,fp);
 637   rd += fread(&s->hi_len,sizeof(uint),1,fp);
 638   rd += fread(&s->low_len,sizeof(uint),1,fp);
 639   s->hi = new uchar[s->hi_len];
 640   rd += fread(s->hi,sizeof(uchar),s->hi_len,fp);
 641   s->low = new uint[s->low_len];
 642   rd += fread(s->low,sizeof(uint),s->low_len,fp);
 643   if(rd!=(6+s->hi_len+s->low_len))
 644     return 1;
 645   s->sd0 = new selectd2;
 646   if(selectd2_load(s->sd0,fp)) return 2;
 647   s->sd1 = new selectd2;
 648   if(selectd2_load(s->sd1,fp)) return 3;
 649   delete [] s->sd0->buf;
 650   delete [] s->sd1->buf;
 651   s->sd0->buf = s->hi;
 652   s->sd1->buf = s->hi;
 653   return 0;
 654 }
 655
 656
 657 void selects3_free(selects3 * s) {
 658   delete [] s->hi;
 659   delete [] s->low;
 660   //delete [] s->sd0->buf;
 661   selectd2_free(s->sd0);
 662   delete s->sd0;
 663   selectd2_free(s->sd1);
 664   delete s->sd1;
 665 }
 666
 667
 668 int selects3_construct(selects3 *select, int n, uint *buf) {
 669   int i,m;
 670   int d,mm;
 671   uint *low;
 672   uchar *buf2;
 673   selectd2 *sd0,*sd1;
 674
 675   m = 0;
 676   for (i=0; i<n; i++) m += __getbit(buf,i);
 677   select->n = n;
 678   select->m = m;
 679
 680   if (m == 0) return 0;
 681
 682   mm = m;
 683   d = 0;
 684   while (mm < n) {
 685     mm <<= 1;
 686     d++;
 687   }
 688
 689   select->d = d;
 690
 691   buf2 = new uchar[(2*m+8-1)/8+1];
 692   for(int k=0;k<(2*m+8-1)/8+1;k++) buf2[k]=0;
 693   select->hi_len = (2*m+8-1)/8+1;
 694   low = new uint[(d*m+PBS-1)/PBS+1];
 695   for(uint k=0;k<(d*m+PBS-1)/PBS+1;k++) low[k]=0;
 696   select->low_len = (d*m+PBS-1)/PBS+1;
 697
 698   select->hi = buf2;
 699   select->low = low;
 700   select->size = sizeof(uchar)*((2*m+8-1)/8+1) + sizeof(uint)*((d*m+PBS-1)/PBS+1);
 701
 702   for (i=0; i<m*2; i++) __setbit2(buf2,i,0);
 703
 704   m = 0;
 705   for (i=0; i<n; i++) {
 706     if (__getbit(buf,i)) {
 707       __setbit2(buf2,(i>>d)+m,1);
 708       __setbits(low,m*d,d,i & ((1<<d)-1));
 709       m++;
 710     }
 711   }
 712
 713   sd1 = new selectd2;
 714   sd0 = new selectd2;
 715   select->size += 2*sizeof(selectd2);
 716
 717   selectd2_construct(sd1,m*2,buf2);
 718   select->sd1 = sd1;
 719
 720   for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
 721   selectd2_construct(sd0,m*2,buf2);
 722   select->sd0 = sd0;
 723
 724   for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
 725   return 0;
 726 }
 727
 728 //selects3 * lasts3=NULL;
 729 //int lasti=0;
 730 //int lasts=0;
 731
 732 int selects3_select(selects3 *select, int i) {
 733   int d,x;
 734
 735   #if 0
 736   if (i > select->m) {
 737     printf("ERROR: m=%d i=%d\n",select->m,i);
 738     exit(1);
 739   }
 740   #endif
 741
 742   if (i == 0) return -1;
 743
 744   d = select->d;
 745         /*if(select->lasti==(uint)i-1) {
 746                 while(!__getbit2(select->sd1->buf,++select->lasts));
 747         }
 748         else {
 749           select->lasts = selectd2_select(select->sd1,i,1);
 750         }
 751         select->lasti = i;
 752         //lasts3 = select; */
 753         x = selectd2_select1(select->sd1,i) - (i-1);
 754   //x = (select->lasts-(i-1)) << d;
 755   x <<= d;
 756   x += __getbits(select->low,(i-1)*d,d);
 757   return x;
 758 }
 759
 760 void pr_byte(FILE* fp, uchar b)
 761 {
 762   uchar * buff = &b;
 763   for(int i = 0; i < 8; i++){
 764     fprintf(stderr, "%i", __getbit2(buff, i));
 765   };
 766 }
 767
 768 int selects3_selectnext(selects3 *select, int i) {
 769         //return selects3_select(select,selects3_rank(select,i)+1);
 770   int d,x,w,y;
 771   int r,j;
 772   int z,ii;
 773   int xoffset;
 774   uint *q;
 775   d = select->d;
 776   q = select->low;
 777   ii = i>>d;
 778   y = selectd2_select0(select->sd0,ii)+1;
 779   int k2=y-ii;
 780   x = y - ii;
 781   int x_orig = x;
 782   j = i - (ii<<d);
 783   r = y & 7;
 784   y >>= 3;
 785   z = select->hi[y];
 786   xoffset = x * d;
 787   while (1) {
 788     if (((z << r) & 0x80) == 0) {
 789       k2 += (x!=x_orig);
 790       break;
 791     };
 792
 793     w = __getbits(q,xoffset,d);
 794     if (w >= j) {
 795       bool t1 = (w == j);
 796       bool t2 = (__getbit2(select->hi,((y << 3)+r)));
 797       if (t2)  k2 += (t1);
 798       x  += t1;
 799       r  += t1;
 800       break;
 801     };
 802
 803     x++;
 804     r++;
 805     xoffset += d;
 806     if(__getbit2(select->hi,( (y << 3)+r))) k2++;
 807     if (r == 8) {
 808       r = 0;
 809       y++;
 810       z = select->hi[y];
 811     }
 812   };
 813
 814   if(x==select->m) return (uint)-1;
 815
 816
 817   int c=(y << 3)+r;
 818   unsigned int mask = 0x00ffffff;
 819   unsigned int tmp = (select->hi[y] << (24+r)) | mask;
 820   unsigned int c_incr = __builtin_clz(tmp);
 821   c += (c_incr & 7);
 822   if (c_incr == 8) {
 823     c += (8-r);
 824     int pp = c >> 3;
 825     while(select->hi[pp]==0) {
 826       pp++;
 827       c += 8;
 828     };
 829     while(!__getbit2(select->hi,c)) c++;
 830
 831   };
 832   c -= (k2);
 833
 834   return __getbits(q,x*d,d)+((c)<<d);
 835 }
 836
 837 int selects3_rank(selects3 *select, int i) {
 838   int d,x,w,y;
 839   int r,j;
 840   int z,ii;
 841   int xoffset;
 842   uint *q;
 843
 844   d = select->d;
 845   q = select->low;
 846
 847   ii = i>>d;
 848
 849   y = selectd2_select0(select->sd0, ii)+1;
 850   //  selectd2_select2(select->sd0,ii,0,&y1,&y2);
 851   //y1++;  y2++;
 852   //printf("y %d y1 %d  %d\n",y,y1,y2-y1);
 853
 854   x = y - ii;
 855
 856   j = i - (ii<<d);
 857
 858   r = y & 7;
 859   y >>= 3;
 860   z = select->hi[y];
 861   xoffset = x * d;
 862   while (1) {
 863     if (((z << r) & 0x80) == 0) break;
 864     w = __getbits(q, xoffset, d);
 865     if (w >= j) {
 866       x += (w == j);
 867       break;
 868     }
 869     x++;
 870     r++;
 871     xoffset += d;
 872     if (r == 8) {
 873       r = 0;
 874       y++;
 875       z = select->hi[y];
 876     }
 877   }
 878
 879         return x;
 880 }
 881