src/static_bitsequence/sdarray.cpp

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