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