tools/bwt/dcs-bwt/src/difference_cover.cc

   1 // Copyright 2007 Google Inc.
   2
   3 // This program is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU General Public License
   5 // as published by the Free Software Foundation; either version 2
   6 // of the License, or (at your option) any later version.
   7
   8 // This program is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 // GNU General Public License for more details.
  12
  13 // You should have received a copy of the GNU General Public License
  14 // along with this program; if not, write to the Free Software
  15 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  16
  17 #include "difference_cover-inl.h"
  18 #include "difference_cover.h"
  19
  20 #include <cassert>     // for assert
  21 #include <iterator>    // for back_insert_iterator
  22 #include <algorithm>   // for copy, fill_n, adjacent_find, find_if
  23 #include <functional>  // for bind2nd, greater_equal, not_equal_to
  24 #include <numeric>     // for partial_sum
  25
  26 namespace dcsbwt {
  27
  28 namespace {
  29
  30 // Precomputed difference covers
  31 // They were computed using a branch-and-bound search.
  32 // All are minimal difference covers for these periods except the one
  33 // for period 256, for which the search was not completed.
  34
  35 unsigned int const cover1[] = {0};
  36 unsigned int const cover2[] = {0,1};
  37 unsigned int const cover4[] = {0,1,3};
  38 unsigned int const cover8[] = {0,1,3,4};
  39 unsigned int const cover16[] = {0,1,3,7,8};
  40 unsigned int const cover32[] = {0,1,3,7,12,17,25};
  41 unsigned int const cover64[] = {0,1,3,8,18,34,40,44,53};
  42 unsigned int const cover128[] = {0,1,3,7,17,40,55,64,75,85,104,109,117};
  43 unsigned int const cover256[] = {
  44   0,1,3,7,12,20,30,44,65,80,89,96,114,122,128,150,196,197,201,219};
  45 unsigned int const * const precomputed_covers[] = {
  46   cover1,
  47   cover2,
  48   cover4,
  49   cover8,
  50   cover16,
  51   cover32,
  52   cover64,
  53   cover128,
  54   cover256
  55 };
  56 unsigned int const precomputed_sizes[] = {
  57   1, 2, 3, 4, 5, 7, 9, 13, 20
  58 };
  59
  60
  61 // Base two logarithm of a power of two.
  62 // Fails if the input is not a power of two.
  63 unsigned int Log2(unsigned int power_of_two) {
  64   unsigned int log = 0;
  65   unsigned int p = (power_of_two >> 1);
  66   while (p) {
  67     ++log;
  68     p >>= 1;
  69   }
  70   // fail if the input was not a power of two
  71   assert((1U << log) == power_of_two);
  72   return log;
  73 }
  74
  75 // find if there is a precomputed cover
  76 bool ExistsPrecomputedCover(unsigned int period) {
  77   return Log2(period) <= Log2(256);
  78 }
  79
  80 // compute a parameter for the Colbourn-Ling algorithm (see below
  81 unsigned int ColbournLingDegree(unsigned int period) {
  82   unsigned int r = 0;
  83   while (24*r*r+36*r+13 < period) ++r;
  84   return r;
  85 }
  86
  87 } // unnamed namespace
  88
  89 //----------------------------------
  90 // DifferenceCover member functions
  91 //----------------------------------
  92
  93 // static member function for computing the cover size without construction
  94 unsigned int DifferenceCover::CoverSizeForPeriod(unsigned int period) {
  95   if (ExistsPrecomputedCover(period)) {
  96     return precomputed_sizes[Log2(period)];
  97   } else {
  98     return 4 + 6 * ColbournLingDegree(period);
  99   }
 100 }
 101
 102 // static member function for computing the object size without construction
 103 size_t DifferenceCover::SizeInBytesForPeriod(unsigned int period) {
 104   return 3 * sizeof(std::vector<unsigned int>) +
 105       (2 * static_cast<size_t>(period) +
 106        CoverSizeForPeriod(period) + 2) * sizeof(unsigned int);
 107 }
 108
 109 // constructor
 110 DifferenceCover::DifferenceCover(unsigned int period)
 111     : logperiod_(Log2(period)), mask_(period-1) {
 112   ComputeCover();
 113   ComputeRanks();
 114   ComputeCoverers();
 115   IsCorrect();
 116 }
 117
 118 // compute the difference cover
 119 // for small periods use precomputed cover
 120 // for large periods use the algorithm described in
 121 //    C.J. Colbourn and A.C.H. Ling:
 122 //    Quorums from Difference Covers.
 123 //    Information Processing Letters 75(1-2):9-12, 2000.
 124 void DifferenceCover::ComputeCover() {
 125   if (ExistsPrecomputedCover(period())) {
 126     // use precomputed cover
 127     unsigned int size = precomputed_sizes[logperiod_];
 128     cover_.reserve(size);
 129     unsigned int const * cp = precomputed_covers[logperiod_];
 130     std::back_insert_iterator<std::vector<unsigned int> > it(cover_);
 131     std::copy(cp, cp+size, it);
 132   } else {
 133     // use the Colbourn-Ling algorithm
 134     unsigned int r = ColbournLingDegree(period());
 135     unsigned int size = 4 + 6*r;
 136     cover_.reserve(size);
 137     std::back_insert_iterator<std::vector<unsigned int> > it(cover_);
 138     *it++ = 0;
 139     std::fill_n(it, r, 1);
 140     *it++ = r+1;
 141     std::fill_n(it, r, 2*r+1);
 142     std::fill_n(it, 2*r+1, 4*r+3);
 143     std::fill_n(it, r+1, 2*r+2);
 144     std::fill_n(it, r, 1);
 145     std::partial_sum(cover_.begin(), cover_.end(), cover_.begin());
 146
 147     // In general, at this point the largest (last) elements could
 148     // be larger than the period length, and then one would have to
 149     // take modulo period, re-sort the elements and remove duplicates.
 150     // However, if period > 90 then all elements are in the period.
 151     // We start using Colbourn-Ling from period length 256, so just
 152     // do a check here to be safe.
 153     assert(cover_.back() < period());
 154   }
 155 }
 156
 157 // Compute the lookup table for finding ranks.
 158 // (See IsCorrect below.)
 159 void DifferenceCover::ComputeRanks() {
 160   ranks_.reserve(period());
 161   std::back_insert_iterator<std::vector<unsigned int> > inserter(ranks_);
 162   std::vector<unsigned int>::const_iterator it;
 163   *inserter++ = 0;
 164   unsigned int pos = 0;
 165   for (it = cover_.begin(); it != cover_.end(); ++it) {
 166     std::fill_n(inserter, (*it)-pos, it-cover_.begin());
 167     pos = *it;
 168   }
 169   std::fill_n(inserter, period()-pos-1, size());
 170   assert(ranks_.size() == period());
 171 }
 172
 173 // Compute lookup table for finding cover elements for a given difference.
 174 // (See IsCorrect below.)
 175 void DifferenceCover::ComputeCoverers() {
 176   coverers_.resize(period());
 177   coverers_[0] = cover_.front();
 178   std::vector<unsigned int>::const_iterator it1, it2;
 179   for (it1 = cover_.begin(); it1 != cover_.end(); ++it1) {
 180     for (it2 = cover_.begin(); it2 != it1; ++it2) {
 181       assert(*it1 > *it2);
 182       unsigned int diff = (*it1 - *it2);
 183       coverers_[diff] = *it2;
 184       coverers_[period()-diff] = *it1;
 185     }
 186   }
 187 }
 188
 189 // Verify the correctness of the computed tables.
 190 bool DifferenceCover::IsCorrect() const {
 191
 192   // cover_ should contain distinct integers from [0,period)
 193   // in increasing order.
 194   // Check increasing order.
 195   assert(cover_.end() ==
 196            std::adjacent_find(cover_.begin(), cover_.end(),
 197                               std::greater_equal<unsigned int>() ));
 198   // Check that they are in [0,period)
 199   assert(cover_.back() < period());
 200
 201   // ranks_ should contain 0s up to and including the first cover element,
 202   // 1s up to and including the second element, etc.
 203   // For example, if period==16:
 204   //   cover_ 0 1   3       7 8
 205   //   ranks_ 0 1 2 2 3 3 3 3 4 5 5 5 5 5 5 5
 206   assert(ranks_.size() == period());
 207   std::vector<unsigned int>::const_iterator it = ranks_.begin();
 208   for (unsigned int i = 0; i < size(); ++i) {
 209     it = std::find_if(it, ranks_.end(),
 210                       std::bind2nd(std::not_equal_to<unsigned int>(), i) );
 211     assert(static_cast<unsigned int>(it-ranks_.begin()) == cover_[i]+1);
 212   }
 213   it = std::find_if(it, ranks_.end(),
 214                     std::bind2nd(std::not_equal_to<unsigned int>(), size()) );
 215   assert(ranks_.end() == it);
 216
 217   // coverers_ should contain at each position diff in [0,period)
 218   // an integer i such that both i and (i+diff) mod period are in
 219   // the cover.
 220   // For example, if period==16 and _cover=={0,1,3,7,8}
 221   //                 diff  0 1 2 3 4 5 6 7 8  9 10 11 12 13 14 15
 222   //      _coverers[diff]  0 0 1 0 3 3 1 0 0  7  7  8  7  3  3  1
 223   // _coverers[diff]+diff  0 1 3 3 7 8 7 7 8 16 17 19 19 16 17 16
 224   // (_coverers[diff]+diff) mod period        0  1  3  3  0  1  0
 225   assert(coverers_.size() == period());
 226   for (unsigned int diff = 0; diff < period(); ++diff) {
 227     unsigned int i = coverers_[diff];
 228     assert(i < period());
 229     unsigned int j = ModuloPeriod(i + diff);
 230     assert(Contains(i));
 231     assert(Contains(j));
 232   }
 233
 234   return true;
 235 }
 236
 237 }  // namespace dcsbwt