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