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 1999, 2000, 2001 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. */
119
lc (mp_ptr rp, gmp_randstate_t rstate)
121
mp_ptr tp, seedp, ap;
123
mp_size_t tn, seedn, an;
126
unsigned long int m2exp;
130
m2exp = rstate->_mp_algdata._mp_lc->_mp_m2exp;
131
c = (mp_limb_t) rstate->_mp_algdata._mp_lc->_mp_c;
133
seedp = PTR (rstate->_mp_seed);
134
seedn = SIZ (rstate->_mp_seed);
138
/* Seed is 0. Result is C % M. */
143
/* M is a power of 2. */
144
if (m2exp < BITS_PER_MP_LIMB)
146
/* Only necessary when M may be smaller than C. */
147
*rp &= (((mp_limb_t) 1 << m2exp) - 1);
152
/* M is not a power of 2. */
153
ASSERT_ALWAYS (0); /* FIXME. */
156
/* Save result as next seed. */
158
SIZ (rstate->_mp_seed) = 1;
159
return BITS_PER_MP_LIMB;
162
ap = PTR (rstate->_mp_algdata._mp_lc->_mp_a);
163
an = SIZ (rstate->_mp_algdata._mp_lc->_mp_a);
165
/* Allocate temporary storage. Let there be room for calculation of
166
(A * seed + C) % M, or M if bigger than that. */
168
ASSERT_ALWAYS (m2exp != 0); /* FIXME. */
172
tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
177
mpn_mul (tp, seedp, seedn, ap, an);
179
mpn_mul (tp, ap, an, seedp, seedn);
183
MPN_INCR_U (tp, tn, c);
188
/* M is a power of 2. The mod operation is trivial. */
190
tp[m2exp / BITS_PER_MP_LIMB] &= ((mp_limb_t) 1 << m2exp % BITS_PER_MP_LIMB) - 1;
191
tn = (m2exp + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
195
ASSERT_ALWAYS (0); /* FIXME. */
198
/* Save result as next seed. */
199
MPN_COPY (PTR (rstate->_mp_seed), tp, tn);
200
SIZ (rstate->_mp_seed) = tn;
204
/* Discard the lower half of the result. */
205
unsigned long int discardb = m2exp / 2;
206
mp_size_t discardl = discardb / BITS_PER_MP_LIMB;
211
if (discardb % BITS_PER_MP_LIMB != 0)
213
mpn_rshift (tp, tp + discardl, tn, discardb % BITS_PER_MP_LIMB);
214
MPN_COPY (rp, tp, (discardb + BITS_PER_MP_LIMB -1) / BITS_PER_MP_LIMB);
216
else /* Even limb boundary. */
217
MPN_COPY_INCR (rp, tp + discardl, tn);
222
MPN_COPY (rp, tp, tn);
227
/* Return number of valid bits in the result. */
229
retval = (m2exp + 1) / 2;
231
retval = SIZ (rstate->_mp_algdata._mp_lc->_mp_m) * BITS_PER_MP_LIMB - shiftcount;
236
/* Set even bits to EVENBITS and odd bits to ! EVENBITS in RP.
237
Number of bits is m2exp in state. */
240
lc_test (mp_ptr rp, gmp_randstate_t s, const int evenbits)
242
unsigned long int rn, nbits;
245
nbits = s->_mp_algdata._mp_lc->_mp_m2exp / 2;
246
rn = nbits / BITS_PER_MP_LIMB + (nbits % BITS_PER_MP_LIMB != 0);
249
for (f = 0; f < nbits; f++)
251
mpn_lshift (rp, rp, rn, 1);
252
if (f % 2 == ! evenbits)
258
#endif /* RAWRANDEBUG */
261
_gmp_rand (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
263
mp_size_t rn; /* Size of R. */
265
rn = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
267
switch (rstate->_mp_alg)
269
case GMP_RAND_ALG_LC:
271
unsigned long int rbitpos;
279
chunk_nbits = rstate->_mp_algdata._mp_lc->_mp_m2exp / 2;
280
tn = (chunk_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
282
tp = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB);
285
while (rbitpos + chunk_nbits <= nbits)
287
mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
289
if (rbitpos % BITS_PER_MP_LIMB != 0)
291
mp_limb_t savelimb, rcy;
292
/* Target of of new chunk is not bit aligned. Use temp space
293
and align things by shifting it up. */
296
rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
298
/* bogus */ if ((chunk_nbits % BITS_PER_MP_LIMB + rbitpos % BITS_PER_MP_LIMB)
304
/* Target of of new chunk is bit aligned. Let `lc' put bits
305
directly into our target variable. */
308
rbitpos += chunk_nbits;
311
/* Handle last [0..chunk_nbits) bits. */
312
if (rbitpos != nbits)
314
mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
315
int last_nbits = nbits - rbitpos;
316
tn = (last_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
318
if (rbitpos % BITS_PER_MP_LIMB != 0)
320
mp_limb_t savelimb, rcy;
321
/* Target of of new chunk is not bit aligned. Use temp space
322
and align things by shifting it up. */
324
rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
326
if (rbitpos + tn * BITS_PER_MP_LIMB - rbitpos % BITS_PER_MP_LIMB < nbits)
331
MPN_COPY (r2p, tp, tn);
333
/* Mask off top bits if needed. */
334
if (nbits % BITS_PER_MP_LIMB != 0)
335
rp[nbits / BITS_PER_MP_LIMB]
336
&= ~ ((~(mp_limb_t) 0) << nbits % BITS_PER_MP_LIMB);
added
removed
gmp3/randraw.c
