swcsa/utils/bitmap.c

   1
   2 // Implements operations over a bitmap
   3
   4 #include "bitmap.h"
   5
   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 //unsigned char OnesInByte[] =
  19 const unsigned char popc[] =   //number of ones in one byte value [0..255].
  20 {
  21 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,
  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 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,
  24 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,
  25 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,
  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 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,
  28 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,
  29 };
  30
  31 uint popcount (register uint x)
  32
  33    { return popc[x&0xFF] + popc[(x>>8)&0xFF] + popc[(x>>16)&0xFF] + popc[x>>24];
  34    }
  35
  36
  37 /******************************************************************/
  38 // FUNCIONS DE EDU ...
  39 /******************************************************************/
  40 /*
  41         Creates a bitmap and structures to rank and select
  42 */
  43
  44 //bitmap createBitmapEdu (uint *string, uint n){  return createBitmap(string,n);}
  45
  46 bitmap createBitmap (uint *string, uint n){
  47     bitmap B;
  48
  49         unsigned int nb;
  50         unsigned int ns;
  51         unsigned int countB, countS, blockIndex, superblockIndex;
  52     register unsigned int block;
  53
  54         B =(struct sbitmap *) malloc (sizeof(struct sbitmap));
  55         B->data = string;
  56         B->n = n;
  57         ns = (n/256)+1;
  58         nb = (n/32)+1;
  59
  60     B->bSize = nb;
  61     B->sSize = ns;
  62         B->bdata =(byte *)malloc(nb*sizeof(byte));  //  Db = (unsigned char *)malloc(nb*sizeof(unsigned char));
  63         B->sdata = (uint*)malloc(ns*sizeof(uint));  //  Ds = (unsigned int *)malloc(ns*sizeof(unsigned int));
  64
  65         B->mem_usage = (ns*sizeof(uint)) + (nb*sizeof(byte)) + (sizeof(struct sbitmap));
  66         /* Ahora construimos los bloques */
  67     blockIndex = 0;
  68     superblockIndex = 0;
  69     countB = 0;
  70     countS = 0;
  71
  72     while(blockIndex < nb) {
  73
  74        if(!(blockIndex%8)) {
  75           countS += countB;
  76           B->sdata[superblockIndex++] = countS;
  77           countB = 0;
  78        }
  79
  80        B->bdata[blockIndex] = countB;
  81            block = string[blockIndex++];
  82
  83            countB += popcount(block);
  84     }
  85
  86         B->pop = countS+countB;
  87
  88 //      {int i;     //fprintf(stderr,"\n");
  89 //     for (i=0;i<ns;i++) {//fprintf(stderr,"%d ",B->sdata[i]);
  90 //      }
  91 //     //fprintf(stderr,"\n");
  92 //     for (i=0;i<8;i++) {//fprintf(stderr,"%d ",B->bdata[i]);
  93 //      }
  94 //      }
  95     return B;
  96 }
  97
  98
  99 /*
 100   Number of 1s in range [0,posicion]
 101 */
 102 //uint rank1Edu(bitmap B, unsigned int position) {
 103 //uint rank1Edu(bitmap B, unsigned int position) { return rank(B,position);}
 104 uint rank(bitmap B, unsigned int position) {
 105     register unsigned int block;
 106     if (position > B->n) return B->pop;
 107         //position -=1;
 108
 109         block = B->data[position/32] << (31-position%32);
 110
 111     return B->sdata[position/256] + B->bdata[position/32] +
 112            popc[block & 0xff] + popc[(block>>8) & 0xff] +
 113            popc[(block>>16) & 0xff] + popc[block>>24];
 114 }
 115
 116
 117 /********************************************************************************************/
 118 /**************************************************************************************/
 119
 120 static uint binsearch (uint *data, uint size, uint val)
 121
 122    { uint i,j,m;
 123      i = 0; j = size;
 124      while (i+1 < j)
 125         { m = (i+j)/2;
 126           if (data[m] >= val) j = m;
 127           else i = m;
 128         }
 129      return i;
 130    }
 131
 132 uint bselect (bitmap B, uint j)
 133
 134    { uint spos,bpos,pos,word,x;
 135      byte *blk;
 136      if (j > B->pop) return B->n;
 137      spos = binsearch(B->sdata,(B->n+256-1)/256,j);
 138
 139      //fprintf(stderr,"\n SPOS IS %d, and B->sdata[pos] = %d",spos,B->sdata[spos]);
 140      j -= B->sdata[spos];
 141      pos = spos<<8;
 142      blk = B->bdata + (pos>>5);
 143      bpos = 0;
 144
 145     //while ((bpos < (1<<3)-1) && (blk[bpos+1] < j)) bpos++;
 146     while ( ((spos*8+bpos) < ((B->n-1)/W)) && (bpos < (1<<3)-1) && (blk[bpos+1] < j)) bpos++;
 147
 148
 149       //fprintf(stderr,"\n BPOS  = %d",bpos);
 150      pos += bpos<<5;
 151      word = B->data[pos>>5];
 152      j -= blk[bpos];
 153      //fprintf(stderr,"\n pos>>5 = %d ... pasou XXX con word = %d, and j= %d",pos>>5,word,j);
 154      while (1)
 155     { x = popc[word & ((1<<8)-1)];
 156         //fprintf(stderr,"\n word = %u popc vale %u",word & ((1<<8)-1),x);
 157           if (j <= x) break;
 158           j -= x; pos += 8;
 159           word >>= 8;
 160
 161         }
 162
 163      while (j) { if (word & 1) j--; word >>= 1; pos++; }
 164
 165      // fprintf(stderr,"\n\nBSELECT::: POSICIÓN FINAL = %u",pos-1);
 166      return pos-1;
 167
 168    }
 169
 170
 171 // destroys the bitmap, freeing the original bitstream
 172 void destroyBitmap (bitmap B)
 173
 174    { //free (B->data);
 175      free (B->bdata);
 176      free (B->sdata);
 177      free (B);
 178    }
 179
 180
 181 // Prints the bit vector
 182 void showBitVector(uint * V, int vectorSize) {
 183      uint bitIndex=0;
 184      while(bitIndex<vectorSize) {
 185         fprintf(stderr,"%d",bitget(V,bitIndex));
 186         bitIndex++;
 187      }
 188 }
 189
 190 void saveBitmap (char *filename, bitmap b) {
 191         int file;
 192         unlink(filename);
 193         if( (file = open(filename, O_WRONLY|O_CREAT,S_IRWXG | S_IRWXU)) < 0) {
 194                 printf("Cannot open file %s\n", filename);
 195                 exit(0);
 196         }
 197         write(file, &(b->sSize), sizeof(uint));
 198         write(file, b->sdata, sizeof(int) * (b->sSize));
 199         write(file, &(b->bSize), sizeof(uint));
 200         write(file, b->bdata, sizeof(byte) * (b->bSize));
 201
 202         write(file, &(b->pop), sizeof(uint));
 203         write(file, &(b->n), sizeof(uint));
 204         close(file);
 205 }
 206
 207 /* loads the Rank structures from disk, and sets Bitmap->data ptr to "string"
 208 */
 209 bitmap loadBitmap (char *filename, uint *string, uint n) {
 210         bitmap B;
 211         int file;
 212
 213         if( (file = open(filename, O_RDONLY)) < 0) {
 214                 printf("Cannot read file %s\n", filename);
 215                 exit(0);
 216         }
 217
 218         B = (struct sbitmap *) malloc (sizeof(struct sbitmap));
 219         B->data = string;
 220
 221         read(file, &(B->sSize), sizeof(uint));
 222         B->sdata = (uint *) malloc(sizeof(uint) * B->sSize);
 223         read(file, B->sdata, sizeof(uint) * B->sSize);
 224
 225         read(file, &(B->bSize), sizeof(uint));
 226         B->bdata = (byte *) malloc(sizeof(byte) * B->bSize);
 227         read(file, B->bdata, sizeof(byte) * B->bSize);
 228
 229         read(file, &(B->pop), sizeof(uint));
 230         read(file, &(B->n), sizeof(uint));
 231         close(file);
 232         B->mem_usage = (sizeof(uint) * B->sSize) + (sizeof(byte) * B->bSize) + (sizeof(struct sbitmap));
 233
 234         if (n != B->n) {printf("\n LoadBitmap failed: %u distinto de %u",n,B->n); exit(0);}
 235         return B;
 236
 237 }
 238
 239
 240 /********************************************************************************************/
 241 /********************************************************************************************/
 242
 243
 244
 245
 246         // creates a bitmap structure from a bitstring, which is shared
 247
 248 bitmap createBitmapGONZA (uint *string, uint n)
 249 //bitmap createBitmap (uint *string, uint n)
 250
 251    { bitmap B;
 252      uint i,j,pop,bpop,pos;
 253      uint s,nb,ns,words;
 254      B = (struct sbitmap *) malloc (sizeof(struct sbitmap));
 255      B->data = string;
 256
 257
 258          B->n = n; words = (n+W-1)/W;
 259      ns = (n+256-1)/256; nb = 256/W; // adjustments
 260
 261      B->bSize = ns*nb;
 262      B->bdata = (byte *) malloc (ns*nb*sizeof(byte));
 263      B->sSize = ns;
 264      B->sdata = (uint *)malloc (ns*sizeof(int));
 265
 266          B->mem_usage = (ns*sizeof(int)) + (ns*nb*sizeof(byte)) + (sizeof(struct sbitmap));
 267 #ifdef INDEXREPORT
 268      printf ("     Bitmap over %i bits took %i bits\n", n,n+ns*nb*8+ns*32);
 269 #endif
 270           //fprintf (stderr,"     Bitmap over %i bits took %i bits\n", n,n+ns*nb*8+ns*32);
 271      pop = 0; pos = 0;
 272      for (i=0;i<ns;i++)
 273         { bpop = 0;
 274           B->sdata[i] = pop;
 275           for (j=0;j<nb;j++)
 276              { if (pos == words) break;
 277                B->bdata[pos++] = bpop;
 278                bpop += popcount(*string++);
 279              }
 280           pop += bpop;
 281         }
 282      B->pop = pop;
 283
 284                 //     //fprintf(stderr,"\n");
 285                 //     for (i=0;i<ns;i++) {//fprintf(stderr,"%d ",B->sdata[i]);
 286                 //      }
 287                 //     //fprintf(stderr,"\n");
 288                 //     for (i=0;i<ns*nb;i++) {//fprintf(stderr,"%d ",B->bdata[i]);
 289                 //      }
 290
 291      return B;
 292    }
 293
 294         // rank(i): how many 1's are there before position i, not included
 295
 296 //uint rank (bitmap B, uint i)
 297 uint rankGONZA (bitmap B, uint i)
 298
 299    {
 300         i++;
 301         if (i > B->n) return B->pop;
 302         return B->sdata[i>>8] + B->bdata[i>>5] +
 303             popcount (B->data[i>>5] & ((1<<(i&0x1F))-1));
 304    }
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316