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