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 static uint __getbits(uint *B, int i, int d)
 133 {
 134   ulong x;
 135   B += (i >> logD);
 136   i &= (D-1);
 137   x = ((ulong *) B)[0];
 138   x = (x << 32)|(x >> 32);
 139   x = (x << i) >> (2*D - d);
 140   return x;
 141 }
 142
 143 #endif
 144
 145 #if HAS_NATIVE_POPCOUNT
 146 static inline unsigned int popcount(unsigned int n){
 147   asm ("popcnt %1, %0" : "=r" (n) : "0" (n));
 148   return n;
 149 }
 150
 151 // static inline unsigned int popcount8(unsigned int n) {
 152 //   return __builtin_popcount(n & 0xff);
 153 // }
 154
 155 static inline unsigned int _fast_popcount(uchar x)
 156 {
 157   return __builtin_popcount(x);
 158 }
 159
 160
 161 #else
 162 static const unsigned int _popCount[] = {
 163   0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,
 164   1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
 165   1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
 166   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 167   1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
 168   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 169   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 170   3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
 171   1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
 172   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 173   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 174   3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
 175   2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
 176   3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
 177   3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
 178   4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8
 179 };
 180
 181 static inline unsigned int
 182 _fast_popcount(int x) {
 183   return _popCount[x & 0xff];
 184 }
 185 #endif
 186
 187
 188
 189 //static unsigned int __selecttbl[8*256];
 190 static int built = 0;
 191
 192 void make___selecttbl(void) {
 193   if(built) return;
 194   built = 1;
 195   int i,x,r;
 196   uint buf[1] = { 0 };
 197
 198   for (x = 0; x < 256; x++) {
 199     __setbits(buf,0,8,x);
 200 //    for (r=0; r<8; r++) __selecttbl[(r<<8)+x] = -1;
 201     r = 0;
 202     for (i=0; i<8; i++) {
 203       if (__getbit(buf,i)) {
 204         //      __selecttbl[(r<<8)+x] = i;
 205         r++;
 206       }
 207     }
 208   }
 209 }
 210
 211 int selectd2_save(selectd2 * s, FILE * fp) {
 212   uint wr = 0;
 213   wr += fwrite(&s->n,sizeof(uint),1,fp);
 214   wr += fwrite(&s->m,sizeof(uint),1,fp);
 215   wr += fwrite(&s->size,sizeof(uint),1,fp);
 216   wr += fwrite(&s->ss_len,sizeof(uint),1,fp);
 217   wr += fwrite(&s->sl_len,sizeof(uint),1,fp);
 218   wr += fwrite(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
 219   uint nl = (s->m-1) / L + 1;
 220   wr += fwrite(s->lp,sizeof(uint),nl+1,fp);
 221   wr += fwrite(s->p,sizeof(uint),nl+1,fp);
 222   wr += fwrite(s->ss,sizeof(ushort),s->ss_len,fp);
 223   wr += fwrite(s->sl,sizeof(uint),s->sl_len,fp);
 224   if(wr!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
 225     return 1;
 226   return 0;
 227 }
 228
 229
 230 int selectd2_load(selectd2 * s, FILE * fp) {
 231   uint rd = 0;
 232   rd += fread(&s->n,sizeof(uint),1,fp);
 233   rd += fread(&s->m,sizeof(uint),1,fp);
 234   rd += fread(&s->size,sizeof(uint),1,fp);
 235   rd += fread(&s->ss_len,sizeof(uint),1,fp);
 236   rd += fread(&s->sl_len,sizeof(uint),1,fp);
 237   s->buf = new uchar[(s->n+7)/8+1];
 238   rd += fread(s->buf,sizeof(uchar),(s->n+7)/8+1,fp);
 239   uint nl = (s->m-1) / L + 1;
 240   s->lp = new uint[nl+1];
 241   rd += fread(s->lp,sizeof(uint),nl+1,fp);
 242   s->p = new uint[nl+1];
 243   rd += fread(s->p,sizeof(uint),nl+1,fp);
 244   s->ss = new ushort[s->ss_len];
 245   rd += fread(s->ss,sizeof(ushort),s->ss_len,fp);
 246   s->sl = new uint[s->sl_len];
 247   rd += fread(s->sl,sizeof(uint),s->sl_len,fp);
 248   if(rd!=s->sl_len+s->ss_len+2*(nl+1)+(s->n+7)/8+1+5)
 249     return 1;
 250   return 0;
 251 }
 252
 253
 254 void selectd2_free(selectd2 * s) {
 255   //delete [] s->buf;
 256   delete [] s->lp;
 257   delete [] s->p;
 258   delete [] s->ss;
 259   delete [] s->sl;
 260 }
 261
 262
 263 int selectd2_construct(selectd2 *select, int n, uchar *buf) {
 264   int i,m;
 265   int nl;
 266   int p,pp;
 267   int il,is,ml,ms;
 268   int r;
 269   uint *s;
 270
 271   //make___selecttbl();
 272
 273   if (L/LLL == 0) {
 274     printf("ERROR: L=%d LLL=%d\n",L,LLL);
 275     exit(1);
 276   }
 277
 278   m = 0;
 279   for (i=0; i<n; i++) m += __getbit2(buf,i);
 280   select->n = n;
 281   select->m = m;
 282   //printf("n=%d m=%d\n",n,m);
 283
 284   select->buf = buf;
 285
 286   s = new uint[m];
 287   m = 0;
 288   for (i=0; i<n; i++) {
 289     if (__getbit2(buf,i)) {
 290       m++;
 291       s[m-1] = i;
 292     }
 293   }
 294
 295   nl = (m-1) / L + 1;
 296   select->size = 0;              //ignoring buf, shared with selects3
 297   select->lp = new uint[nl+1];
 298   for(int k=0;k<nl+1;k++) select->lp[k]=0;
 299   select->size += (nl+1)*sizeof(uint);
 300   select->p = new uint[nl+1];
 301   for(int k=0;k<nl+1;k++) select->p[k]=0;
 302   select->size += (nl+1)*sizeof(uint);
 303
 304   for (r = 0; r < 2; r++) {
 305     ml = ms = 0;
 306     for (il = 0; il < nl; il++) {
 307       pp = s[il*L];
 308       select->lp[il] = pp;
 309       i = min((il+1)*L-1,m-1);
 310       p = s[i];
 311       //printf("%d ",p-pp);
 312       if (p - pp >= LL) {
 313         if (r == 1) {
 314           for (is = 0; is < L; is++) {
 315             if (il*L+is >= m) break;
 316             select->sl[ml*L+is] = s[il*L+is];
 317           }
 318         }
 319         select->p[il] = -((ml<<logL)+1);
 320         ml++;
 321       }
 322       else {
 323         if (r == 1) {
 324           for (is = 0; is < L/LLL; is++) {
 325             if (il*L+is*LLL >= m) break;
 326             select->ss[ms*(L/LLL)+is] = s[il*L+is*LLL] - pp;
 327           }
 328         }
 329         select->p[il] = ms << (logL-logLLL);
 330         ms++;
 331       }
 332     }
 333     if (r == 0) {
 334       select->sl = new uint[ml*L+1];
 335       for(int k=0;k<ml*L+1;k++) select->sl[k]=0;
 336       select->size += sizeof(uint)*(ml*L+1);
 337       select->sl_len = ml*L+1;
 338       select->ss = new ushort[ms*(L/LLL)+1];
 339       for(int k=0;k<ms*(L/LLL)+1;k++) select->ss[k]=0;
 340       select->ss_len = ms*(L/LLL)+1;
 341       select->size += sizeof(ushort)*(ms*(L/LLL)+1);
 342     }
 343   }
 344   delete [] s;
 345   return 0;
 346 }
 347
 348
 349 int selectd2_select1(selectd2 *select, int i) {
 350   int p,r;
 351   int il;
 352   int rr;
 353   uchar *q;
 354   if (i <= 0) return -1;
 355   i--;
 356
 357   il = select->p[i>>logL];
 358   if (il < 0) {
 359     il = -il-1;
 360     //p = select->sl[(il<<logL)+(i & (L-1))];
 361     p = select->sl[il+(i & (L-1))];
 362   }
 363   else {
 364     p = select->lp[i>>logL];
 365     p += select->ss[il+((i & (L-1))>>logLLL)];
 366     r = i - (i & (LLL-1));
 367
 368     q = &(select->buf[p>>3]);
 369
 370     rr = p & (8-1);
 371     r -= _fast_popcount(*q >> (8-1-rr));
 372
 373     while (1) {
 374         //rr = _popCount[*q];
 375       rr = _fast_popcount(*q);
 376       if (r + rr >= i) break;
 377       r += rr;
 378       //p += 8;
 379       q++;
 380     }
 381       p = (q - select->buf) << 3;
 382       p += selecttbl[((i-r-1)<<8)+(*q)];
 383   }
 384   return p;
 385 }
 386
 387 int selectd2_select0(selectd2 *select, int i) {
 388   int p,r;
 389   int il;
 390   int rr;
 391   uchar *q;
 392
 393   if (i <= 0) return -1;
 394   i--;
 395
 396   il = select->p[i>>logL];
 397   if (il < 0) {
 398     il = -il-1;
 399     return select->sl[il+(i & (L-1))];
 400
 401   } else {
 402     p = select->lp[i>>logL];
 403
 404     p += select->ss[il+((i & (L-1))>>logLLL)];
 405     r = i - (i & (LLL-1));
 406
 407     q = &(select->buf[p>>3]);
 408
 409     rr = p & 7;
 410
 411     uint masked_q = *q ^ 0xff;
 412
 413     r -= _fast_popcount(masked_q >> (7-rr));
 414
 415     for(;;) {
 416       rr = _fast_popcount(masked_q);
 417       if (rr >= i-r) {
 418         p = (q - select->buf) << 3;
 419         p += selecttbl[((i-r-1)<<8)+masked_q];
 420         return p;
 421       };
 422       r += rr;
 423       q++;
 424       masked_q = *q ^ 0xff;
 425     };
 426   }
 427 }
 428
 429 int selectd2_select(selectd2 *select, int i,int f) {
 430   return f ? selectd2_select1(select, i) :
 431     selectd2_select0(select, i);
 432 }
 433
 434
 435 int selectd2_select2(selectd2 *select, int i,int f, int *st, int *en) {
 436   int p,r,p2;
 437   int il;
 438   int rr;
 439   uchar *q;
 440
 441   if (i == 0) {
 442     *st = -1;
 443     return -1;
 444   }
 445
 446   #if 0
 447   if (i > select->m) {
 448     printf("ERROR: m=%d i=%d\n",select->m,i);
 449     exit(1);
 450   }
 451   #endif
 452
 453   i--;
 454
 455   il = select->p[i>>logL];
 456   if (il < 0) {
 457     il = -il-1;
 458     //p = select->sl[(il<<logL)+(i & (L-1))];
 459     p = select->sl[il+(i & (L-1))];
 460
 461     if ((i>>logL) == ((i+1)>>logL)) {
 462       p2 = select->sl[il+((i+1) & (L-1))];
 463     }
 464     else {
 465       p2 = selectd2_select(select,i+2,f);
 466     }
 467   }
 468   else {
 469     p = select->lp[i>>logL];
 470     //p += select->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];
 471     p += select->ss[il+((i & (L-1))>>logLLL)];
 472     r = i - (i & (LLL-1));
 473
 474     q = &(select->buf[p>>3]);
 475
 476     if (f == 1) {
 477       rr = p & (8-1);
 478       //r -= _popCount[*q >> (8-1-rr)];
 479       r -= _fast_popcount(*q >> (8-1-rr));
 480       //p = p - rr;
 481
 482       while (1) {
 483         //rr = _popCount[*q];
 484         rr = _fast_popcount(*q);
 485         if (r + rr >= i) break;
 486         r += rr;
 487         //p += 8;
 488         q++;
 489       }
 490       p = (q - select->buf) << 3;
 491       p += selecttbl[((i-r-1)<<8)+(*q)];
 492
 493       if ((i>>logL) == ((i+1)>>logL)) {
 494         i++;
 495         while (1) {
 496           //rr = _popCount[*q];
 497           r = _fast_popcount(*q);
 498           if (r + rr >= i) break;
 499           r += rr;
 500           q++;
 501         }
 502         p2 = (q - select->buf) << 3;
 503         p2 += selecttbl[((i-r-1)<<8)+(*q)];
 504       }
 505       else {
 506         p2 = selectd2_select(select,i+2,f);
 507       }
 508
 509     }
 510     else {
 511       rr = p & (8-1);
 512       //r -= _popCount[(*q ^ 0xff) >> (8-1-rr)];
 513       r -= _fast_popcount((*q ^ 0xff) >> (8-1-rr));
 514       //p = p - rr;
 515
 516       while (1) {
 517         //rr = _popCount[*q ^ 0xff];
 518         rr = _fast_popcount(*q ^ 0xff);
 519         if (r + rr >= i) break;
 520         r += rr;
 521         //p += 8;
 522         q++;
 523       }
 524       p = (q - select->buf) << 3;
 525       p += selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
 526
 527       if ((i>>logL) == ((i+1)>>logL)) {
 528         i++;
 529         while (1) {
 530           //rr = _popCount[*q ^ 0xff];
 531           rr = _fast_popcount(*q ^ 0xff);
 532           if (r + rr >= i) break;
 533           r += rr;
 534           q++;
 535         }
 536         p2 = (q - select->buf) << 3;
 537         p2 += selecttbl[((i-r-1)<<8)+(*q ^ 0xff)];
 538       }
 539       else {
 540         p2 = selectd2_select(select,i+2,f);
 541       }
 542     }
 543   }
 544   *st = p;
 545   *en = p2;
 546   return p;
 547 }
 548
 549
 550 int selects3_save(selects3 * s, FILE * fp) {
 551   uint wr = 0;
 552   wr += fwrite(&s->n,sizeof(uint),1,fp);
 553   wr += fwrite(&s->m,sizeof(uint),1,fp);
 554   wr += fwrite(&s->size,sizeof(uint),1,fp);
 555   wr += fwrite(&s->d,sizeof(uint),1,fp);
 556   wr += fwrite(&s->hi_len,sizeof(uint),1,fp);
 557   wr += fwrite(&s->low_len,sizeof(uint),1,fp);
 558   wr += fwrite(s->hi,sizeof(uchar),s->hi_len,fp);
 559   wr += fwrite(s->low,sizeof(uint),s->low_len,fp);
 560   if(wr!=(6+s->hi_len+s->low_len))
 561     return 1;
 562   if(selectd2_save(s->sd0,fp)) return 2;
 563   if(selectd2_save(s->sd1,fp)) return 3;
 564   return 0;
 565 }
 566
 567
 568 int selects3_load(selects3 * s, FILE * fp) {
 569   uint rd = 0;
 570   rd += fread(&s->n,sizeof(uint),1,fp);
 571   rd += fread(&s->m,sizeof(uint),1,fp);
 572   rd += fread(&s->size,sizeof(uint),1,fp);
 573   rd += fread(&s->d,sizeof(uint),1,fp);
 574   rd += fread(&s->hi_len,sizeof(uint),1,fp);
 575   rd += fread(&s->low_len,sizeof(uint),1,fp);
 576   s->hi = new uchar[s->hi_len];
 577   rd += fread(s->hi,sizeof(uchar),s->hi_len,fp);
 578   s->low = new uint[s->low_len];
 579   rd += fread(s->low,sizeof(uint),s->low_len,fp);
 580   if(rd!=(6+s->hi_len+s->low_len))
 581     return 1;
 582   s->sd0 = new selectd2;
 583   if(selectd2_load(s->sd0,fp)) return 2;
 584   s->sd1 = new selectd2;
 585   if(selectd2_load(s->sd1,fp)) return 3;
 586   delete [] s->sd0->buf;
 587   delete [] s->sd1->buf;
 588   s->sd0->buf = s->hi;
 589   s->sd1->buf = s->hi;
 590   return 0;
 591 }
 592
 593
 594 void selects3_free(selects3 * s) {
 595   delete [] s->hi;
 596   delete [] s->low;
 597   //delete [] s->sd0->buf;
 598   selectd2_free(s->sd0);
 599   delete s->sd0;
 600   selectd2_free(s->sd1);
 601   delete s->sd1;
 602 }
 603
 604
 605 int selects3_construct(selects3 *select, int n, uint *buf) {
 606   int i,m;
 607   int d,mm;
 608   uint *low;
 609   uchar *buf2;
 610   selectd2 *sd0,*sd1;
 611
 612   m = 0;
 613   for (i=0; i<n; i++) m += __getbit(buf,i);
 614   select->n = n;
 615   select->m = m;
 616
 617   if (m == 0) return 0;
 618
 619   mm = m;
 620   d = 0;
 621   while (mm < n) {
 622     mm <<= 1;
 623     d++;
 624   }
 625
 626   select->d = d;
 627
 628   buf2 = new uchar[(2*m+8-1)/8+1];
 629   for(int k=0;k<(2*m+8-1)/8+1;k++) buf2[k]=0;
 630   select->hi_len = (2*m+8-1)/8+1;
 631   low = new uint[(d*m+PBS-1)/PBS+1];
 632   for(uint k=0;k<(d*m+PBS-1)/PBS+1;k++) low[k]=0;
 633   select->low_len = (d*m+PBS-1)/PBS+1;
 634
 635   select->hi = buf2;
 636   select->low = low;
 637   select->size = sizeof(uchar)*((2*m+8-1)/8+1) + sizeof(uint)*((d*m+PBS-1)/PBS+1);
 638
 639   for (i=0; i<m*2; i++) __setbit2(buf2,i,0);
 640
 641   m = 0;
 642   for (i=0; i<n; i++) {
 643     if (__getbit(buf,i)) {
 644       __setbit2(buf2,(i>>d)+m,1);
 645       __setbits(low,m*d,d,i & ((1<<d)-1));
 646       m++;
 647     }
 648   }
 649
 650   sd1 = new selectd2;
 651   sd0 = new selectd2;
 652   select->size += 2*sizeof(selectd2);
 653
 654   selectd2_construct(sd1,m*2,buf2);
 655   select->sd1 = sd1;
 656
 657   for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
 658   selectd2_construct(sd0,m*2,buf2);
 659   select->sd0 = sd0;
 660
 661   for (i=0; i<m*2; i++) __setbit2(buf2,i,1-__getbit2(buf2,i));
 662   return 0;
 663 }
 664
 665 //selects3 * lasts3=NULL;
 666 //int lasti=0;
 667 //int lasts=0;
 668
 669 int selects3_select(selects3 *select, int i) {
 670   int d,x;
 671
 672   #if 0
 673   if (i > select->m) {
 674     printf("ERROR: m=%d i=%d\n",select->m,i);
 675     exit(1);
 676   }
 677   #endif
 678
 679   if (i == 0) return -1;
 680
 681   d = select->d;
 682         /*if(select->lasti==(uint)i-1) {
 683                 while(!__getbit2(select->sd1->buf,++select->lasts));
 684         }
 685         else {
 686           select->lasts = selectd2_select(select->sd1,i,1);
 687         }
 688         select->lasti = i;
 689         //lasts3 = select; */
 690         x = selectd2_select1(select->sd1,i) - (i-1);
 691   //x = (select->lasts-(i-1)) << d;
 692   x <<= d;
 693   x += __getbits(select->low,(i-1)*d,d);
 694   return x;
 695 }
 696
 697 void pr_byte(FILE* fp, uchar b)
 698 {
 699   uchar * buff = &b;
 700   for(int i = 0; i < 8; i++){
 701     fprintf(stderr, "%i", __getbit2(buff, i));
 702   };
 703 }
 704
 705 int selects3_selectnext(selects3 *select, int i) {
 706   //return selects3_select(select,selects3_rank(select,i)+1);
 707   int d,x,w,y;
 708   int r,j;
 709   int z,ii;
 710   int xoffset;
 711   uint *q;
 712   d = select->d;
 713   q = select->low;
 714   ii = i>>d;
 715   y = selectd2_select0(select->sd0,ii)+1;
 716   int k2=y-ii;
 717   x = y - ii;
 718   int x_orig = x;
 719   j = i - (ii<<d);
 720   r = y & 7;
 721   y >>= 3;
 722   z = select->hi[y];
 723   xoffset = x * d;
 724   while (1) {
 725     if (((z << r) & 0x80) == 0) {
 726       k2 += (x!=x_orig);
 727       break;
 728     };
 729
 730     w = __getbits(q,xoffset,d);
 731     if (w >= j) {
 732       if (w == j) {
 733         if (__getbit2(select->hi,((y << 3)+r)))  k2 ++;
 734         x++;
 735         r++;
 736         xoffset += d;
 737       };
 738       break;
 739     };
 740
 741     x++;
 742     r++;
 743     xoffset += d;
 744     if(__getbit2(select->hi,( (y << 3)+r))) k2++;
 745     if (r == 8) {
 746       r = 0;
 747       y++;
 748       z = select->hi[y];
 749     }
 750   };
 751
 752   if(x==select->m) return (uint)-1;
 753
 754
 755   int c = (y << 3)+r;
 756
 757   uint mask = ~(0xffu << (8*sizeof(uint) - 8));
 758   uint tmp = (select->hi[y] << ((8*sizeof(uint) - 8)+r)) | mask;
 759   uint c_incr = __builtin_clz(tmp);
 760   c += (c_incr & 7);
 761   if (c_incr == 8) {
 762     c += (8-r);
 763
 764     uchar * pp = &(select->hi[c >> 3]);
 765
 766     while(*pp == 0) {
 767       pp++;
 768       c += 8;
 769     };
 770
 771     c += __builtin_clz(*pp) - 24;
 772
 773   };
 774   c -= k2;
 775
 776   return __getbits(q,xoffset,d)+((c)<<d);
 777 }
 778
 779 int selects3_rank(selects3 *select, int i) {
 780   int d,x,w,y;
 781   int r,j;
 782   int z,ii;
 783   int xoffset;
 784   uint *q;
 785
 786   d = select->d;
 787   q = select->low;
 788
 789   ii = i >> d;
 790
 791   y = selectd2_select0(select->sd0, ii)+1;
 792
 793   x = y - ii;
 794
 795   j = i - (ii << d);
 796
 797   r = y & 7;
 798   y >>= 3;
 799   z = select->hi[y];
 800   xoffset = x * d;
 801
 802   while (1) {
 803     if (((z << r) & 0x80) == 0) return x;
 804
 805     w = __getbits(q, xoffset, d);
 806
 807     if (w >= j) return x + (w == j);
 808
 809     x++;
 810     r++;
 811     xoffset += d;
 812     if (r == 8) {
 813       r = 0;
 814       y++;
 815       z = select->hi[y];
 816     }
 817   }
 818 }
 819