BSGAP.cpp

   1 /**
   2  * Binary-searchable gap encoding scheme (BSGAP)
   3  * Niko Välimäki
   4  *
   5  * Compile: g++ -Wall -o testBSGAP BSGAP.cpp Tools.cpp
   6  *
   7  * Based on:
   8  *
   9  * Ankur Gupta and Wing-Kai Hon and Rahul Shah and Jeffrey Scott Vitter. Compressed data structures:
  10  * Dictionaries and data-aware measures, Theor. Comput. Sci., Volume 387, Issue 3 (November 2007).
  11  */
  12
  13 #include "BSGAP.h"
  14 #include <stdexcept>
  15 #include <cassert>
  16 using namespace std;
  17
  18
  19 const uchar BSGAP::bit_table[] = {
  20     0,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,2,0,
  21     1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,6,0,1,0,2,0,1,0,3,0,1,0,
  22     2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,
  23     1,0,2,0,1,0,3,0,1,0,2,0,1,0,7,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,
  24     3,0,1,0,2,0,1,0,5,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,
  25     1,0,6,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1,0,5,0,1,0,
  26     2,0,1,0,3,0,1,0,2,0,1,0,4,0,1,0,2,0,1,0,3,0,1,0,2,0,1};
  27
  28 /**
  29  * samplerate == number of keys in each gap encoded binary search tree
  30  * default value is log^2 n where n is the number of keys.
  31  */
  32 BSGAP::BSGAP(ulong *B, ulong u_, bool freeB, ulong sampleRate)
  33     : u(u_), n(0), topCount(0), bitsInP(0), P(0), b(0), offsetA(0), firstKeyA(0)
  34 {
  35     // Count keys
  36     for (ulong i = 0; i < u; i ++)
  37         if (Tools::GetField(B, 1, i))
  38             this->n ++;
  39
  40     if (n == 0) // Sanity check
  41     {
  42         if (freeB)
  43             delete [] B;
  44
  45         return;
  46     }
  47
  48     unsigned log2n = Tools::FloorLog2(this->n);
  49     if (log2n < 2)
  50         log2n = 2;
  51     if (sampleRate == 0)
  52         this->b = log2n * log2n;
  53     else
  54         this->b = sampleRate;
  55
  56     ulong firstKey = 0, lastKey;
  57     // Find the first key
  58     while (!Tools::GetField(B, 1, firstKey))
  59         firstKey ++;
  60
  61     // Temp array of BSGAP structures
  62     ulong **tempP = new ulong * [n / b + 1];
  63     ulong *bitsInBSD = new ulong [n / b + 1];
  64     ulong totalBits = 0;
  65     this->firstKeyA = new BlockArray(n/b + 1, Tools::CeilLog2(this->u));
  66
  67     for (ulong i = 0; i < n/b; i ++)
  68     {
  69         // Find the last key
  70         lastKey = firstKey;
  71         for (ulong k = 0; k < b + 1 && lastKey < u; lastKey ++)
  72             if (Tools::GetField(B, 1, lastKey))
  73                 k ++;
  74         lastKey --;
  75         // Sanity check
  76         if (lastKey > u)
  77         {
  78             cout << "error: lastKey = " << lastKey << ", u = " << u << endl;
  79             exit(0);
  80         }
  81
  82         // lastKey of the last substructure is this->u
  83         if (i == n/b - 1 && n%b == 0)
  84             lastKey = u;
  85
  86         tempP[i] = GetSubtree(B, firstKey, firstKey, lastKey, b, true, bitsInBSD[i]);
  87         totalBits += bitsInBSD[i];
  88         (*firstKeyA)[i] = firstKey;
  89
  90         firstKey = lastKey;
  91     }
  92
  93     if (n - (n/b) * b != 0)
  94     {
  95         // Find the first key
  96         while (!Tools::GetField(B, 1, firstKey))
  97             firstKey ++;
  98
  99         tempP[n/b] = GetSubtree(B, firstKey, firstKey, u, n - (n/b) * b, true, bitsInBSD[n/b]);
 100         totalBits += bitsInBSD[n/b];
 101         (*firstKeyA)[n/b] = firstKey;
 102     }
 103
 104     // Number of nodes:
 105     this->topCount = n%b ? n/b + 1: n/b;
 106
 107     // Catenate binary trees into one bitvector
 108     this->P = new ulong[totalBits/W + 1];
 109     this->bitsInP = totalBits;
 110     ulong offset = 0;
 111     for (ulong i = 0; i < this->topCount ; ++i)
 112     {
 113         BitCopy(P, tempP[i], offset, bitsInBSD[i]);
 114         delete [] tempP[i];
 115         offset += bitsInBSD[i];
 116     }
 117     delete [] tempP;
 118
 119     // Pointers to bit vector P  (binary tree starting offsets)
 120     // FIXME Use more succinct representation?
 121     offset = 0;
 122     this->offsetA = new BlockArray(this->topCount, Tools::CeilLog2(totalBits));
 123     for (ulong i = 0; i < this->topCount ; ++i)
 124     {
 125         (*offsetA)[i] = offset;
 126         offset += bitsInBSD[i];
 127     }
 128     delete [] bitsInBSD;
 129
 130     if (freeB)
 131         delete [] B;
 132 }
 133
 134 BSGAP::~BSGAP()
 135 {
 136         delete [] P;
 137         delete offsetA;
 138         delete firstKeyA;
 139 }
 140
 141 ulong BSGAP::rank(ulong i)
 142 {
 143     if (n == 0) // Trivial case
 144         return 0;
 145     if (i >= u)
 146         throw std::out_of_range("BSGAP::rank(): Parameter i was out of range");
 147     if ((*firstKeyA)[0] > i)
 148         return 0;
 149     ulong l = 0, r = topCount - 1;
 150     while (l < r)   // Binary search on the array firstKeyA
 151     {
 152         ulong mid = l + (r - l)/2lu;
 153         if ((*firstKeyA)[mid] < i)
 154             l = mid + 1;
 155         else
 156             r = mid;
 157     }
 158     if ((*firstKeyA)[l] > i && l == 0)
 159         return 0;
 160     if ((*firstKeyA)[l] == i)
 161         return 1 + l * b;
 162     if ((*firstKeyA)[l] > i)
 163         l --;
 164
 165     ulong lastKey, firstKey = (*firstKeyA)[l];
 166     if (l == topCount - 1)
 167         lastKey = u;
 168     else
 169         lastKey = (*firstKeyA)[l + 1];
 170     ulong nSub = b;
 171     if (l == topCount - 1 && n%b)
 172         nSub = n%b;
 173
 174     ulong k = rank(P, (*offsetA)[l], firstKey, i, firstKey, lastKey, nSub, true);
 175     // Add rank of substructures P[0], P[1], ..., P[j - 2]
 176     return l * b + k;
 177 }
 178
 179 ulong BSGAP::rank0(ulong i)
 180 {
 181     return i - rank(i) + 1;
 182 }
 183
 184 bool BSGAP::IsBitSet(ulong i)
 185 {
 186     if (n == 0) // Trivial case
 187         return false;
 188     if (i >= u)
 189         return false; // FIXME throw std::out_of_range("BSGAP::rank(): Parameter i was out of range");
 190     if ((*firstKeyA)[0] > i)
 191         return false;
 192     ulong l = 0, r = topCount - 1;
 193     while (l < r)   // Binary search on the array firstKeyA
 194     {
 195         ulong mid = l + (r - l)/2lu;
 196         if ((*firstKeyA)[mid] < i)
 197             l = mid + 1;
 198         else
 199             r = mid;
 200     }
 201     if ((*firstKeyA)[l] > i && l == 0)
 202         return false;
 203     if ((*firstKeyA)[l] == i)
 204         return true;
 205     if ((*firstKeyA)[l] > i)
 206         l --;
 207
 208     ulong lastKey, firstKey = (*firstKeyA)[l];
 209     if (l == topCount - 1)
 210         lastKey = u;
 211     else
 212         lastKey = (*firstKeyA)[l + 1];
 213     ulong nSub = b;
 214     if (l == topCount - 1 && n%b)
 215         nSub = n%b;
 216
 217     return IsBitSet(P, (*offsetA)[l], firstKey, i, firstKey, lastKey, nSub, true);
 218 }
 219
 220 /**
 221  * Returns also rank_1(i) via the second parameter.
 222  */
 223 bool BSGAP::IsBitSet(ulong i, ulong *resultRank)
 224 {
 225     *resultRank = 0;
 226     if (n == 0) // Trivial case
 227         return false;
 228     if (i >= u)
 229     {
 230         *resultRank = n;
 231         return false; // FIXME throw std::out_of_range("BSGAP::rank(): Parameter i was out of range");
 232     }
 233     if ((*firstKeyA)[0] > i)
 234         return false;
 235     ulong l = 0, r = topCount - 1;
 236     while (l < r)   // Binary search on the array firstKeyA
 237     {
 238         ulong mid = l + (r - l)/2lu;
 239         if ((*firstKeyA)[mid] < i)
 240             l = mid + 1;
 241         else
 242             r = mid;
 243     }
 244     if ((*firstKeyA)[l] > i && l == 0)
 245         return false;
 246     if ((*firstKeyA)[l] == i)
 247     {
 248         *resultRank =  1 + l * b;
 249         return true;
 250     }
 251     if ((*firstKeyA)[l] > i)
 252         l --;
 253
 254     ulong lastKey, firstKey = (*firstKeyA)[l];
 255     if (l == topCount - 1)
 256         lastKey = u;
 257     else
 258         lastKey = (*firstKeyA)[l + 1];
 259     ulong nSub = b;
 260     if (l == topCount - 1 && n%b)
 261         nSub = n%b;
 262
 263     bool result = IsBitSet(P, (*offsetA)[l], firstKey, i, firstKey, lastKey, nSub, true, resultRank);
 264     *resultRank += l * b;
 265     return result;
 266 }
 267
 268
 269 ulong BSGAP::select(ulong i)
 270 {
 271     if (n == 0)
 272         return 0;
 273     if (i == 0)
 274         return 0;
 275     if (i > n)
 276         throw std::out_of_range("BSGAP::select(): Parameter i was out of range");
 277
 278     ulong j = i/b;
 279     if (i%b == 0)
 280         j --;
 281     ulong lastKey, firstKey = (*firstKeyA)[j];
 282     if (j == topCount - 1)
 283         lastKey = u;
 284     else
 285         lastKey = (*firstKeyA)[j + 1];
 286     ulong nSub = b;
 287     if (j == topCount - 1 && n%b)
 288         nSub = n%b;
 289
 290     return select(P, (*offsetA)[j], firstKey, i - j * b, firstKey, lastKey, nSub, true);
 291 }
 292
 293 ulong BSGAP::select0(ulong i)
 294 {
 295     if (n == 0) // Trivial case
 296         return i-1;
 297     if (i > u - n)
 298         throw std::out_of_range("BSGAP::select0(): Parameter i was out of range");
 299
 300     ulong l = 0, r = topCount - 1;
 301     while (l < r)   // Binary search on the array firstKeyA
 302     {
 303         ulong mid = l + (r - l)/2lu;
 304         if ((*firstKeyA)[mid] - mid*b < i) // Number of zeros before [mid]
 305             l = mid + 1;
 306         else
 307             r = mid;
 308     }
 309     if (l == 0 && (*firstKeyA)[l] >= i)
 310         return i-1;
 311     if (l > 0 && (*firstKeyA)[l] - l*b >= i)
 312         --l; // FIXME: TEST
 313
 314     i -= (*firstKeyA)[l] - l*b;
 315
 316     ulong lastKey, firstKey = (*firstKeyA)[l];
 317     if (l == topCount - 1)
 318         lastKey = u;
 319     else
 320         lastKey = (*firstKeyA)[l + 1];
 321     ulong nSub = b;
 322     if (l == topCount - 1 && n%b)
 323         nSub = n%b;
 324
 325     return select0(P, (*offsetA)[l], firstKey, i, firstKey, lastKey, nSub, true);
 326 }
 327
 328 ulong BSGAP::select0(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild)
 329 {
 330     while (1)
 331     {
 332         if (n == 0)
 333             return l+i;
 334
 335         // Check if subtree is full  --- FIXME Make this test an assert()!
 336         if (IsSubtreeFull(firstKey, l, r, n))
 337             throw runtime_error("BSGAP::select0(): subtree was full"); // Should not happen // (l == firstKey ? l + i - 1: l + i);
 338
 339         // Parse key value
 340         ulong key;
 341         ulong x = DeltaDecode(B, offset);
 342         if (leftChild)
 343             key = r - x;
 344         else
 345             key = l + x;
 346
 347         if (numberOfZeros(firstKey, l, key, n) < i)
 348         {
 349             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 350             i = i - numberOfZeros(firstKey, l, key, n);
 351             l = key;
 352             n = n - n/2 - 1;
 353             leftChild = false;
 354         }
 355         else
 356         {
 357             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 358             r = key;
 359             n = n/2;
 360             leftChild = true;
 361         }
 362     }
 363 }
 364
 365
 366 ulong BSGAP::rank(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild)
 367 {
 368     // Number of keys in subtrees on left-side
 369     ulong result = 0;
 370
 371     while (1)
 372     {
 373
 374         if (n == 0)
 375             return result;
 376
 377         // Check if subtree is full
 378         if (IsSubtreeFull(firstKey, l, r, n))
 379             return result + (l == firstKey ? i - l + 1 : i - l);
 380
 381         // Parse key value
 382         ulong key;
 383         ulong x = DeltaDecode(B, offset);
 384         if (leftChild)
 385             key = r - x;
 386         else
 387             key = l + x;
 388         if (key == i)
 389             return result + n/2 + 1;
 390
 391         if (key < i)
 392         {
 393             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 394             result += n/2 + 1;
 395             l = key;
 396             n = n - n/2 - 1;
 397             leftChild = false;
 398         }
 399         else
 400         {
 401             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 402             r = key;
 403             n = n/2;
 404             leftChild = true;
 405         }
 406     }
 407 }
 408
 409
 410 bool BSGAP::IsBitSet(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild)
 411 {
 412     // Number of keys in subtrees on left-side
 413     ulong result = 0;
 414
 415     while (1)
 416     {
 417
 418         if (n == 0)
 419             return false;
 420
 421         // Check if subtree is full
 422         if (IsSubtreeFull(firstKey, l, r, n))
 423             return result + (l == firstKey ? i - l + 1 : i - l);
 424
 425         // Parse key value
 426         ulong key;
 427         ulong x = DeltaDecode(B, offset);
 428         if (leftChild)
 429             key = r - x;
 430         else
 431             key = l + x;
 432         if (key == i)
 433             return true;
 434
 435         if (key < i)
 436         {
 437             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 438             result += n/2 + 1;
 439             l = key;
 440             n = n - n/2 - 1;
 441             leftChild = false;
 442         }
 443         else
 444         {
 445             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 446             r = key;
 447             n = n/2;
 448             leftChild = true;
 449         }
 450     }
 451 }
 452
 453 // Returns also the rank of i
 454 bool BSGAP::IsBitSet(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild, ulong *resultRank)
 455 {
 456     // Number of keys in subtrees on left-side
 457     ulong result = 0;
 458
 459     while (1)
 460     {
 461
 462         if (n == 0)
 463         {
 464             *resultRank = result;
 465             return false;
 466         }
 467
 468         // Check if subtree is full
 469         if (IsSubtreeFull(firstKey, l, r, n))
 470         {
 471             *resultRank = result + (l == firstKey ? i - l + 1 : i - l);
 472             return true;
 473         }
 474
 475         // Parse key value
 476         ulong key;
 477         ulong x = DeltaDecode(B, offset);
 478         if (leftChild)
 479             key = r - x;
 480         else
 481             key = l + x;
 482         if (key == i)
 483         {
 484             *resultRank = result + n/2 + 1;
 485             return true;
 486         }
 487
 488         if (key < i)
 489         {
 490             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 491             result += n/2 + 1;
 492             l = key;
 493             n = n - n/2 - 1;
 494             leftChild = false;
 495         }
 496         else
 497         {
 498             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 499             r = key;
 500             n = n/2;
 501             leftChild = true;
 502         }
 503     }
 504 }
 505
 506
 507 ulong BSGAP::select(ulong *B, ulong offset, ulong firstKey, ulong i, ulong l, ulong r, ulong n, bool leftChild)
 508 {
 509     while (1)
 510     {
 511         if (n == 0)
 512             return 0;
 513
 514         // Check if subtree is full
 515         if (IsSubtreeFull(firstKey, l, r, n))
 516             return (l == firstKey ? l + i - 1: l + i);
 517
 518         // Parse key value
 519         ulong key;
 520         ulong x = DeltaDecode(B, offset);
 521         if (leftChild)
 522             key = r - x;
 523         else
 524             key = l + x;
 525         if (n/2 + 1 == i)
 526             return key;
 527         if (n/2 + 1 < i)
 528         {
 529             offset = FindRightSubtree(B, offset, firstKey, key, l, r, n);
 530             i = i - n/2 - 1;
 531             l = key;
 532             n = n - n/2 - 1;
 533             leftChild = false;
 534         }
 535         else
 536         {
 537             offset = FindLeftSubtree(B, offset, firstKey, key, l, r, n);
 538             r =  key;
 539             n = n/2;
 540             leftChild = true;
 541         }
 542     }
 543 }
 544
 545 ulong * BSGAP::GetSubtree(ulong *B, ulong firstKey, ulong l, ulong r, ulong n, bool leftChild, ulong &bits)
 546 {
 547     bits = 0;
 548     if (n == 0)
 549         return 0;
 550
 551     // Check if subtree is full
 552     if (IsSubtreeFull(firstKey, l, r, n))
 553         return 0;
 554
 555     // Pick the median
 556     ulong key = FindMedian(B, firstKey, l, r, n);
 557
 558     ulong *leftSub, *rightSub, leftBits, rightBits;
 559     leftSub = GetSubtree(B, firstKey, l, key, n/2, true, leftBits);
 560     rightSub = GetSubtree(B, firstKey, key, r, n - n/2 - 1, false, rightBits);
 561
 562     // Encode key value
 563     ulong *keyDelta, keyBits;
 564     if (leftChild)
 565         keyDelta = DeltaEncode(r - key, keyBits, true);
 566     else
 567         keyDelta = DeltaEncode(key - l, keyBits, true);
 568
 569     // Encode jump offset if left and right subtrees exists
 570     ulong *jumpOffset = 0, jumpBits = 0;
 571     if (leftBits != 0 && rightBits != 0)
 572         jumpOffset = DeltaEncode(leftBits, jumpBits);
 573
 574     // bits is the sum of keyBits, jumpBits, leftBits and rightBits
 575     bits = keyBits + jumpBits + leftBits + rightBits;
 576     ulong *output = new ulong[bits / W + 1];
 577
 578     /**
 579      * Write "output"
 580      */
 581     BitCopy(output, keyDelta, 0, keyBits);
 582     delete [] keyDelta;
 583     ulong offset = keyBits;
 584
 585     if (jumpBits != 0)
 586     {
 587         BitCopy(output, jumpOffset, offset, jumpBits);
 588         offset += jumpBits;
 589         delete [] jumpOffset;
 590     }
 591
 592     BitCopy(output, leftSub, offset, leftBits);
 593     offset += leftBits;
 594     BitCopy(output, rightSub, offset, rightBits);
 595     offset += rightBits;
 596
 597     assert(offset == bits);
 598     if (leftBits != 0)
 599         delete [] leftSub;
 600     if (rightBits != 0)
 601         delete [] rightSub;
 602     return output;
 603 }
 604
 605 void BSGAP::BitCopy(ulong *dest, ulong *src, ulong offset, ulong len)
 606 {
 607     if (len == 0)
 608         return;
 609     ulong i;
 610
 611     if (len <= W)
 612     {
 613         Tools::SetVariableField(dest, len, offset, *src);
 614         return;
 615     }
 616
 617     for (i = 0; i < len/W; i ++)
 618         Tools::SetVariableField(dest, W, offset + i * W, src[i]);
 619
 620     if (i * W < len)
 621         Tools::SetVariableField(dest, len - i*W, offset + i * W, src[i]);
 622 }
 623
 624 ulong BSGAP::FindMedian(ulong *B, ulong firstKey, ulong l, ulong r, ulong n)
 625 {
 626     // Linear scan: slow but affects only the construction time
 627     n = n/2 + 1;
 628     if (l != firstKey)
 629         l ++;   // Skip left boundary
 630
 631     while (n && l <= r)
 632     {
 633         if (Tools::GetField(B, 1, l))
 634             n --;
 635         l ++;
 636     }
 637     return l - 1;
 638 }
 639
 640 ulong BSGAP::DeltaDecode(ulong *B, ulong &offset)
 641 {
 642     // Dense coding
 643     // http://www.dcc.uchile.cl/~gnavarro/ps/ir06.pdf
 644     const unsigned s = 180; // FIXME These should be class variables
 645     const unsigned c = 256 - s; // FIXME Choose <s,c> values!
 646     const unsigned b = 8; // one byte.
 647
 648     ulong i = 0;
 649     ulong base = 0;
 650     ulong tot = s;
 651     while (Tools::GetVariableField(B, b, offset) < c)
 652     {
 653         i = i*c + Tools::GetVariableField(B, b, offset);
 654         offset += b;
 655         base += tot;
 656         tot = tot * c;
 657     }
 658     i = i * s + (Tools::GetVariableField(B, b, offset) - c);
 659     offset += b;
 660     return i + base;
 661 }
 662
 663
 664 ulong * BSGAP::DeltaEncode(ulong value, ulong &bits, bool onlyPositive, bool negative)
 665 {
 666     // Dense coding
 667     // http://www.dcc.uchile.cl/~gnavarro/ps/ir06.pdf
 668     const unsigned s = 180;  // FIXME These should be class variables
 669     const unsigned c = 256 - s; // FIXME Choose <s,c> values!
 670     const unsigned b = 8; // one byte
 671
 672     // Calculate size first:
 673     bits = b;
 674     ulong x = value / s;
 675     while (x > 0)
 676     {
 677         --x;
 678         bits += b;
 679         x = x/c;
 680     }
 681
 682     // bits is now the length of the whole codeword
 683     ulong *B = new ulong[bits/W + 1];
 684
 685     // output codeword (backwards):
 686     unsigned offset = bits - b;
 687     ulong output = c + (value % s);
 688     Tools::SetVariableField(B, b, offset, output);
 689
 690     x = value / s;
 691     while (x > 0)
 692     {
 693         --x;
 694         offset -= b;
 695
 696         output = x % c;
 697         Tools::SetVariableField(B, b, offset, output);
 698
 699         x = x/c;
 700     }
 701
 702     return B;
 703 }
 704
 705
 706 /**
 707  * Saving the following data fields:
 708  *   ulong u, n;             // Universe size, number of 1-bits in B
 709  *   ulong topCount;   // Top structure and the number of substructures
 710  *   ulong bitsInP;
 711  *   ulong *P;               // Pointer to BSGAP structures
 712  *   unsigned b;             // Subdictionary size (\log^2 n)
 713  *   BlockArray *offsetA;    // Array of pointers (into bitvector P)
 714  *   BlockArray *firstKeyA;  // Array of first key positions of the substructures
 715  */
 716 void BSGAP::Save(FILE *file) const
 717 {
 718     if (std::fwrite(&(this->u), sizeof(ulong), 1, file) != 1)
 719         throw std::runtime_error("BSGAP::Save(): file write error (u).");
 720
 721     if (std::fwrite(&(this->n), sizeof(ulong), 1, file) != 1)
 722         throw std::runtime_error("BSGAP::Save(): file write error (n).");
 723
 724     if (n == 0)
 725         return; // Done.
 726
 727     if (std::fwrite(&(this->topCount), sizeof(ulong), 1, file) != 1)
 728         throw std::runtime_error("BSGAP::Save(): file write error (topCount).");
 729
 730     if (std::fwrite(&(this->bitsInP), sizeof(ulong), 1, file) != 1)
 731         throw std::runtime_error("BSGAP::Save(): file write error (bitsInP).");
 732
 733     for (ulong offset = 0; offset < bitsInP/W+1; offset ++)
 734     {
 735         if (std::fwrite(this->P+offset, sizeof(ulong), 1, file) != 1)
 736             throw std::runtime_error("BSGAP::Save(): file write error (P).");
 737     }
 738
 739     if (std::fwrite(&(this->b), sizeof(unsigned), 1, file) != 1)
 740         throw std::runtime_error("BSGAP::Save(): file write error (b).");
 741
 742     offsetA->Save(file);
 743     firstKeyA->Save(file);
 744 }
 745
 746 /**
 747  * Load from file
 748  */
 749 BSGAP::BSGAP(FILE *file)
 750     : u(0), n(0), topCount(0), bitsInP(0), P(0), b(0), offsetA(0), firstKeyA(0)
 751 {
 752     if (std::fread(&(this->u), sizeof(ulong), 1, file) != 1)
 753         throw std::runtime_error("BSGAP::Load(): file read error (u).");
 754
 755     if (std::fread(&(this->n), sizeof(ulong), 1, file) != 1)
 756         throw std::runtime_error("BSGAP::Load(): file read error (n).");
 757
 758     if (n == 0)
 759         return; // Done.
 760
 761     if (std::fread(&(this->topCount), sizeof(ulong), 1, file) != 1)
 762         throw std::runtime_error("BSGAP::Load(): file read error (topCount).");
 763
 764     if (std::fread(&(this->bitsInP), sizeof(ulong), 1, file) != 1)
 765         throw std::runtime_error("BSGAP::Load(): file read error (bitsInP).");
 766
 767     P = new ulong[bitsInP/W+1];
 768     for (ulong offset = 0; offset < bitsInP/W+1; offset ++)
 769     {
 770         if (std::fread(this->P+offset, sizeof(ulong), 1, file) != 1)
 771             throw std::runtime_error("BSGAP::Load(): file read error (P).");
 772     }
 773
 774     if (std::fread(&(this->b), sizeof(unsigned), 1, file) != 1)
 775         throw std::runtime_error("BSGAP::Load(): file read error (b).");
 776
 777     offsetA = new BlockArray(file);
 778     firstKeyA = new BlockArray(file);
 779 }
 780