Annotation of OpenXM_contrib/gmp/randraw.c, Revision 1.1.1.2
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:
1.1.1.2 ! ohara 5: Copyright 1999, 2000, 2001, 2002 Free Software Foundation, Inc.
1.1 maekawa 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:
86: /* lc (rp, state) -- Generate next number in LC sequence. Return the
87: number of valid bits in the result. NOTE: If 'm' is a power of 2
88: (m2exp != 0), discard the lower half of the result. */
89:
90: static
91: unsigned long int
92: lc (mp_ptr rp, gmp_randstate_t rstate)
93: {
94: mp_ptr tp, seedp, ap;
95: mp_size_t ta;
96: mp_size_t tn, seedn, an;
97: unsigned long int m2exp;
98: mp_limb_t c;
99: TMP_DECL (mark);
100:
1.1.1.2 ! ohara 101: m2exp = rstate->_mp_algdata._mp_lc->_mp_m2exp;
1.1 maekawa 102:
1.1.1.2 ! ohara 103: /* The code below assumes the mod part is a power of two. Make sure
! 104: that is the case. */
! 105: ASSERT_ALWAYS (m2exp != 0);
! 106:
! 107: c = (mp_limb_t) rstate->_mp_algdata._mp_lc->_mp_c;
! 108:
! 109: seedp = PTR (rstate->_mp_seed);
! 110: seedn = SIZ (rstate->_mp_seed);
1.1 maekawa 111:
112: if (seedn == 0)
113: {
1.1.1.2 ! ohara 114: /* Seed is 0. Result is C % M. Assume table is sensibly stored,
! 115: with C smaller than M*/
1.1 maekawa 116: *rp = c;
117:
1.1.1.2 ! ohara 118: *seedp = c;
! 119: SIZ (rstate->_mp_seed) = 1;
! 120: return m2exp;
1.1 maekawa 121: }
122:
1.1.1.2 ! ohara 123: ap = PTR (rstate->_mp_algdata._mp_lc->_mp_a);
! 124: an = SIZ (rstate->_mp_algdata._mp_lc->_mp_a);
1.1 maekawa 125:
126: /* Allocate temporary storage. Let there be room for calculation of
127: (A * seed + C) % M, or M if bigger than that. */
128:
129: TMP_MARK (mark);
130: ta = an + seedn + 1;
131: tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
132:
133: /* t = a * seed */
134: if (seedn >= an)
1.1.1.2 ! ohara 135: mpn_mul (tp, seedp, seedn, ap, an);
1.1 maekawa 136: else
1.1.1.2 ! ohara 137: mpn_mul (tp, ap, an, seedp, seedn);
1.1 maekawa 138: tn = an + seedn;
139:
140: /* t = t + c */
1.1.1.2 ! ohara 141: tp[tn] = 0; /* sentinel, stops MPN_INCR_U */
! 142: MPN_INCR_U (tp, tn, c);
1.1 maekawa 143:
1.1.1.2 ! ohara 144: ASSERT_ALWAYS (m2exp / GMP_NUMB_BITS < ta);
1.1 maekawa 145:
1.1.1.2 ! ohara 146: /* t = t % m */
! 147: tp[m2exp / GMP_NUMB_BITS] &= ((mp_limb_t) 1 << m2exp % GMP_NUMB_BITS) - 1;
! 148: tn = (m2exp + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
1.1 maekawa 149:
150: /* Save result as next seed. */
1.1.1.2 ! ohara 151: MPN_COPY (PTR (rstate->_mp_seed), tp, tn);
! 152: SIZ (rstate->_mp_seed) = tn;
1.1 maekawa 153:
1.1.1.2 ! ohara 154: {
! 155: /* Discard the lower m2exp/2 bits of result. */
! 156: unsigned long int bits = m2exp / 2;
! 157: mp_size_t xn = bits / GMP_NUMB_BITS;
! 158:
! 159: tn -= xn;
! 160: if (tn > 0)
! 161: {
! 162: unsigned int cnt = bits % GMP_NUMB_BITS;
! 163: if (cnt != 0)
! 164: {
! 165: mpn_rshift (tp, tp + xn, tn, cnt);
! 166: MPN_COPY_INCR (rp, tp, xn + 1);
! 167: }
! 168: else /* Even limb boundary. */
! 169: MPN_COPY_INCR (rp, tp + xn, tn);
! 170: }
! 171: }
1.1 maekawa 172:
173: TMP_FREE (mark);
174:
175: /* Return number of valid bits in the result. */
1.1.1.2 ! ohara 176: return (m2exp + 1) / 2;
1.1 maekawa 177: }
178:
179: #ifdef RAWRANDEBUG
180: /* Set even bits to EVENBITS and odd bits to ! EVENBITS in RP.
181: Number of bits is m2exp in state. */
182: /* FIXME: Remove. */
183: unsigned long int
184: lc_test (mp_ptr rp, gmp_randstate_t s, const int evenbits)
185: {
186: unsigned long int rn, nbits;
187: int f;
188:
1.1.1.2 ! ohara 189: nbits = s->_mp_algdata._mp_lc->_mp_m2exp / 2;
! 190: rn = nbits / GMP_NUMB_BITS + (nbits % GMP_NUMB_BITS != 0);
1.1 maekawa 191: MPN_ZERO (rp, rn);
192:
193: for (f = 0; f < nbits; f++)
194: {
195: mpn_lshift (rp, rp, rn, 1);
196: if (f % 2 == ! evenbits)
197: rp[0] += 1;
198: }
199:
200: return nbits;
201: }
202: #endif /* RAWRANDEBUG */
203:
204: void
205: _gmp_rand (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
206: {
207: mp_size_t rn; /* Size of R. */
208:
1.1.1.2 ! ohara 209: rn = (nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
1.1 maekawa 210:
1.1.1.2 ! ohara 211: switch (rstate->_mp_alg)
1.1 maekawa 212: {
213: case GMP_RAND_ALG_LC:
214: {
215: unsigned long int rbitpos;
216: int chunk_nbits;
217: mp_ptr tp;
218: mp_size_t tn;
219: TMP_DECL (lcmark);
220:
221: TMP_MARK (lcmark);
222:
1.1.1.2 ! ohara 223: chunk_nbits = rstate->_mp_algdata._mp_lc->_mp_m2exp / 2;
! 224: tn = (chunk_nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
1.1 maekawa 225:
226: tp = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB);
227:
228: rbitpos = 0;
229: while (rbitpos + chunk_nbits <= nbits)
230: {
1.1.1.2 ! ohara 231: mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS;
1.1 maekawa 232:
1.1.1.2 ! ohara 233: if (rbitpos % GMP_NUMB_BITS != 0)
1.1 maekawa 234: {
235: mp_limb_t savelimb, rcy;
236: /* Target of of new chunk is not bit aligned. Use temp space
237: and align things by shifting it up. */
238: lc (tp, rstate);
239: savelimb = r2p[0];
1.1.1.2 ! ohara 240: rcy = mpn_lshift (r2p, tp, tn, rbitpos % GMP_NUMB_BITS);
1.1 maekawa 241: r2p[0] |= savelimb;
1.1.1.2 ! ohara 242: /* bogus */ if ((chunk_nbits % GMP_NUMB_BITS + rbitpos % GMP_NUMB_BITS)
! 243: > GMP_NUMB_BITS)
1.1 maekawa 244: r2p[tn] = rcy;
245: }
246: else
247: {
248: /* Target of of new chunk is bit aligned. Let `lc' put bits
249: directly into our target variable. */
250: lc (r2p, rstate);
251: }
252: rbitpos += chunk_nbits;
253: }
254:
255: /* Handle last [0..chunk_nbits) bits. */
256: if (rbitpos != nbits)
257: {
1.1.1.2 ! ohara 258: mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS;
1.1 maekawa 259: int last_nbits = nbits - rbitpos;
1.1.1.2 ! ohara 260: tn = (last_nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;
1.1 maekawa 261: lc (tp, rstate);
1.1.1.2 ! ohara 262: if (rbitpos % GMP_NUMB_BITS != 0)
1.1 maekawa 263: {
264: mp_limb_t savelimb, rcy;
265: /* Target of of new chunk is not bit aligned. Use temp space
266: and align things by shifting it up. */
267: savelimb = r2p[0];
1.1.1.2 ! ohara 268: rcy = mpn_lshift (r2p, tp, tn, rbitpos % GMP_NUMB_BITS);
1.1 maekawa 269: r2p[0] |= savelimb;
1.1.1.2 ! ohara 270: if (rbitpos + tn * GMP_NUMB_BITS - rbitpos % GMP_NUMB_BITS < nbits)
1.1 maekawa 271: r2p[tn] = rcy;
272: }
273: else
274: {
275: MPN_COPY (r2p, tp, tn);
276: }
277: /* Mask off top bits if needed. */
1.1.1.2 ! ohara 278: if (nbits % GMP_NUMB_BITS != 0)
! 279: rp[nbits / GMP_NUMB_BITS]
! 280: &= ~ ((~(mp_limb_t) 0) << nbits % GMP_NUMB_BITS);
1.1 maekawa 281: }
282:
283: TMP_FREE (lcmark);
284: break;
285: }
286:
287: default:
1.1.1.2 ! ohara 288: ASSERT (0);
1.1 maekawa 289: break;
290: }
291: }
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