lzindex/bitmap.c

   1
   2
   3 // Implements operations over a bitmap
   4
   5 #include "bitmap.h"
   6
   7   // In theory, we should have superblocks of size s=log^2 n divided into
   8   // blocks of size b=(log n)/2. This takes
   9   // O(n log n / log^2 n + n log log n / log n + log n sqrt n log log n) bits
  10   // In practice, we can have any s and b, and the needed amount of bits is
  11   // (n/s) log n + (n/b) log s + b 2^b log b bits
  12   // Optimizing it turns out that s should be exactly s = b log n
  13   // Optimizing b is more difficult but could be done numerically.
  14   // However, the exponential table does no more than popcounting, so why not
  15   // setting up a popcount algorithm tailored to the computer register size,
  16   // defining that size as b, and proceeding.
  17
  18 const unsigned char popc[] =
  19  {
  20     0,1,1,2,1,2,2,3,1,2,2,3,2,3,3,4,1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,
  21     1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
  22     1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
  23     2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
  24     1,2,2,3,2,3,3,4,2,3,3,4,3,4,4,5,2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,
  25     2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
  26     2,3,3,4,3,4,4,5,3,4,4,5,4,5,5,6,3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,
  27     3,4,4,5,4,5,5,6,4,5,5,6,5,6,6,7,4,5,5,6,5,6,6,7,5,6,6,7,6,7,7,8,
  28  };
  29
  30 uint popcount (register uint x)
  31  {
  32     return popc[x&0xFF] + popc[(x>>8)&0xFF] + popc[(x>>16)&0xFF] + popc[x>>24];
  33  }
  34
  35 uint popcount8(register int x)
  36  {
  37     return popc[x & 0xff];
  38  }
  39
  40
  41         // creates a bitmap structure from a bitstring, which is shared
  42
  43 bitmap createBitmap (uint *string, uint n, bool isfullbmap)
  44  {
  45     bitmap B;
  46     uint i,j,pop,bpop,pos;
  47     uint nb,ns,words;
  48     B = malloc (sizeof(struct sbitmap));
  49     B->data = string;
  50     B->n = n; words = (n+W-1)/W;
  51     ns = (n+256-1)/256; nb = 256/W; // adjustments
  52     B->bdata = malloc (ns*nb*sizeof(byte));
  53     B->sdata = malloc (ns*sizeof(int));
  54 #ifdef INDEXREPORT
  55     printf ("     Bitmap over %i bits took %i bits\n", n,n+ns*nb*8+ns*32);
  56 #endif
  57     pop = 0; pos = 0;
  58     for (i=0;i<ns;i++) {
  59        bpop = 0;
  60        B->sdata[i] = pop;
  61        for (j=0;j<nb;j++) {
  62           if (pos == words) break;
  63           B->bdata[pos++] = bpop;
  64           bpop += popcount(*string++);
  65        }
  66        pop += bpop;
  67     }
  68     if (isfullbmap) {// creates the data structure to solve select_0() queries
  69        B->bdata_0 = malloc(ns*nb*sizeof(byte));
  70        B->sdata_0 = malloc(ns*sizeof(int));
  71        string = B->data;
  72        pop = 0; pos = 0;
  73        for (i = 0; i < ns; i++) {
  74           bpop = 0;
  75           B->sdata_0[i] = pop;
  76           for (j = 0; j < nb; j++) {
  77              if (pos == words) break;
  78              B->bdata_0[pos++] = bpop;
  79              bpop += popcount(~(*string++));
  80           }
  81           pop += bpop;
  82        }
  83     }
  84     else { B->bdata_0 = NULL; B->sdata_0 = NULL;}
  85     return B;
  86  }
  87
  88         // rank(i): how many 1's are there before position i, not included
  89
  90 uint rank (bitmap B, uint i)
  91  {
  92     return B->sdata[i>>8] + B->bdata[i>>5] +
  93            popcount (B->data[i>>5] & ((1<<(i&0x1F))-1));
  94  }
  95
  96
  97         // select_1(x): returns the position i of the bitmap B such that
  98         // the number of ones up to position i is x.
  99
 100 uint select_1(bitmap B, uint x)
 101  {
 102     uint s = 256;
 103     uint l = 0, n = B->n, r = n/s,left;
 104     uint mid = (l+r)/2;
 105     uint ones, j;
 106     uint rankmid = B->sdata[mid];
 107     // first, binary search on the superblock array
 108     while (l <= r) {
 109        if (rankmid < x)
 110           l = mid+1;
 111        else
 112           r = mid-1;
 113        mid = (l+r)/2;
 114        rankmid = B->sdata[mid];
 115     }
 116     // sequential search using popcount over an int
 117     left = mid*8;
 118     x -= rankmid;
 119     j = B->data[left];
 120     ones = popcount(j);
 121     while (ones < x) {
 122        x -= ones;
 123        left++;
 124        if (left > (n+W-1)/W) return n;
 125        j = B->data[left];
 126        ones = popcount(j);
 127     }
 128     // sequential search using popcount over a char
 129     left = left*32;
 130     rankmid = popcount8(j);
 131     if (rankmid < x) {
 132        j = j>>8;
 133        x -= rankmid;
 134        left += 8;
 135        rankmid = popcount8(j);
 136        if (rankmid < x) {
 137           j = j>>8;
 138           x -= rankmid;
 139           left += 8;
 140           rankmid = popcount8(j);
 141           if (rankmid < x) {
 142              j = j>>8;
 143              x -= rankmid;
 144              left += 8;
 145           }
 146        }
 147     }
 148   // finally sequential search bit per bit
 149     while (x > 0) {
 150        if  (j&1) x--;
 151        j = j>>1;
 152        left++;
 153     }
 154     return left-1;
 155 }
 156
 157
 158         // select_0(x): returns the position i of the bitmap B such that
 159         // the number of zeros up to position i is x.
 160
 161 uint select_0(bitmap B, uint x)
 162  {
 163     uint s = 256;
 164     uint l = 0, n = B->n, r = n/s, left;
 165     uint mid = (l+r)/2;
 166     uint ones, j;
 167     uint rankmid = B->sdata_0[mid];
 168     // first, binary search on the superblock array
 169     while (l <= r) {
 170        if (rankmid < x)
 171           l = mid+1;
 172        else
 173           r = mid-1;
 174        mid = (l+r)/2;
 175        rankmid = B->sdata_0[mid];
 176     }
 177     // sequential search using popcount over an int
 178     left = mid*8;
 179     x -= rankmid;
 180     j = B->data[left];
 181     ones = popcount(~j);
 182     while (ones < x) {
 183        x -= ones;
 184        left++;
 185        if (left > (n+W-1)/W) return n;
 186        j = B->data[left];
 187        ones = popcount(~j);
 188     }
 189     // sequential search using popcount over a char
 190     left = left*32;
 191     rankmid = popcount8(~j);
 192     if (rankmid < x) {
 193        j = j>>8;
 194        x -= rankmid;
 195        left += 8;
 196        rankmid = popcount8(~j);
 197        if (rankmid < x) {
 198           j = j>>8;
 199           x -= rankmid;
 200           left += 8;
 201           rankmid = popcount8(~j);
 202           if (rankmid < x) {
 203              j = j>>8;
 204              x -= rankmid;
 205              left += 8;
 206           }
 207        }
 208     }
 209   // finally sequential search bit per bit
 210     while (x > 0) {
 211        if  ((j&1)==0) x--;
 212        j = j>>1;
 213        left++;
 214     }
 215     return left-1;
 216  }
 217
 218
 219
 220
 221         // destroys the bitmap, freeing the original bitstream
 222
 223 void destroyBitmap (bitmap B)
 224
 225    { free (B->data);
 226      free (B->bdata);
 227      free (B->sdata);
 228      free (B);
 229    }
 230
 231
 232 uint sizeofBitmap(bitmap B)
 233  {
 234     return sizeof(struct sbitmap) +
 235            ((B->n+W-1)/W)*sizeof(uint) + // data
 236            ((B->n+256-1)/256)*sizeof(int) + // sdata
 237            ((B->n+256-1)/256)*(256/W)*sizeof(byte); // bdata
 238            //((B->sdata_0)?((B->n+256-1)/256)*sizeof(int):0) +
 239            //((B->bdata_0)?((B->n+256-1)/256)*(256/W)*sizeof(byte):0);
 240  }
 241