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