1 /* Copyright (C) 2007,2010 LinBox 2 * Written by <Jean-Guillaume.Dumas@imag.fr> 3 * Modified by Brice Boyer (briceboyer) <boyer.brice@gmail.com> (RandomPrimeIter) 4 * Modified by Clement Pernet <clement.pernet@imag.fr> (PrimeIter) 5 * 6 * 7 * 8 * ========LICENCE======== 9 * This file is part of the library LinBox. 10 * 11 * LinBox is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU Lesser General Public 13 * License as published by the Free Software Foundation; either 14 * version 2.1 of the License, or (at your option) any later version. 15 * 16 * This library is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 19 * Lesser General Public License for more details. 20 * 21 * You should have received a copy of the GNU Lesser General Public 22 * License along with this library; if not, write to the Free Software 23 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 24 * ========LICENCE======== 25 */ 26 27 /*! @file randiter/random-prime.h 28 * @ingroup randiter 29 * @brief Generates random positive prime \ref integers. 30 */ 31 32 #ifndef __LINBOX_random_prime_iterator_H 33 #define __LINBOX_random_prime_iterator_H 34 35 36 #include <givaro/givinteger.h> 37 #include <givaro/givintprime.h> 38 #include "linbox/util/timer.h" 39 #include "linbox/util/debug.h" 40 #include "linbox/field/field-traits.h" 41 namespace LinBox 42 { 43 /*! \brief Information about the type of Prime Iterator 44 */ 45 namespace IteratorCategories { 46 //! Iterator following a deterministic sequence of primes (from the largest one, in decreasing order 47 struct DeterministicTag{}; 48 //! Iterator sampling uniformly from all primes of given bitsize 49 struct UniformTag{}; 50 //! Iterator sampling randomly (no distribution guaranteed whatsoever) from all primes of given bitsize 51 struct HeuristicTag{}; 52 }; 53 54 /*! \brief Whether a prime generator generates a sequence with non repeating 55 * numbers 56 */ 57 template<class IteratorTrait> 58 struct UniqueSamplingTrait 59 :public std::false_type { }; /* default: no guarantee of uniqueness */ 60 61 template<> 62 struct UniqueSamplingTrait<IteratorCategories::DeterministicTag> 63 :public std::true_type { }; 64 65 66 /*! @brief Prime Iterator. 67 * @ingroup primes 68 * @ingroup randiter 69 * 70 * Generates prime of specified length using a heuristically random distribution 71 * (no guarantee whatsoever). 72 * @internal 73 * It is given by <code>nextprime(2^_bits-p)</code> where <code>size(p) < _bits</code>. 74 */ 75 template<class Trait = IteratorCategories::HeuristicTag> 76 class PrimeIterator{ 77 protected: 78 uint64_t _bits; //!< common lenght of all primes 79 integer _prime; //!< the generated prime. 80 Givaro::IntPrimeDom _IPD; //!< empty struct dealing with primality. 81 82 virtual void generatePrime(); 83 public: 84 typedef integer Prime_Type ; 85 typedef UniqueSamplingTrait<Trait> UniqueSamplingTag; //!< whether a prime can be picked more than once 86 typedef Trait IteratorTag; 87 88 /*! Constructor. 89 * @param bits size of primes (in bits). Default is 23 so it 90 * can fit in a <code>Linbox::Modular<double></code>. 91 * @param seed if \c 0 a seed will be generated, otherwise, the 92 * provided seed will be use. 93 */ 94 PrimeIterator(uint64_t bits = 23, uint64_t seed = 0) : 95 _bits(bits) 96 { 97 linbox_check(bits >1); 98 if (! seed) 99 seed = BaseTimer::seed(); 100 setSeed (seed); 101 _prime = Prime_Type(1)<<bits; 102 generatePrime(); 103 } 104 105 /** @brief operator++() (prefix ++ operator) 106 * creates a new random prime. 107 */ 108 inline PrimeIterator<Trait> &operator ++ () 109 { 110 generatePrime(); 111 return *this; 112 } 113 114 /** @brief get the random prime. 115 * returns the actual prime. 116 * @warning a new prime is not generated. 117 */ 118 const Prime_Type &operator * () const {return _prime;} 119 120 /** @brief Sets the seed. 121 * Set the random seed to be \p ul. 122 * @param ul the new seed. 123 */ 124 void static setSeed(uint64_t ul){integer::seeding(ul);} 125 126 /** @brief Sets the bit size. 127 * @param bits the new bit size. 128 */ 129 void setBits(uint64_t bits) { 130 _bits = bits; 131 linbox_check(bits >1); 132 generatePrime(); 133 } 134 135 uint64_t getBits() const {return _bits;} 136 }; 137 138 template<> 139 void PrimeIterator<IteratorCategories::HeuristicTag>::generatePrime(){ 140 integer::random_exact_2exp(_prime,_bits); 141 _IPD.nextprimein(_prime); 142 while (_prime.bitsize()>_bits) 143 _IPD.prevprimein(_prime); 144 } 145 146 template<> 147 void PrimeIterator<IteratorCategories::DeterministicTag>::generatePrime(){ 148 if (_prime < 3) { 149 throw LinboxError("LinBox ERROR: Ran out of primes in Deterministic Prime Iterator.\n"); 150 } 151 _IPD.prevprimein(_prime); 152 } 153 154 template<> 155 void PrimeIterator<IteratorCategories::UniformTag>::generatePrime(){ 156 do{ 157 integer::random_exact_2exp(_prime,_bits); 158 switch (_prime %6){ 159 case 0: _prime++; break; 160 case 4: _prime++; break; 161 case 2: _prime--; break; 162 case 3: _prime+=2; break; 163 } 164 } while(!_IPD.isprime(_prime)); 165 } 166 167 168 169 /* ================ 170 * Find the log base 2 of an N-bit integer 171 * From Bit Twiddling Hacks 172 * By Sean Eron Anderson 173 * seander@cs.stanford.edu 174 * https://graphics.stanford.edu/~seander/bithacks.html 175 */ 176 static const uint32_t MultiplyDeBruijnBitPosition[32] = 177 { 178 0U, 9U, 1U, 10U, 13U, 21U, 2U, 29U, 11U, 14U, 16U, 179 18U, 22U, 25U, 3U, 30U, 8U, 12U, 20U, 28U, 15U, 17U, 180 24U, 7U, 19U, 27U, 23U, 6U, 26U, 5U, 4U, 31U 181 }; 182 183 uint32_t MultiplyDeBruijnHighestBit(uint32_t v) 184 { 185 v |= v >> 1; // first round down to one less than a power of 2 186 v |= v >> 2; 187 v |= v >> 4; 188 v |= v >> 8; 189 v |= v >> 16; 190 191 return MultiplyDeBruijnBitPosition[(uint32_t)(v * 0x07C4ACDDU) >> 27]; 192 } 193 /* ================ */ 194 195 /*! @brief Adaptor class to make a fixed-length sequence behave like a PrimeIterator. 196 */ 197 template <typename IteratorT, class UniqueTrait = std::true_type> 198 class PrimeSequence { 199 private: 200 IteratorT _cur; 201 const IteratorT _end; 202 203 public: 204 using Prime_Type = typename std::remove_cv<typename std::remove_reference<decltype(*_cur)>::type>::type; 205 using UniqueSamplingTag = UniqueTrait; 206 using IteratorTag = IteratorCategories::DeterministicTag; 207 208 PrimeSequence(IteratorT begin, IteratorT end) : 209 _cur(begin), _end(end) 210 { } 211 212 /** @brief operator++() (prefix ++ operator) 213 * moves to the next prime in the sequence. 214 */ 215 inline PrimeSequence& operator++ () { 216 if (_cur == _end) { 217 throw LinboxError("LinBox ERROR: Ran out of primes in PrimeSequence.\n"); 218 } 219 ++_cur; 220 return *this; 221 } 222 223 /** @brief get the prime. 224 * returns the actual prime. 225 * @warning a new prime is not generated. 226 */ 227 const Prime_Type &operator * () const {return *_cur;} 228 }; 229 230 /*! convenience factory to create a PrimeSequence from an STL-like container. */ 231 template <typename Container> 232 PrimeSequence<typename Container::const_iterator> 233 create_prime_sequence(const Container& cont) 234 { 235 return PrimeSequence<typename Container::const_iterator>(cont.begin(), cont.end()); 236 } 237 238 239 /*! @brief Masked Prime Iterator. 240 * @ingroup primes 241 * @ingroup randiter 242 * 243 * Generates prime of specified length with fixed lower bits 244 * @internal 245 */ 246 template<class Trait = IteratorCategories::HeuristicTag> 247 class MaskedPrimeIterator : public PrimeIterator<Trait> { 248 private: 249 const uint32_t _shift; 250 const uint32_t _fffff; 251 const uint32_t _nmask; 252 253 public: 254 typedef MaskedPrimeIterator<Trait> Self_t; 255 typedef PrimeIterator<Trait> Father_t; 256 void generatePrime(); 257 258 inline Self_t &operator ++ () 259 { 260 Father_t::operator++(); 261 return *this; 262 } 263 264 // Makes mask odd: implicit _mask=2*mask+1 265 MaskedPrimeIterator(uint32_t mask, uint32_t max, uint64_t bits = 23, uint64_t seed = 0) : 266 Father_t(bits,seed), 267 _shift( MultiplyDeBruijnHighestBit(max) + 2), 268 _fffff((1U<<_shift)-1), 269 _nmask( ~((mask<<1) | 0x1) & _fffff ) // NOT _mask 270 { 271 this->_prime |= _fffff; // set lowest bits to 11111 272 this->_prime ^= _nmask; // set lowest bits to _mask 273 generatePrime(); 274 } 275 276 const uint32_t getMask() const { return ( ~(_nmask) & _fffff ); } 277 const uint32_t getShift() const { return _shift; } 278 }; 279 280 template<> 281 void MaskedPrimeIterator<IteratorCategories::HeuristicTag>::generatePrime(){ 282 integer::random_exact_2exp(_prime,_bits); 283 284 _prime |= _fffff; // set lowest bits to 11111 285 _prime ^= _nmask; // set lowest bits to _mask 286 287 while(! _IPD.isprime(_prime) ) { 288 _prime += (1<<_shift); 289 } 290 } 291 292 template<> 293 void MaskedPrimeIterator<IteratorCategories::DeterministicTag>::generatePrime(){ 294 do { 295 _prime -= (1<<_shift); 296 if (_prime < 2) { 297 throw LinboxError("LinBox ERROR: Ran out of primes in Masked Deterministic Prime Iterator.\n"); 298 } 299 } while(! _IPD.isprime(_prime) ); 300 } 301 302 303 template<> 304 void MaskedPrimeIterator<IteratorCategories::UniformTag>::generatePrime(){ 305 do{ 306 integer::random_exact_2exp(_prime,_bits); 307 switch (_prime %6){ 308 case 0: _prime++; break; 309 case 4: _prime++; break; 310 case 2: _prime--; break; 311 case 3: _prime+=2; break; 312 } 313 314 _prime |= _fffff; // set lowest bits to 11111 315 _prime ^= _nmask; // set lowest bits to _mask 316 317 } while(!_IPD.isprime(_prime)); 318 } 319 } 320 321 #endif //__LINBOX_random_prime_iterator_H 322 323 // Local Variables: 324 // mode: C++ 325 // tab-width: 4 326 // indent-tabs-mode: nil 327 // c-basic-offset: 4 328 // End: 329 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 330