darray.c

   1 #include <stdio.h>\r
   2 #include <stdlib.h>\r
   3 #include "darray.h"\r
   4 \r
   5 //typedef unsigned char byte;\r
   6 //typedef unsigned short word;\r
   7 //typedef unsigned int dword;\r
   8 \r
   9 //typedef dword pb;\r
  10 #define logD 5\r
  11 \r
  12 #define PBS (sizeof(pb)*8)\r
  13 #define D (1<<logD)\r
  14 #define logM 5\r
  15 #define M (1<<logM)\r
  16 #define logP 8\r
  17 #define P (1<<logP)\r
  18 #define logLL 16    // size of word\r
  19 #define LL (1<<logLL)\r
  20 #define logLLL 7\r
  21 //#define logLLL 5\r
  22 #define LLL (1<<logLLL)\r
  23 //#define logL 10\r
  24 //#define logL (logLL-3)\r
  25 #define logL (logLL-1-5)\r
  26 #define L (1<<logL)\r
  27 \r
  28 #ifndef min\r
  29  #define min(x,y) ((x)<(y)?(x):(y))\r
  30 #endif\r
  31 \r
  32 #define mymalloc(p,n,f) {p = (__typeof__(p)) malloc((n)*sizeof(*p)); if ((p)==NULL) {printf("not enough memory\n"); exit(1);}}\r
  33 \r
  34 int setbit(pb *B, int i,int x)\r
  35 {\r
  36   int j,l;\r
  37 \r
  38   j = i / D;\r
  39   l = i % D;\r
  40   if (x==0) B[j] &= (~(1<<(D-1-l)));\r
  41   else if (x==1) B[j] |= (1<<(D-1-l));\r
  42   else {\r
  43     printf("error setbit x=%d\n",x);\r
  44     exit(1);\r
  45   }\r
  46   return x;\r
  47 }\r
  48 \r
  49 int setbits(pb *B, int i, int d, int x)\r
  50 {\r
  51   int j;\r
  52 \r
  53   for (j=0; j<d; j++) {\r
  54     setbit(B,i+j,(x>>(d-j-1))&1);\r
  55   }\r
  56   return x;\r
  57 }\r
  58 \r
  59 int getbit(pb *B, int i)\r
  60 {\r
  61   int j,l;\r
  62 \r
  63   //j = i / D;\r
  64   //l = i % D;\r
  65   j = i >> logD;\r
  66   l = i & (D-1);\r
  67   return (B[j] >> (D-1-l)) & 1;\r
  68 }\r
  69 \r
  70 dword getbits(pb *B, int i, int d)\r
  71 {\r
  72   dword j,x;\r
  73 \r
  74   x = 0;\r
  75   for (j=0; j<d; j++) {\r
  76     x <<= 1;\r
  77     x += getbit(B,i+j);\r
  78   }\r
  79   return x;\r
  80 }\r
  81 \r
  82 int getpattern(pb *B, int i, int k, pb pat)\r
  83 {\r
  84   int j;\r
  85   int x;\r
  86   x = 1;\r
  87   for (j=0; j<k; j++) {\r
  88     x &= getbit(B,i+j) ^ (~(pat>>(k-1-j)));\r
  89   }\r
  90   //printf("getpattern(%d) = %d\n",i,x);\r
  91   return x;\r
  92 }\r
  93 \r
  94 \r
  95 static const unsigned int popCount[] = {\r
  96 0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,\r
  97 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,\r
  98 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,\r
  99 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
 100 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,\r
 101 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
 102 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
 103 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,\r
 104 1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,\r
 105 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
 106 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
 107 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,\r
 108 2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,\r
 109 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,\r
 110 3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,\r
 111 4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8\r
 112 };\r
 113 \r
 114 static int selecttbl[8*256];\r
 115 \r
 116 unsigned int popcount(pb x)\r
 117 {\r
 118   pb r;\r
 119 #if 0\r
 120   r = x;\r
 121   r = r - ((r>>1) & 0x77777777) - ((r>>2) & 0x33333333) - ((r>>3) & 0x11111111);\r
 122   r = ((r + (r>>4)) & 0x0f0f0f0f) % 0xff;\r
 123 #elif 1\r
 124   r = x;\r
 125   r = ((r & 0xaaaaaaaa)>>1) + (r & 0x55555555);\r
 126   r = ((r & 0xcccccccc)>>2) + (r & 0x33333333);\r
 127   //r = ((r & 0xf0f0f0f0)>>4) + (r & 0x0f0f0f0f);\r
 128   r = ((r>>4) + r) & 0x0f0f0f0f;\r
 129   //r = ((r & 0xff00ff00)>>8) + (r & 0x00ff00ff);\r
 130   r = (r>>8) + r;\r
 131   //r = ((r & 0xffff0000)>>16) + (r & 0x0000ffff);\r
 132   r = ((r>>16) + r) & 63;\r
 133 #else\r
 134   r = popCount[x & 0xff];\r
 135   x >>= 8;\r
 136   r += popCount[x & 0xff];\r
 137   x >>= 8;\r
 138   r += popCount[x & 0xff];\r
 139   x >>= 8;\r
 140   r += popCount[x & 0xff];\r
 141 #endif\r
 142   return r;\r
 143 }\r
 144 \r
 145 unsigned int popcount8(pb x)\r
 146 {\r
 147   dword r;\r
 148 #if 1\r
 149   r = x;\r
 150   r = ((r & 0xaa)>>1) + (r & 0x55);\r
 151   r = ((r & 0xcc)>>2) + (r & 0x33);\r
 152   r = ((r>>4) + r) & 0x0f;\r
 153 #else\r
 154   r = popCount[x & 0xff];\r
 155 #endif\r
 156   return r;\r
 157 }\r
 158 \r
 159 void make_selecttbl(void)\r
 160 {\r
 161   int i,x,r;\r
 162   pb buf[1];\r
 163 \r
 164   for (x = 0; x < 256; x++) {\r
 165     setbits(buf,0,8,x);\r
 166     for (r=0; r<8; r++) selecttbl[(r<<8)+x] = -1;\r
 167     r = 0;\r
 168     for (i=0; i<8; i++) {\r
 169       if (getbit(buf,i)) {\r
 170         selecttbl[(r<<8)+x] = i;\r
 171         r++;\r
 172       }\r
 173     }\r
 174   }\r
 175 }\r
 176 \r
 177 \r
 178 int darray_construct(darray *da, int n, pb *buf, int opt)\r
 179 {\r
 180   int i,j,k,m;\r
 181   int nl;\r
 182   int p,pp;\r
 183   int il,is,ml,ms;\r
 184   int r,m2;\r
 185   dword *s;\r
 186   int p1,p2,p3,p4,s1,s2,s3,s4;\r
 187 \r
 188   da->idx_size = 0;\r
 189 \r
 190   make_selecttbl();\r
 191 \r
 192 \r
 193   if (L/LLL == 0) {\r
 194     printf("ERROR: L=%d LLL=%d\n",L,LLL);\r
 195     exit(1);\r
 196   }\r
 197 \r
 198   m = 0;\r
 199   for (i=0; i<n; i++) m += getbit(buf,i);\r
 200   da->n = n;\r
 201   da->m = m;\r
 202   //printf("n=%d m=%d\n",n,m);\r
 203 \r
 204   da->buf = buf;\r
 205 \r
 206   if (opt & (~OPT_NO_RANK)) {  // construct select table\r
 207 #if 0\r
 208     mymalloc(s,m,0);\r
 209     m = 0;\r
 210     for (i=0; i<n; i++) {\r
 211       if (getbit(buf,i)) {\r
 212         m++;\r
 213         s[m-1] = i;\r
 214       }\r
 215     }\r
 216 #endif    \r
 217     nl = (m-1) / L + 1;\r
 218     mymalloc(da->lp,nl+1,1);  da->idx_size += (nl+1)*sizeof(*da->lp);\r
 219     mymalloc(da->p,nl+1,1);  da->idx_size += (nl+1)*sizeof(*da->p);\r
 220 #if 0\r
 221     printf("lp table: %d bytes (%1.2f bpc)\n",(nl+1)*sizeof(*da->lp), (double)(nl+1)*sizeof(*da->lp) * 8/n);\r
 222     printf("p table: %d bytes (%1.2f bpc)\n",(nl+1)*sizeof(*da->p), (double)(nl+1)*sizeof(*da->p) * 8/n);\r
 223 #endif    \r
 224 \r
 225     for (r = 0; r < 2; r++) {\r
 226       s1 = s2 = s3 = s4 = 0;\r
 227       p1 = p2 = p3 = p4 = -1;\r
 228     \r
 229       ml = ms = 0;\r
 230       for (il = 0; il < nl; il++) {\r
 231         //pp = s[il*L];\r
 232         while (s1 <= il*L) {\r
 233           if (getbit(buf,p1+1)) s1++;\r
 234           p1++;\r
 235         }\r
 236         pp = p1;\r
 237         da->lp[il] = pp;\r
 238         i = min((il+1)*L-1,m-1);\r
 239         //p = s[i];\r
 240         while (s2 <= i) {\r
 241           if (getbit(buf,p2+1)) s2++;\r
 242           p2++;\r
 243         }\r
 244         p = p2;\r
 245         //printf("%d ",p-pp);\r
 246         if (p - pp >= LL) {\r
 247           if (r == 1) {\r
 248             for (is = 0; is < L; is++) {\r
 249               if (il*L+is >= m) break;\r
 250               //da->sl[ml*L+is] = s[il*L+is];\r
 251               while (s3 <= il*L+is) {\r
 252                 if (getbit(buf,p3+1)) s3++;\r
 253                 p3++;\r
 254               }\r
 255               da->sl[ml*L+is] = p3;\r
 256             }\r
 257           }\r
 258           da->p[il] = -(ml+1);\r
 259           ml++;\r
 260         } else {\r
 261           if (r == 1) {\r
 262             for (is = 0; is < L/LLL; is++) {\r
 263               if (il*L+is*LLL >= m) break;\r
 264               while (s4 <= il*L+is*LLL) {\r
 265                 if (getbit(buf,p4+1)) s4++;\r
 266                 p4++;\r
 267               }\r
 268               //da->ss[ms*(L/LLL)+is] = s[il*L+is*LLL] - pp;\r
 269               da->ss[ms*(L/LLL)+is] = p4 - pp;\r
 270             }\r
 271           }\r
 272           da->p[il] = ms;\r
 273           ms++;\r
 274         }\r
 275       }\r
 276       if (r == 0) {\r
 277         mymalloc(da->sl,ml*L+1,1);  da->idx_size += (ml*L+1)*sizeof(*da->sl);\r
 278         mymalloc(da->ss,ms*(L/LLL)+1,1);  da->idx_size += (ms*(L/LLL)+1)*sizeof(*da->ss);\r
 279 #if 0\r
 280         printf("sl table: %d bytes (%1.2f bpc)\n",(ml*L+1)*sizeof(*da->sl), (double)(ml*L+1)*sizeof(*da->sl) * 8/n);\r
 281         printf("ss table: %d bytes (%1.2f bpc)\n",(ms*(L/LLL)+1)*sizeof(*da->ss), (double)(ms*(L/LLL)+1)*sizeof(*da->ss) * 8/n);\r
 282 #endif\r
 283       }\r
 284     }\r
 285     //free(s);\r
 286   } else { // no select table\r
 287     da->lp = NULL;\r
 288     da->p = da->sl = NULL;\r
 289     da->ss = NULL;\r
 290   }\r
 291 \r
 292   // construct rank table\r
 293 \r
 294   if ((opt & OPT_NO_RANK) == 0) {\r
 295     mymalloc(da->rl,n/R1+2,1);  da->idx_size += (n/R1+2)*sizeof(*da->rl);\r
 296     mymalloc(da->rm,n/RR+2,1);  da->idx_size += (n/RR+2)*sizeof(*da->rm);\r
 297     mymalloc(da->rs,n/RRR+2,1);  da->idx_size += (n/RRR+2)*sizeof(*da->rs);\r
 298 #if 0\r
 299     printf("rl table: %d bytes (%1.2f bpc)\n",(n/R1+2)*sizeof(*da->rl), (double)(n/R1+2)*sizeof(*da->rl) * 8/n);\r
 300     printf("rm table: %d bytes (%1.2f bpc)\n",(n/RR+2)*sizeof(*da->rm), (double)(n/RR+2)*sizeof(*da->rm) * 8/n);\r
 301     printf("rs table: %d bytes (%1.2f bpc)\n",(n/RRR+2)*sizeof(*da->rs), (double)(n/RRR+2)*sizeof(*da->rs) * 8/n);\r
 302 #endif\r
 303     r = 0;\r
 304     for (k=0; k<=n; k+=R1) {\r
 305       da->rl[k/R1] = r;\r
 306       m2 = 0;\r
 307       for (i=0; i<R1; i+=RR) {\r
 308         if (k+i <= n) da->rm[(k+i)/RR] = m2;\r
 309         m = 0;\r
 310         for (j=0; j<RR; j++) {\r
 311           if (k+i+j < n && getbit(buf,k+i+j)==1) m++;\r
 312           if (j % RRR == RRR-1) {\r
 313             if (k+i+j <= n) da->rs[(k+i+j)/RRR] = m;\r
 314             m2 += m;\r
 315             m = 0;\r
 316           }\r
 317         }\r
 318         if (m != 0) {\r
 319           printf("???\n");\r
 320         }\r
 321         //m2 += m;\r
 322       }\r
 323       r += m2;\r
 324     }\r
 325   }\r
 326 \r
 327   return 0;\r
 328 }\r
 329 \r
 330 // destroyDarray: frees the memory of darray\r
 331 // Added by Diego Arroyuelo\r
 332 void destroyDarray(darray *da) {\r
 333    if (!da) return;  // nothing to free\r
 334    \r
 335    if (da->buf) free(da->buf);\r
 336    if (da->lp) free(da->lp);\r
 337    if (da->sl) free(da->sl);\r
 338    if (da->ss) free(da->ss);\r
 339    if (da->p) free(da->p);\r
 340    if (da->rl) free(da->rl);\r
 341    if (da->rm) free(da->rm);\r
 342    if (da->rs) free(da->rs);\r
 343 }\r
 344 \r
 345 int darray_rank0(darray *da, int i)\r
 346 {\r
 347   int r,j;\r
 348   pb *p;\r
 349 \r
 350 #if (RRR == D*2)\r
 351   r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];\r
 352   p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
 353   j = i & (RRR-1);\r
 354   if (j < D) r += popcount(*p >> (D-1-j));\r
 355   else r += popcount(*p) + popcount(p[1] >> (D-1-(j-D)));\r
 356 #else\r
 357 \r
 358   j = i & (RRR-1);\r
 359   if (j < RRR/2) {\r
 360     r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];\r
 361     p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
 362     while (j >= D) {\r
 363       r += popcount(*p++);\r
 364       j -= D;\r
 365     }\r
 366     r += popcount(*p >> (D-1-j));\r
 367   } else {\r
 368     j = RRR-1 - (i & (RRR-1));\r
 369     i += j+1;\r
 370     r = da->rl[i>>logR] + da->rm[i>>logRR] + da->rs[i>>logRRR];\r
 371     p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
 372     while (j >= D) {\r
 373       r -= popcount(*--p);\r
 374       j -= D;\r
 375     }\r
 376     if (j > 0) r -= popcount(*--p << (D-j));\r
 377   }\r
 378 \r
 379 #endif\r
 380 \r
 381   return r;\r
 382 }\r
 383 \r
 384 int darray_rank(darray *da, int i)\r
 385 {\r
 386   int r,j;\r
 387   pb *p;\r
 388 \r
 389   r = da->rl[i>>logR] + da->rm[i>>logRR];\r
 390   j = (i>>logRRR) & (RR/RRR-1);\r
 391   while (j > 0) {\r
 392     r += da->rs[((i>>logRR)<<(logRR-logRRR))+j-1];\r
 393     j--;\r
 394   }\r
 395 \r
 396   p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
 397   j = i & (RRR-1);\r
 398   while (j >= D) {\r
 399     r += popcount(*p++);\r
 400     j -= D;\r
 401   }\r
 402   r += popcount(*p >> (D-1-j));\r
 403 \r
 404   return r;\r
 405 }\r
 406 \r
 407 int darray_select_bsearch(darray *da, int i, pb (*getpat)(pb *))\r
 408 {\r
 409   int j;\r
 410   int l,r,m,n;\r
 411   int s;\r
 412   int t,x,rr;\r
 413   pb *p,w;\r
 414 \r
 415   // for debug\r
 416   //s = darray_select(da,i,1);\r
 417   //\r
 418   //printf("select(%d)=%d\n",i,s);\r
 419 \r
 420 \r
 421 \r
 422   if (i == 0) return -1;\r
 423 \r
 424   if (i > da->m) {\r
 425     return -1;\r
 426   }\r
 427 \r
 428   i--;\r
 429 \r
 430 \r
 431 \r
 432   n = da->m;\r
 433 \r
 434   t = i;\r
 435 \r
 436   l = 0;  r = (n-1)>>logR;\r
 437   // find the smallest index x s.t. rl[x] >= t\r
 438   while (l < r) {\r
 439     m = (l+r)/2;\r
 440     //printf("m=%d rl[m+1]=%d t=%d\n",m,da->rl[m+1],t);\r
 441     if (da->rl[m+1] > t) { // m+1 is out of range\r
 442       r = m;  // new r = m >= l\r
 443     } else {\r
 444       l = m+1; // new l = m+1 <= r\r
 445     }\r
 446   }\r
 447   x = l;\r
 448   t -= da->rl[x];\r
 449 \r
 450   x <<= logR;\r
 451 \r
 452   l = x >> logRR;  r = (min(x+R1-1,n))>>logRR;\r
 453   while (l < r) {\r
 454     m = (l+r)/2;\r
 455     //printf("m=%d rm[m+1]=%d t=%d\n",m,da->rm[m+1],t);\r
 456     if (da->rm[m+1] > t) { // m+1 is out of range\r
 457       r = m;\r
 458     } else {\r
 459       l = m+1; // new l = m+1 <= r\r
 460     }\r
 461   }\r
 462   x = l;\r
 463   t -= da->rm[x];\r
 464 \r
 465   x <<= logRR;\r
 466 \r
 467 #if 0\r
 468   l = x >> logRRR;  r = (min(x+RR-1,n))>>logRRR;\r
 469   while (l < r) {\r
 470     m = (l+r)/2;\r
 471     //printf("m=%d rs[m+1]=%d t=%d\n",m,da->rs[m+1],t);\r
 472     if (da->rs[m+1] > t) { // m+1 is out of range\r
 473       r = m;\r
 474     } else {\r
 475       l = m+1; // new l = m+1 <= r\r
 476     }\r
 477   }\r
 478   x = l;\r
 479   t -= da->rs[x];\r
 480 #else\r
 481   l = x >> logRRR;\r
 482   while (t > da->rs[l]) {\r
 483     t -= da->rs[l];\r
 484     l++;\r
 485   }\r
 486   x = l;\r
 487 #endif\r
 488 \r
 489   x <<= logRRR;\r
 490 \r
 491   p = &da->buf[x >> logD];\r
 492   while (1) {\r
 493     m = popcount(getpat(p));\r
 494     if (m > t) break;\r
 495     t -= m;\r
 496     x += D;\r
 497     p++;\r
 498   }\r
 499 \r
 500   w = getpat(p);\r
 501   while (1) {\r
 502     rr = popCount[w >> (D-8)];\r
 503     if (rr > t) break;\r
 504     t -= rr;\r
 505     x += 8;\r
 506     w <<= 8;\r
 507   }\r
 508   x += selecttbl[((t-0)<<8)+(w>>(D-8))];\r
 509 \r
 510 #if 0\r
 511   if (x != s) {\r
 512     printf("error x=%d s=%d\n",x,s);\r
 513   }\r
 514 #endif\r
 515   return x;\r
 516 }\r
 517 \r
 518 int darray_pat_rank(darray *da, int i, pb (*getpat)(pb *))\r
 519 {\r
 520   int r,j;\r
 521   pb *p;\r
 522 \r
 523   r = da->rl[i>>logR] + da->rm[i>>logRR];\r
 524   j = (i>>logRRR) & (RR/RRR-1);\r
 525   while (j > 0) {\r
 526     r += da->rs[((i>>logRR)<<(logRR-logRRR))+j-1];\r
 527     j--;\r
 528   }\r
 529 \r
 530   p = da->buf + ((i>>logRRR)<<(logRRR-logD));\r
 531   j = i & (RRR-1);\r
 532   while (j >= D) {\r
 533     r += popcount(getpat(p));\r
 534     p++;\r
 535     j -= D;\r
 536   }\r
 537   r += popcount(getpat(p) >> (D-1-j));\r
 538 \r
 539   return r;\r
 540 }\r
 541 \r
 542 \r
 543 int darray_select(darray *da, int i,int f)\r
 544 {\r
 545   int p,r;\r
 546   int il;\r
 547   int rr;\r
 548   pb x;\r
 549   pb *q;\r
 550 \r
 551   if (i == 0) return -1;\r
 552 \r
 553   if (i > da->m) {\r
 554     return -1;\r
 555     //printf("ERROR: m=%d i=%d\n",da->m,i);\r
 556     //exit(1);\r
 557   }\r
 558 \r
 559   i--;\r
 560 \r
 561   il = da->p[i>>logL];\r
 562   if (il < 0) {\r
 563     il = -il-1;\r
 564     p = da->sl[(il<<logL)+(i & (L-1))];\r
 565   } else {\r
 566     p = da->lp[i>>logL];\r
 567     p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];\r
 568     r = i - (i & (LLL-1));\r
 569 \r
 570     q = &(da->buf[p>>logD]);\r
 571 \r
 572     if (f == 1) {\r
 573       rr = p & (D-1);\r
 574       r -= popcount(*q >> (D-1-rr));\r
 575       p = p - rr;\r
 576       \r
 577       while (1) {\r
 578         rr = popcount(*q);\r
 579         if (r + rr >= i) break;\r
 580         r += rr;\r
 581         p += D;\r
 582         q++;\r
 583       }\r
 584       \r
 585       x = *q;\r
 586       while (1) {\r
 587         //rr = popcount(x >> (D-8));\r
 588         rr = popCount[x >> (D-8)];\r
 589         //rr = popcount8(x >> (D-8));\r
 590         if (r + rr >= i) break;\r
 591         r += rr;\r
 592         p += 8;\r
 593         x <<= 8;\r
 594       }\r
 595       p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];\r
 596     } else {\r
 597       rr = p & (D-1);\r
 598       r -= popcount((~(*q))  >> (D-1-rr));\r
 599       p = p - rr;\r
 600       \r
 601       while (1) {\r
 602         rr = popcount(~(*q));\r
 603         if (r + rr >= i) break;\r
 604         r += rr;\r
 605         p += D;\r
 606         q++;\r
 607       }\r
 608       \r
 609       x = ~(*q);\r
 610 \r
 611       while (1) {\r
 612         //rr = popcount(x >> (D-8));\r
 613         rr = popCount[x >> (D-8)];\r
 614         //rr = popcount8(x >> (D-8));\r
 615         if (r + rr >= i) break;\r
 616         r += rr;\r
 617         p += 8;\r
 618         x <<= 8;\r
 619       }\r
 620       p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];\r
 621     }\r
 622   }\r
 623   return p;\r
 624 }\r
 625 \r
 626 int darray_pat_select(darray *da, int i, pb (*getpat)(pb *))\r
 627 {\r
 628   int p,r;\r
 629   int il;\r
 630   int rr;\r
 631   pb x;\r
 632   pb *q;\r
 633 \r
 634   if (i == 0) return -1;\r
 635 \r
 636   if (i > da->m) {\r
 637     return -1;\r
 638     //printf("ERROR: m=%d i=%d\n",da->m,i);\r
 639     //exit(1);\r
 640   }\r
 641 \r
 642   i--;\r
 643 \r
 644   il = da->p[i>>logL];\r
 645   if (il < 0) {\r
 646     il = -il-1;\r
 647     p = da->sl[(il<<logL)+(i & (L-1))];\r
 648   } else {\r
 649     p = da->lp[i>>logL];\r
 650     p += da->ss[(il<<(logL-logLLL))+(i & (L-1))/LLL];\r
 651     r = i - (i & (LLL-1));\r
 652 \r
 653     q = &(da->buf[p>>logD]);\r
 654 \r
 655     rr = p & (D-1);\r
 656     r -= popcount(getpat(q) >> (D-1-rr));\r
 657     p = p - rr;\r
 658     \r
 659     while (1) {\r
 660       rr = popcount(getpat(q));\r
 661       if (r + rr >= i) break;\r
 662       r += rr;\r
 663       p += D;\r
 664       q++;\r
 665     }\r
 666     \r
 667     x = getpat(q);\r
 668     while (1) {\r
 669       //rr = popcount(x >> (D-8));\r
 670       rr = popCount[x >> (D-8)];\r
 671       //rr = popcount8(x >> (D-8));\r
 672       if (r + rr >= i) break;\r
 673       r += rr;\r
 674       p += 8;\r
 675       x <<= 8;\r
 676     }\r
 677     p += selecttbl[((i-r-1)<<8)+(x>>(D-8))];\r
 678   }\r
 679   return p;\r
 680 }\r
 681 \r
 682 int darray_pat_construct(darray *da, int n, pb *buf, int k, pb pat, int opt)\r
 683 {\r
 684   int i;\r
 685   pb *b;\r
 686   mymalloc(b,(n+D-1)/D,0);\r
 687 \r
 688   for (i=0; i<n-k+1; i++) {\r
 689     setbit(b,i,getpattern(buf,i,k,pat));\r
 690   }\r
 691   for (i=n-k+1; i<n; i++) {\r
 692     setbit(b,i,0);\r
 693   }\r
 694 \r
 695   darray_construct(da,n,b,opt);\r
 696   da->buf = buf;\r
 697 \r
 698   free(b);\r
 699   \r
 700   return 0;\r
 701 }\r