bp-darray.c

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