1 /* _gmp_rand (rp, state, nbits) -- Generate a random bitstream of
2 length NBITS in RP. RP must have enough space allocated to hold
5 Copyright (C) 1999, 2000 Free Software Foundation, Inc.
7 This file is part of the GNU MP Library.
9 The GNU MP Library is free software; you can redistribute it and/or modify
10 it under the terms of the GNU Lesser General Public License as published by
11 the Free Software Foundation; either version 2.1 of the License, or (at your
12 option) any later version.
14 The GNU MP Library is distributed in the hope that it will be useful, but
15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17 License for more details.
19 You should have received a copy of the GNU Lesser General Public License
20 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
21 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
22 MA 02111-1307, USA. */
28 /* For linear congruential (LC), we use one of algorithms (1) or (2).
29 (gmp-3.0 uses algorithm (1) with 'm' as a power of 2.)
35 [D. Knuth, "The Art of Computer Programming: Volume 2, Seminumerical Algorithms",
36 Third Edition, Addison Wesley, 1998, pp. 184-185.]
38 X is the seed and the result
40 a mod 8 = 5 [3.2.1.2] and [3.2.1.3]
42 its binary or decimal digits is not a simple, regular pattern
43 it has no large quotients when Euclid's algorithm is used to find
45 it passes the spectral test [3.3.4]
46 it passes several tests of [3.3.2]
47 c has no factor in common with m (c=1 or c=a can be good)
49 is a power of 2 [3.2.1.1]
51 The least significant digits of the generated number are not very
52 random. It should be regarded as a random fraction X/m. To get a
53 random integer between 0 and n-1, multiply X/m by n and truncate.
54 (Don't use X/n [ex 3.4.1-3])
56 The ``accuracy'' in t dimensions is one part in ``the t'th root of m'' [3.3.4].
58 Don't generate more than about m/1000 numbers without changing a, c, or m.
60 The sequence length depends on chosen a,c,m.
65 X = a * (X mod q) - r * (long) (X/q)
70 X is the seed and the result
71 as a seed is nonzero and less than m
72 a is a primitive root of m (which means that a^2 <= m)
75 m is a prime number near the largest easily computed integer
79 X = a * (X % ((long) m / a)) -
80 (M % a) * ((long) (X / ((long) m / a)))
82 Since m is prime, the least-significant bits of X are just as random as
83 the most-significant bits. */
85 /* Blum, Blum, and Shub.
87 [Bruce Schneier, "Applied Cryptography", Second Edition, John Wiley
88 & Sons, Inc., 1996, pp. 417-418.]
90 "Find two large prime numbers, p and q, which are congruent to 3
91 modulo 4. The product of those numbers, n, is a blum integer.
92 Choose another random integer, x, which is relatively prime to n.
95 That's the seed for the generator."
97 To generate a random bit, compute
99 The least significant bit of x[i] is the one we want.
101 We can use more than one bit from x[i], namely the
102 log2(bitlength of x[i])
103 least significant bits of x[i].
105 So, for a 32-bit seed we get 5 bits per computation.
107 The non-predictability of this generator is based on the difficulty
111 /* -------------------------------------------------- */
113 /* lc (rp, state) -- Generate next number in LC sequence. Return the
114 number of valid bits in the result. NOTE: If 'm' is a power of 2
115 (m2exp != 0), discard the lower half of the result. */
120 lc (mp_ptr rp, gmp_randstate_t rstate)
124 gmp_randstate_t rstate;
127 mp_ptr tp, seedp, ap;
129 mp_size_t tn, seedn, an;
132 unsigned long int m2exp;
136 m2exp = rstate->algdata.lc->m2exp;
137 c = (mp_limb_t) rstate->algdata.lc->c;
139 seedp = PTR (rstate->seed);
140 seedn = SIZ (rstate->seed);
144 /* Seed is 0. Result is C % M. */
149 /* M is a power of 2. */
150 if (m2exp < BITS_PER_MP_LIMB)
152 /* Only necessary when M may be smaller than C. */
153 *rp &= (((mp_limb_t) 1 << m2exp) - 1);
158 /* M is not a power of 2. */
159 abort (); /* FIXME. */
162 /* Save result as next seed. */
164 SIZ (rstate->seed) = 1;
165 return BITS_PER_MP_LIMB;
168 ap = PTR (rstate->algdata.lc->a);
169 an = SIZ (rstate->algdata.lc->a);
171 /* Allocate temporary storage. Let there be room for calculation of
172 (A * seed + C) % M, or M if bigger than that. */
174 ASSERT_ALWAYS (m2exp != 0); /* FIXME. */
178 tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
183 mpn_mul_basecase (tp, seedp, seedn, ap, an);
185 mpn_mul_basecase (tp, ap, an, seedp, seedn);
194 /* M is a power of 2. The mod operation is trivial. */
196 tp[m2exp / BITS_PER_MP_LIMB] &= ((mp_limb_t) 1 << m2exp % BITS_PER_MP_LIMB) - 1;
197 tn = (m2exp + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
201 abort (); /* FIXME. */
204 /* Save result as next seed. */
205 MPN_COPY (PTR (rstate->seed), tp, tn);
206 SIZ (rstate->seed) = tn;
210 /* Discard the lower half of the result. */
211 unsigned long int discardb = m2exp / 2;
212 mp_size_t discardl = discardb / BITS_PER_MP_LIMB;
217 if (discardb % BITS_PER_MP_LIMB != 0)
219 mpn_rshift (tp, tp + discardl, tn, discardb % BITS_PER_MP_LIMB);
220 MPN_COPY (rp, tp, (discardb + BITS_PER_MP_LIMB -1) / BITS_PER_MP_LIMB);
222 else /* Even limb boundary. */
223 MPN_COPY_INCR (rp, tp + discardl, tn);
228 MPN_COPY (rp, tp, tn);
233 /* Return number of valid bits in the result. */
235 retval = (m2exp + 1) / 2;
237 retval = SIZ (rstate->algdata.lc->m) * BITS_PER_MP_LIMB - shiftcount;
242 /* Set even bits to EVENBITS and odd bits to ! EVENBITS in RP.
243 Number of bits is m2exp in state. */
246 lc_test (mp_ptr rp, gmp_randstate_t s, const int evenbits)
248 unsigned long int rn, nbits;
251 nbits = s->algdata.lc->m2exp / 2;
252 rn = nbits / BITS_PER_MP_LIMB + (nbits % BITS_PER_MP_LIMB != 0);
255 for (f = 0; f < nbits; f++)
257 mpn_lshift (rp, rp, rn, 1);
258 if (f % 2 == ! evenbits)
264 #endif /* RAWRANDEBUG */
268 _gmp_rand (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
270 _gmp_rand (rp, rstate, nbits)
272 gmp_randstate_t rstate;
273 unsigned long int nbits;
276 mp_size_t rn; /* Size of R. */
278 rn = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
282 case GMP_RAND_ALG_LC:
284 unsigned long int rbitpos;
292 chunk_nbits = rstate->algdata.lc->m2exp / 2;
293 tn = (chunk_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
295 tp = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB);
298 while (rbitpos + chunk_nbits <= nbits)
300 mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
302 if (rbitpos % BITS_PER_MP_LIMB != 0)
304 mp_limb_t savelimb, rcy;
305 /* Target of of new chunk is not bit aligned. Use temp space
306 and align things by shifting it up. */
309 rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
311 /* bogus */ if ((chunk_nbits % BITS_PER_MP_LIMB + rbitpos % BITS_PER_MP_LIMB)
317 /* Target of of new chunk is bit aligned. Let `lc' put bits
318 directly into our target variable. */
321 rbitpos += chunk_nbits;
324 /* Handle last [0..chunk_nbits) bits. */
325 if (rbitpos != nbits)
327 mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
328 int last_nbits = nbits - rbitpos;
329 tn = (last_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
331 if (rbitpos % BITS_PER_MP_LIMB != 0)
333 mp_limb_t savelimb, rcy;
334 /* Target of of new chunk is not bit aligned. Use temp space
335 and align things by shifting it up. */
337 rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
339 if (rbitpos + tn * BITS_PER_MP_LIMB - rbitpos % BITS_PER_MP_LIMB < nbits)
344 MPN_COPY (r2p, tp, tn);
346 /* Mask off top bits if needed. */
347 if (nbits % BITS_PER_MP_LIMB != 0)
348 rp[nbits / BITS_PER_MP_LIMB]
349 &= ~ ((~(mp_limb_t) 0) << nbits % BITS_PER_MP_LIMB);
357 gmp_errno |= GMP_ERROR_UNSUPPORTED_ARGUMENT;