darray.c

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