libcds/src/static_bitsequence/sdarray.cpp

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