Annotation of OpenXM_contrib/gmp/randraw.c, Revision 1.1.1.1
1.1 maekawa 1: /* _gmp_rand (rp, state, nbits) -- Generate a random bitstream of
2: length NBITS in RP. RP must have enough space allocated to hold
3: NBITS.
4:
5: Copyright (C) 1999, 2000 Free Software Foundation, Inc.
6:
7: This file is part of the GNU MP Library.
8:
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.
13:
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.
18:
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. */
23:
24: #include "gmp.h"
25: #include "gmp-impl.h"
26: #include "longlong.h"
27:
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.)
30:
31: LC algorithm (1).
32:
33: X = (aX + c) mod m
34:
35: [D. Knuth, "The Art of Computer Programming: Volume 2, Seminumerical Algorithms",
36: Third Edition, Addison Wesley, 1998, pp. 184-185.]
37:
38: X is the seed and the result
39: a is chosen so that
40: a mod 8 = 5 [3.2.1.2] and [3.2.1.3]
41: .01m < a < .99m
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
44: gcd(a, m) [3.3.3]
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)
48: m is large (2^30)
49: is a power of 2 [3.2.1.1]
50:
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])
55:
56: The ``accuracy'' in t dimensions is one part in ``the t'th root of m'' [3.3.4].
57:
58: Don't generate more than about m/1000 numbers without changing a, c, or m.
59:
60: The sequence length depends on chosen a,c,m.
61:
62:
63: LC algorithm (2).
64:
65: X = a * (X mod q) - r * (long) (X/q)
66: if X<0 then X+=m
67:
68: [Knuth, pp. 185-186.]
69:
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)
73: q is (long) m / a
74: r is m mod a
75: m is a prime number near the largest easily computed integer
76:
77: which gives
78:
79: X = a * (X % ((long) m / a)) -
80: (M % a) * ((long) (X / ((long) m / a)))
81:
82: Since m is prime, the least-significant bits of X are just as random as
83: the most-significant bits. */
84:
85: /* Blum, Blum, and Shub.
86:
87: [Bruce Schneier, "Applied Cryptography", Second Edition, John Wiley
88: & Sons, Inc., 1996, pp. 417-418.]
89:
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.
93: Compute
94: x[0] = x^2 mod n
95: That's the seed for the generator."
96:
97: To generate a random bit, compute
98: x[i] = x[i-1]^2 mod n
99: The least significant bit of x[i] is the one we want.
100:
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].
104:
105: So, for a 32-bit seed we get 5 bits per computation.
106:
107: The non-predictability of this generator is based on the difficulty
108: of factoring n.
109: */
110:
111: /* -------------------------------------------------- */
112:
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. */
116:
117: static
118: unsigned long int
119: #if __STDC__
120: lc (mp_ptr rp, gmp_randstate_t rstate)
121: #else
122: lc (rp, rstate)
123: mp_ptr rp;
124: gmp_randstate_t rstate;
125: #endif
126: {
127: mp_ptr tp, seedp, ap;
128: mp_size_t ta;
129: mp_size_t tn, seedn, an;
130: mp_size_t retval;
131: int shiftcount = 0;
132: unsigned long int m2exp;
133: mp_limb_t c;
134: TMP_DECL (mark);
135:
136: m2exp = rstate->algdata.lc->m2exp;
137: c = (mp_limb_t) rstate->algdata.lc->c;
138:
139: seedp = PTR (rstate->seed);
140: seedn = SIZ (rstate->seed);
141:
142: if (seedn == 0)
143: {
144: /* Seed is 0. Result is C % M. */
145: *rp = c;
146:
147: if (m2exp != 0)
148: {
149: /* M is a power of 2. */
150: if (m2exp < BITS_PER_MP_LIMB)
151: {
152: /* Only necessary when M may be smaller than C. */
153: *rp &= (((mp_limb_t) 1 << m2exp) - 1);
154: }
155: }
156: else
157: {
158: /* M is not a power of 2. */
159: abort (); /* FIXME. */
160: }
161:
162: /* Save result as next seed. */
163: *seedp = *rp;
164: SIZ (rstate->seed) = 1;
165: return BITS_PER_MP_LIMB;
166: }
167:
168: ap = PTR (rstate->algdata.lc->a);
169: an = SIZ (rstate->algdata.lc->a);
170:
171: /* Allocate temporary storage. Let there be room for calculation of
172: (A * seed + C) % M, or M if bigger than that. */
173:
174: ASSERT_ALWAYS (m2exp != 0); /* FIXME. */
175:
176: TMP_MARK (mark);
177: ta = an + seedn + 1;
178: tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
179: MPN_ZERO (tp, ta);
180:
181: /* t = a * seed */
182: if (seedn >= an)
183: mpn_mul_basecase (tp, seedp, seedn, ap, an);
184: else
185: mpn_mul_basecase (tp, ap, an, seedp, seedn);
186: tn = an + seedn;
187:
188: /* t = t + c */
189: mpn_incr_u (tp, c);
190:
191: /* t = t % m */
192: if (m2exp != 0)
193: {
194: /* M is a power of 2. The mod operation is trivial. */
195:
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;
198: }
199: else
200: {
201: abort (); /* FIXME. */
202: }
203:
204: /* Save result as next seed. */
205: MPN_COPY (PTR (rstate->seed), tp, tn);
206: SIZ (rstate->seed) = tn;
207:
208: if (m2exp != 0)
209: {
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;
213:
214: tn -= discardl;
215: if (tn > 0)
216: {
217: if (discardb % BITS_PER_MP_LIMB != 0)
218: {
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);
221: }
222: else /* Even limb boundary. */
223: MPN_COPY_INCR (rp, tp + discardl, tn);
224: }
225: }
226: else
227: {
228: MPN_COPY (rp, tp, tn);
229: }
230:
231: TMP_FREE (mark);
232:
233: /* Return number of valid bits in the result. */
234: if (m2exp != 0)
235: retval = (m2exp + 1) / 2;
236: else
237: retval = SIZ (rstate->algdata.lc->m) * BITS_PER_MP_LIMB - shiftcount;
238: return retval;
239: }
240:
241: #ifdef RAWRANDEBUG
242: /* Set even bits to EVENBITS and odd bits to ! EVENBITS in RP.
243: Number of bits is m2exp in state. */
244: /* FIXME: Remove. */
245: unsigned long int
246: lc_test (mp_ptr rp, gmp_randstate_t s, const int evenbits)
247: {
248: unsigned long int rn, nbits;
249: int f;
250:
251: nbits = s->algdata.lc->m2exp / 2;
252: rn = nbits / BITS_PER_MP_LIMB + (nbits % BITS_PER_MP_LIMB != 0);
253: MPN_ZERO (rp, rn);
254:
255: for (f = 0; f < nbits; f++)
256: {
257: mpn_lshift (rp, rp, rn, 1);
258: if (f % 2 == ! evenbits)
259: rp[0] += 1;
260: }
261:
262: return nbits;
263: }
264: #endif /* RAWRANDEBUG */
265:
266: void
267: #if __STDC__
268: _gmp_rand (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
269: #else
270: _gmp_rand (rp, rstate, nbits)
271: mp_ptr rp;
272: gmp_randstate_t rstate;
273: unsigned long int nbits;
274: #endif
275: {
276: mp_size_t rn; /* Size of R. */
277:
278: rn = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
279:
280: switch (rstate->alg)
281: {
282: case GMP_RAND_ALG_LC:
283: {
284: unsigned long int rbitpos;
285: int chunk_nbits;
286: mp_ptr tp;
287: mp_size_t tn;
288: TMP_DECL (lcmark);
289:
290: TMP_MARK (lcmark);
291:
292: chunk_nbits = rstate->algdata.lc->m2exp / 2;
293: tn = (chunk_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
294:
295: tp = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB);
296:
297: rbitpos = 0;
298: while (rbitpos + chunk_nbits <= nbits)
299: {
300: mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
301:
302: if (rbitpos % BITS_PER_MP_LIMB != 0)
303: {
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. */
307: lc (tp, rstate);
308: savelimb = r2p[0];
309: rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
310: r2p[0] |= savelimb;
311: /* bogus */ if ((chunk_nbits % BITS_PER_MP_LIMB + rbitpos % BITS_PER_MP_LIMB)
312: > BITS_PER_MP_LIMB)
313: r2p[tn] = rcy;
314: }
315: else
316: {
317: /* Target of of new chunk is bit aligned. Let `lc' put bits
318: directly into our target variable. */
319: lc (r2p, rstate);
320: }
321: rbitpos += chunk_nbits;
322: }
323:
324: /* Handle last [0..chunk_nbits) bits. */
325: if (rbitpos != nbits)
326: {
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;
330: lc (tp, rstate);
331: if (rbitpos % BITS_PER_MP_LIMB != 0)
332: {
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. */
336: savelimb = r2p[0];
337: rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
338: r2p[0] |= savelimb;
339: if (rbitpos + tn * BITS_PER_MP_LIMB - rbitpos % BITS_PER_MP_LIMB < nbits)
340: r2p[tn] = rcy;
341: }
342: else
343: {
344: MPN_COPY (r2p, tp, tn);
345: }
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);
350: }
351:
352: TMP_FREE (lcmark);
353: break;
354: }
355:
356: default:
357: gmp_errno |= GMP_ERROR_UNSUPPORTED_ARGUMENT;
358: break;
359: }
360: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to randraw.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp