bp-darray.c

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