dcover/doubling.hpp

   1 /*
   2  * doubling.hpp
   3  *
   4  * Copyright 2003 Max-Planck-Institut fuer Informatik
   5  *
   6  */
   7
   8 //======================================================================
   9 // Doubling algorithm for suffix array construction.
  10 // Follows closely the algorithm described in
  11 //
  12 //   N.J. Larsson and K. Sadakane. Faster Suffix Sorting.
  13 //   Technical Report LU-CS-TR:99-214, Dept. of Computer Science,
  14 //   Lund University, Sweden, 1999.
  15 //======================================================================
  16
  17 #ifndef DOUBLING_HPP
  18 #define DOUBLING_HPP
  19
  20 #include "partition.hpp"
  21
  22 #include <iterator>
  23 #include <algorithm>
  24 #include <utility>
  25 #include <sstream>
  26 #include <stdexcept>
  27
  28 //----------------------------------------------------------------------
  29 // functor isa_access
  30 //
  31 // Given a suffix array iterator computes the corresponding
  32 // inverse suffix array entry to be used in comparisons.
  33 //----------------------------------------------------------------------
  34 template <typename ISAIterator, typename SAIterator>
  35 class isa_access
  36   : public std::unary_function<SAIterator,
  37            typename std::iterator_traits<ISAIterator>::value_type>
  38 {
  39 private:
  40   typedef typename std::iterator_traits<ISAIterator>::difference_type
  41           difference_type;
  42   const ISAIterator isa;
  43 public:
  44   isa_access(ISAIterator s, difference_type lcp) : isa(s+lcp) {}
  45   typename std::iterator_traits<ISAIterator>::value_type
  46   operator() (SAIterator i) const {
  47     return isa[*i];
  48   }
  49 };
  50
  51 //----------------------------------------------------------------------
  52 // function update_group
  53 //
  54 // Update a group after doubling.
  55 // WARNING: cannot handle an empty group
  56 //----------------------------------------------------------------------
  57 template <typename SAIterator, typename ISAIterator>
  58 void update_group(SAIterator beg, SAIterator end,
  59                   SAIterator sa, ISAIterator isa)
  60 {
  61   typedef typename std::iterator_traits<SAIterator>::value_type
  62     pos_type;
  63   pos_type groupnumber = std::distance(sa, end) - 1;
  64   if (end-beg > 1) {
  65     for (SAIterator i = beg; i != end; ++i) {
  66       isa[*i] = groupnumber;
  67     }
  68   } else {
  69     isa[*beg] = groupnumber;
  70     *beg = -1;
  71   }
  72 }
  73
  74
  75 //----------------------------------------------------------------------
  76 // function selection_partition
  77 //
  78 // Performs doubling partitioning for small groups.
  79 //----------------------------------------------------------------------
  80 template <typename SAIterator, typename ISAIterator>
  81 void selection_partition(SAIterator sa, ISAIterator isa,
  82                          SAIterator beg, SAIterator end,
  83          typename std::iterator_traits<SAIterator>::difference_type lcp)
  84 {
  85   typedef typename std::iterator_traits<SAIterator>::value_type
  86     pos_type;
  87   isa_access<ISAIterator, SAIterator> key(isa, lcp);
  88   while (end-beg > 1) {
  89     SAIterator min_end = beg+1;
  90     pos_type min_key = key(beg);
  91     for (SAIterator i = beg+1; i != end; ++i) {
  92       if (key(i) < min_key) {
  93         std::swap(*beg, *i);
  94         min_end = beg+1;
  95         min_key = key(beg);
  96       } else if (!(min_key < key(i))) {
  97         std::swap(*min_end, *i);
  98         ++min_end;
  99       }
 100     }
 101     update_group(beg, min_end, sa, isa);
 102     beg = min_end;
 103   }
 104   if (beg != end) update_group(beg, end, sa, isa);
 105 }
 106
 107
 108 //----------------------------------------------------------------------
 109 // function doubling_partition
 110 //
 111 // Given a group of suffixes with a common prefix of length lcp
 112 // partition it to subgroups according to the next lcp characters.
 113 // The next lcp characters are represented by a single value
 114 // in the inverse suffix array.
 115 //----------------------------------------------------------------------
 116 template <typename SAIterator, typename ISAIterator>
 117 void doubling_partition(SAIterator sa, ISAIterator isa,
 118                         SAIterator beg, SAIterator end,
 119           typename std::iterator_traits<SAIterator>::difference_type lcp)
 120 {
 121   typedef typename std::iterator_traits<SAIterator>::value_type
 122     pos_type;
 123
 124   while (end-beg > 10) {
 125
 126     // ternary partition
 127     isa_access<ISAIterator, SAIterator> get_isa_entry(isa, lcp);
 128     SAIterator pivot = random_pivot(beg, end, std::less<pos_type>(),
 129                                     get_isa_entry);
 130     std::pair<SAIterator, SAIterator> midrange
 131       = ternary_partition(beg, end, pivot, std::less<pos_type>(),
 132                           get_isa_entry);
 133
 134     // recurse on <-part
 135     if (beg != midrange.first) {
 136       doubling_partition(sa, isa, beg, midrange.first, lcp);
 137     }
 138
 139     // update =-part
 140     update_group(midrange.first, midrange.second, sa, isa);
 141
 142     // loop on >-part
 143     beg = midrange.second;
 144   }
 145
 146   selection_partition(sa, isa, beg, end, lcp);
 147 }
 148
 149
 150
 151 //======================================================================
 152 // function doubling_sort
 153 //
 154 // The main doubling algorithm
 155 //======================================================================
 156 template <typename SAIterator, typename ISAIterator>
 157 void
 158 doubling_sort(SAIterator sa, SAIterator sa_end,
 159               ISAIterator isa, ISAIterator isa_end,
 160          typename std::iterator_traits<ISAIterator>::difference_type lcp)
 161 {
 162   typedef typename std::iterator_traits<SAIterator>::difference_type
 163     difference_type;
 164
 165   difference_type size = std::distance(sa, sa_end);
 166   if (isa_end-isa != size) {
 167     std::ostringstream os;
 168     os << "doubling_sort:"
 169        << " wrong range size";
 170     throw std::logic_error(os.str());
 171   }
 172
 173   while (*sa > -size) {
 174
 175     SAIterator beg = sa;
 176     difference_type nrunlen = 0;
 177     while (beg != sa_end) {
 178       difference_type val = *beg;
 179       if (val < 0) {
 180         // jump over already sorted part
 181         beg -= val;
 182         nrunlen += val;
 183       } else {
 184         // combine preceding sorted parts
 185         if (nrunlen < 0) {
 186           *(beg+nrunlen) = nrunlen;
 187           nrunlen = 0;
 188         }
 189         // partition unsorted group
 190         SAIterator end = sa + isa[val] + 1;
 191         doubling_partition(sa, isa, beg, end, lcp);
 192         beg = end;
 193       }
 194     }
 195     if (nrunlen < 0) {
 196       *(beg+nrunlen) = nrunlen;
 197     }
 198
 199     lcp *= 2;
 200   }
 201 }
 202
 203
 204 #endif // DOUBLING_HPP