libcds/src/static_bitsequence/static_bitsequence_rrr02_light.cpp

   1 /* static_bitsequence_rrr02_light.cpp
   2  * Copyright (C) 2008, Francisco Claude, all rights reserved.
   3  *
   4  * static_bitsequence_rrr02_light definition
   5  *
   6  * This library is free software; you can redistribute it and/or
   7  * modify it under the terms of the GNU Lesser General Public
   8  * License as published by the Free Software Foundation; either
   9  * version 2.1 of the License, or (at your option) any later version.
  10  *
  11  * This library is distributed in the hope that it will be useful,
  12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * Lesser General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU Lesser General Public
  17  * License along with this library; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  19  *
  20  */
  21
  22 #include <static_bitsequence_rrr02_light.h>
  23
  24 #define VARS_NEEDED uint C_len = len/BLOCK_SIZE_LIGHT + (len%BLOCK_SIZE_LIGHT!=0);\
  25 uint C_field_bits = bits(BLOCK_SIZE_LIGHT);\
  26 uint O_len = uint_len(1,O_bits_len);\
  27 uint C_sampling_len = C_len/sample_rate+2;\
  28 uint C_sampling_field_bits = bits(ones);\
  29 uint O_pos_len = C_len/sample_rate+1;\
  30 uint O_pos_field_bits = bits(O_bits_len);
  31
  32
  33 table_offset * static_bitsequence_rrr02_light::E = NULL;
  34
  35 static_bitsequence_rrr02_light::static_bitsequence_rrr02_light() {
  36   ones=0;
  37   len=0;
  38   if(E==NULL) E = new table_offset(BLOCK_SIZE_LIGHT);
  39   E->use();
  40   C = NULL;
  41   O = NULL;
  42   C_sampling = NULL;
  43   O_pos = NULL;
  44   sample_rate = DEFAULT_SAMPLING_LIGHT;
  45   O_bits_len = 0;
  46 }
  47
  48 static_bitsequence_rrr02_light::static_bitsequence_rrr02_light(uint * bitseq, uint len, uint sample_rate) {
  49   ones = 0;
  50   this->len = len;
  51   if(E==NULL) E = new table_offset(BLOCK_SIZE_LIGHT);
  52   E->use();
  53   // Table C
  54   uint C_len = len/BLOCK_SIZE_LIGHT + (len%BLOCK_SIZE_LIGHT!=0);
  55   uint C_field_bits = bits(BLOCK_SIZE_LIGHT);
  56   C = new uint[uint_len(C_len,C_field_bits)];
  57   for(uint i=0;i<uint_len(C_len,C_field_bits);i++)
  58     C[i] = 0;
  59   O_bits_len = 0;
  60   for(uint i=0;i<C_len;i++) {
  61     uint value = popcount(get_var_field(bitseq,i*BLOCK_SIZE_LIGHT,min((uint)len-1,(i+1)*BLOCK_SIZE_LIGHT-1)));
  62     assert(value<=BLOCK_SIZE_LIGHT);
  63     set_field(C,C_field_bits,i,value);
  64     ones += value;
  65     O_bits_len += E->get_log2binomial(BLOCK_SIZE_LIGHT,value);
  66   }
  67   // Table O
  68   uint O_len = uint_len(1,O_bits_len);
  69   O = new uint[O_len];
  70   for(uint i=0;i<O_len;i++)
  71     O[i] = 0;
  72   uint O_pos = 0;
  73   for(uint i=0;i<C_len;i++) {
  74     uint value = (ushort)get_var_field(bitseq,i*BLOCK_SIZE_LIGHT,min((uint)len-1,(i+1)*BLOCK_SIZE_LIGHT-1));
  75     set_var_field(O,O_pos,O_pos+E->get_log2binomial(BLOCK_SIZE_LIGHT,popcount(value))-1,E->compute_offset((ushort)value));
  76     O_pos += E->get_log2binomial(BLOCK_SIZE_LIGHT,popcount(value));
  77   }
  78   C_sampling = NULL;
  79   this->O_pos = NULL;
  80
  81   create_sampling(sample_rate);
  82 }
  83
  84 void static_bitsequence_rrr02_light::create_sampling(uint sample_rate) {
  85   this->sample_rate = sample_rate;
  86 /*  for(uint i=0;i<C_len;i++) {
  87     O_bits_len += E->get_log2binomial(BLOCK_SIZE_LIGHT,get_field(C,C_field_bits,i));
  88   }*/
  89   // Sampling for C
  90   uint C_len = len/BLOCK_SIZE_LIGHT + (len%BLOCK_SIZE_LIGHT!=0);
  91   uint C_field_bits = bits(BLOCK_SIZE_LIGHT);
  92   uint C_sampling_len = C_len/sample_rate+2;
  93   uint C_sampling_field_bits = bits(ones);
  94   if(C_sampling!=NULL) delete [] C_sampling;
  95   C_sampling = new uint[max((uint)1,uint_len(C_sampling_len,C_sampling_field_bits))];
  96   for(uint i=0;i<max((uint)1,uint_len(C_sampling_len,C_sampling_field_bits));i++)
  97     C_sampling[i] = 0;
  98   uint sum = 0;
  99   for(uint i=0;i<C_len;i++) {
 100     if(i%sample_rate==0)
 101       set_field(C_sampling,C_sampling_field_bits,i/sample_rate,sum);
 102     sum += get_field(C,C_field_bits,i);
 103   }
 104   for(uint i=(C_len-1)/sample_rate+1;i<C_sampling_len;i++)
 105     set_field(C_sampling,C_sampling_field_bits,i,sum);
 106   // Sampling for O (table S) (Code separated from previous construction for readability)
 107   uint O_pos_len = C_len/sample_rate+1;
 108   uint O_pos_field_bits = bits(O_bits_len);
 109   if(O_pos!=NULL) delete [] O_pos;
 110   O_pos = new uint[uint_len(O_pos_len,O_pos_field_bits)];
 111   for(uint i=0;i<uint_len(O_pos_len,O_pos_field_bits);i++)
 112     O_pos[i] = 0;
 113   uint pos = 0;
 114   for(uint i=0;i<C_len;i++) {
 115     if(i%sample_rate==0)
 116       set_field(O_pos,O_pos_field_bits,i/sample_rate,pos);
 117     pos += E->get_log2binomial(BLOCK_SIZE_LIGHT,get_field(C,C_field_bits,i));
 118   }
 119 }
 120
 121 bool static_bitsequence_rrr02_light::access(uint i) {
 122   uint C_len = len/BLOCK_SIZE_LIGHT + (len%BLOCK_SIZE_LIGHT!=0);
 123   uint C_field_bits = bits(BLOCK_SIZE_LIGHT);
 124   uint O_pos_field_bits = bits(O_bits_len);
 125   uint nearest_sampled_value = i/BLOCK_SIZE_LIGHT/sample_rate;
 126   uint pos_O = get_field(O_pos,O_pos_field_bits,nearest_sampled_value);
 127   uint pos = i/BLOCK_SIZE_LIGHT;
 128   assert(pos<=C_len);
 129   for(uint k=nearest_sampled_value*sample_rate;k<pos;k++) {
 130     uint aux = get_field(C,C_field_bits,k);
 131     pos_O += E->get_log2binomial(BLOCK_SIZE_LIGHT,aux);
 132   }
 133   uint c = get_field(C,C_field_bits,pos);
 134   return ((1<<(i%BLOCK_SIZE_LIGHT))&E->short_bitmap(c,get_var_field(O,pos_O,pos_O+E->get_log2binomial(BLOCK_SIZE_LIGHT,c)-1)))!=0;
 135 }
 136
 137 uint static_bitsequence_rrr02_light::rank0(uint i) {
 138   if(i+1==0) return 0;
 139   return 1+i-rank1(i);
 140 }
 141
 142 uint static_bitsequence_rrr02_light::rank1(uint i) {
 143   uint C_field_bits = bits(BLOCK_SIZE_LIGHT);
 144   uint C_sampling_field_bits = bits(ones);
 145   uint O_pos_field_bits = bits(O_bits_len);
 146   if(i+1==0) return 0;
 147   uint nearest_sampled_value = i/BLOCK_SIZE_LIGHT/sample_rate;
 148   uint sum = get_field(C_sampling,C_sampling_field_bits,nearest_sampled_value);
 149   uint pos_O = get_field(O_pos,O_pos_field_bits,nearest_sampled_value);
 150   uint pos = i/BLOCK_SIZE_LIGHT;
 151   uint k=nearest_sampled_value*sample_rate;
 152   if(k%2==1 && k<pos) {
 153     uint aux = get_field(C,C_field_bits,k);
 154     sum += aux;
 155     pos_O += E->get_log2binomial(BLOCK_SIZE_LIGHT,aux);
 156     k++;
 157   }
 158   uchar * a = (uchar *)C;
 159   uint mask = 0x0F;
 160   a += k/2;
 161   while(k<(uint)max(0,(int)pos-1)) {
 162     assert(((*a)&mask)==get_field(C,C_field_bits,k));
 163     assert((*a)/16==get_field(C,C_field_bits,k+1));
 164     sum += ((*a)&mask)+(*a)/16;
 165     pos_O += E->get_log2binomial(BLOCK_SIZE_LIGHT,((*a)&mask))+E->get_log2binomial(BLOCK_SIZE_LIGHT,((*a)/16));
 166     a++;
 167     k+=2;
 168   }
 169   if(k<pos) {
 170     uint aux = get_field(C,C_field_bits,k);
 171     sum += aux;
 172     pos_O += E->get_log2binomial(BLOCK_SIZE_LIGHT,aux);
 173     k++;
 174   }
 175   uint c = get_field(C,C_field_bits,pos);
 176   sum += popcount(((2<<(i%BLOCK_SIZE_LIGHT))-1) & E->short_bitmap(c,get_var_field(O,pos_O,pos_O+E->get_log2binomial(BLOCK_SIZE_LIGHT,c)-1)));
 177   return sum;
 178 }
 179
 180 uint static_bitsequence_rrr02_light::select0(uint i) {
 181   uint C_len = len/BLOCK_SIZE_LIGHT + (len%BLOCK_SIZE_LIGHT!=0);
 182   uint C_field_bits = bits(BLOCK_SIZE_LIGHT);
 183   uint C_sampling_len = C_len/sample_rate+2;
 184   uint C_sampling_field_bits = bits(ones);
 185   uint O_pos_field_bits = bits(O_bits_len);
 186   if(i==0) return -1;
 187   if(i>len-ones) return len;
 188   // Search over partial sums
 189   uint start=0;
 190   uint end=C_sampling_len-1;
 191   uint med, acc=0, pos;
 192   while(start<end-1) {
 193     med = (start+end)/2;
 194     acc = med*sample_rate*BLOCK_SIZE_LIGHT-get_field(C_sampling,C_sampling_field_bits,med);
 195     if(acc<i) {
 196       if(med==start) break;
 197       start=med;
 198     }
 199     else  {
 200       if(end==0) break;
 201       end = med-1;
 202     }
 203   }
 204   acc = get_field(C_sampling,C_sampling_field_bits,start);
 205   while(start<C_len-1 && acc+sample_rate*BLOCK_SIZE_LIGHT==get_field(C_sampling,C_sampling_field_bits,start+1)) {
 206     start++;
 207     acc +=sample_rate*BLOCK_SIZE_LIGHT;
 208   }
 209   acc = start*sample_rate*BLOCK_SIZE_LIGHT-acc;
 210   pos = (start)*sample_rate;
 211   uint pos_O = get_field(O_pos,O_pos_field_bits,start);
 212   // Sequential search over C
 213   uint s = 0;
 214   for(;pos<C_len;pos++) {
 215     s = get_field(C,C_field_bits,pos);
 216     if(acc+BLOCK_SIZE_LIGHT-s>=i) break;
 217     pos_O += E->get_log2binomial(BLOCK_SIZE_LIGHT,s);
 218     acc += BLOCK_SIZE_LIGHT-s;
 219   }
 220   pos = (pos)*BLOCK_SIZE_LIGHT;
 221   // Search inside the block
 222
 223   while(acc<i) {
 224     uint new_posO = pos_O+E->get_log2binomial(BLOCK_SIZE_LIGHT,s);
 225     uint block = E->short_bitmap(s,get_var_field(O,pos_O,new_posO-1));
 226     pos_O = new_posO;
 227     new_posO = 0;
 228     while(acc<i && new_posO<BLOCK_SIZE_LIGHT) {
 229       pos++;new_posO++;
 230       acc += (((block&1)==0)?1:0);
 231       block = block/2;
 232     }
 233   }
 234   pos--;
 235   assert(acc==i);
 236   assert(rank0(pos)==i);
 237   assert(!access(pos));
 238   return pos;
 239 }
 240
 241 uint static_bitsequence_rrr02_light::select1(uint i) {
 242   uint C_len = len/BLOCK_SIZE_LIGHT + (len%BLOCK_SIZE_LIGHT!=0);
 243   uint C_field_bits = bits(BLOCK_SIZE_LIGHT);
 244   uint C_sampling_len = C_len/sample_rate+2;
 245   uint C_sampling_field_bits = bits(ones);
 246   uint O_pos_field_bits = bits(O_bits_len);
 247   if(i==0) return -1;
 248   if(i>ones) return len;
 249   // Search over partial sums
 250   uint start=0;
 251   uint end=C_sampling_len-1;
 252   uint med, acc=0, pos;
 253   while(start<end-1) {
 254     med = (start+end)/2;
 255     acc = get_field(C_sampling,C_sampling_field_bits,med);
 256     if(acc<i) {
 257       if(med==start) break;
 258       start=med;
 259     }
 260     else  {
 261       if(end==0) break;
 262       end = med-1;
 263     }
 264   }
 265   acc = get_field(C_sampling,C_sampling_field_bits,start);
 266   while(start<C_len-1 && acc==get_field(C_sampling,C_sampling_field_bits,start+1)) start++;
 267   pos = (start)*sample_rate;
 268   uint pos_O = get_field(O_pos,O_pos_field_bits,start);
 269   acc = get_field(C_sampling,C_sampling_field_bits,start);
 270   // Sequential search over C
 271   uint s = 0;
 272   for(;pos<C_len;pos++) {
 273     s = get_field(C,C_field_bits,pos);
 274     if(acc+s>=i) break;
 275     pos_O += E->get_log2binomial(BLOCK_SIZE_LIGHT,s);
 276     acc += s;
 277   }
 278   pos = (pos)*BLOCK_SIZE_LIGHT;
 279   //cout << "pos=" << pos << endl;
 280   // Search inside the block
 281   while(acc<i) {
 282     uint new_posO = pos_O+E->get_log2binomial(BLOCK_SIZE_LIGHT,s);
 283     uint block = E->short_bitmap(s,get_var_field(O,pos_O,new_posO-1));
 284     pos_O = new_posO;
 285     new_posO = 0;
 286     while(acc<i && new_posO<BLOCK_SIZE_LIGHT) {
 287       pos++;new_posO++;
 288       acc += (((block&1)!=0)?1:0);
 289       block = block/2;
 290     }
 291     //cout << "i=" << i << " len=" << len << " ones=" << ones << " pos=" << pos << " acc=" << acc << " rank=" << rank1(pos) << endl;
 292   }
 293   pos--;
 294   assert(acc==i);
 295   assert(rank1(pos)==i);
 296   assert(access(pos));
 297   return pos;
 298 }
 299
 300 uint static_bitsequence_rrr02_light::size() {
 301   VARS_NEEDED
 302   /*cout << "RRR02 SIZE: " << endl;
 303   cout << "Default: " << 9*sizeof(uint)+sizeof(uint*)*4 << endl;
 304   cout << "Cs:      " << uint_len(C_len,C_field_bits)*sizeof(uint) << endl;
 305   cout << "Os:      " << O_len*sizeof(uint) << endl;
 306   cout << "CSamp:   " << uint_len(C_sampling_len,C_sampling_field_bits)*sizeof(uint) << endl;
 307   cout << "OSamp:   " << uint_len(O_pos_len,O_pos_field_bits)*sizeof(uint) << endl;
 308   cout << "E:       " << E->size() << endl;*/
 309   uint sum = sizeof(uint)*8;//sizeof(static_bitsequence_rrr02_light);
 310   sum += uint_len(C_len,C_field_bits)*sizeof(uint);
 311   sum += O_len*sizeof(uint);
 312   sum += uint_len(C_sampling_len,C_sampling_field_bits)*sizeof(uint);
 313   sum += uint_len(O_pos_len,O_pos_field_bits)*sizeof(uint);
 314   //sum += E->size();
 315   return sum;
 316 }
 317
 318 static_bitsequence_rrr02_light::~static_bitsequence_rrr02_light() {
 319   if(C!=NULL) delete [] C;
 320   if(O!=NULL) delete [] O;
 321   if(C_sampling!=NULL) delete [] C_sampling;
 322   if(O_pos!=NULL) delete [] O_pos;
 323   E = E->unuse();
 324 }
 325
 326 int static_bitsequence_rrr02_light::save(FILE * fp) {
 327   uint C_len = len/BLOCK_SIZE_LIGHT + (len%BLOCK_SIZE_LIGHT!=0);
 328   uint C_field_bits = bits(BLOCK_SIZE_LIGHT);
 329   uint O_len = uint_len(1,O_bits_len);
 330   uint wr = RRR02_LIGHT_HDR;
 331   wr = fwrite(&wr,sizeof(uint),1,fp);
 332   wr += fwrite(&len,sizeof(uint),1,fp);
 333   wr += fwrite(&ones,sizeof(uint),1,fp);
 334   wr += fwrite(&O_bits_len,sizeof(uint),1,fp);
 335   wr += fwrite(&sample_rate,sizeof(uint),1,fp);
 336   if(wr!=5) return -1;
 337   wr = fwrite(C,sizeof(uint),uint_len(C_len,C_field_bits),fp);
 338   if(wr!=uint_len(C_len,C_field_bits)) return -1;
 339   wr = fwrite(O,sizeof(uint),O_len,fp);
 340   if(wr!=O_len) return -1;
 341   return 0;
 342 }
 343
 344 static_bitsequence_rrr02_light * static_bitsequence_rrr02_light::load(FILE * fp) {
 345   static_bitsequence_rrr02_light * ret = new static_bitsequence_rrr02_light();
 346   uint rd = 0, type;
 347   rd += fread(&type,sizeof(uint),1,fp);
 348   rd += fread(&ret->len,sizeof(uint),1,fp);
 349   rd += fread(&ret->ones,sizeof(uint),1,fp);
 350   rd += fread(&ret->O_bits_len,sizeof(uint),1,fp);
 351   rd += fread(&ret->sample_rate,sizeof(uint),1,fp);
 352   uint C_len = ret->len/BLOCK_SIZE_LIGHT + (ret->len%BLOCK_SIZE_LIGHT!=0);
 353   uint C_field_bits = bits(BLOCK_SIZE_LIGHT);
 354   uint O_len = uint_len(1,ret->O_bits_len);
 355   if(rd!=5 || type!=RRR02_LIGHT_HDR) {
 356     delete ret;
 357     return NULL;
 358   }
 359   ret->C = new uint[uint_len(C_len,C_field_bits)];
 360   rd = fread(ret->C,sizeof(uint),uint_len(C_len,C_field_bits),fp);
 361   if(rd!=uint_len(C_len,C_field_bits)) {
 362     ret->C=NULL;
 363     delete ret;
 364     return NULL;
 365   }
 366   ret->O = new uint[O_len];
 367   rd = fread(ret->O,sizeof(uint),O_len,fp);
 368   if(rd!=O_len) {
 369     ret->O=NULL;
 370     delete ret;
 371     return NULL;
 372   }
 373   ret->create_sampling(ret->sample_rate);
 374   return ret;
 375 }