#
source:
trunk/third/gmp/randraw.c
@
22254

Revision 22254, 8.4 KB checked in by ghudson, 19 years ago (diff) |
---|

Rev | Line | |
---|---|---|

[15293] | 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 | ||

[22253] | 5 | Copyright 1999, 2000, 2001, 2002, 2004 Free Software Foundation, Inc. |

[15293] | 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 | ||

[22253] | 101 | /* Zero out the limbs _gmp_rand below expects us to write. This is a hack |

102 | to cover the seedn==0 case, and in case tn < cn due to small "a" and | |

103 | seed. (Incidentally, the "return m2exp" for the seedn==0 case is | |

104 | bogus, _gmp_rand ignores the return, it expects and looks at just | |

105 | "m2exp/2" always.) */ | |

106 | { | |

107 | int chunk_nbits = rstate->_mp_algdata._mp_lc->_mp_m2exp / 2; | |

108 | mp_size_t cn = (chunk_nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; | |

109 | MPN_ZERO (rp, cn); | |

110 | } | |

111 | ||

[18190] | 112 | m2exp = rstate->_mp_algdata._mp_lc->_mp_m2exp; |

[15293] | 113 | |

[18190] | 114 | /* The code below assumes the mod part is a power of two. Make sure |

115 | that is the case. */ | |

116 | ASSERT_ALWAYS (m2exp != 0); | |

[15293] | 117 | |

[18190] | 118 | c = (mp_limb_t) rstate->_mp_algdata._mp_lc->_mp_c; |

119 | ||

120 | seedp = PTR (rstate->_mp_seed); | |

121 | seedn = SIZ (rstate->_mp_seed); | |

122 | ||

[22253] | 123 | ap = PTR (rstate->_mp_algdata._mp_lc->_mp_a); |

124 | an = SIZ (rstate->_mp_algdata._mp_lc->_mp_a); | |

125 | ||

126 | if (seedn == 0 || an == 0) | |

[15293] | 127 | { |

[18190] | 128 | /* Seed is 0. Result is C % M. Assume table is sensibly stored, |

129 | with C smaller than M*/ | |

[15293] | 130 | *rp = c; |

131 | ||

[22253] | 132 | /* Discard the lower m2exp/2 bits of result. */ |

133 | { | |

134 | unsigned long int bits = m2exp / 2; | |

135 | mp_size_t xn = bits / GMP_NUMB_BITS; | |

136 | if (bits >= GMP_LIMB_BITS) | |

137 | *rp = 0; | |

138 | else | |

139 | *rp >>= bits; | |

140 | } | |

141 | ||

[18190] | 142 | *seedp = c; |

143 | SIZ (rstate->_mp_seed) = 1; | |

144 | return m2exp; | |

[15293] | 145 | } |

146 | ||

147 | /* Allocate temporary storage. Let there be room for calculation of | |

148 | (A * seed + C) % M, or M if bigger than that. */ | |

149 | ||

150 | TMP_MARK (mark); | |

151 | ta = an + seedn + 1; | |

152 | tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB); | |

153 | ||

154 | /* t = a * seed */ | |

155 | if (seedn >= an) | |

[18190] | 156 | mpn_mul (tp, seedp, seedn, ap, an); |

[15293] | 157 | else |

[18190] | 158 | mpn_mul (tp, ap, an, seedp, seedn); |

[15293] | 159 | tn = an + seedn; |

160 | ||

161 | /* t = t + c */ | |

[18190] | 162 | tp[tn] = 0; /* sentinel, stops MPN_INCR_U */ |

163 | MPN_INCR_U (tp, tn, c); | |

[15293] | 164 | |

[22253] | 165 | if (tn > m2exp / GMP_NUMB_BITS) |

166 | { | |

[15293] | 167 | /* t = t % m */ |

[18190] | 168 | tp[m2exp / GMP_NUMB_BITS] &= ((mp_limb_t) 1 << m2exp % GMP_NUMB_BITS) - 1; |

169 | tn = (m2exp + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; | |

[22253] | 170 | } |

[15293] | 171 | |

172 | /* Save result as next seed. */ | |

[18190] | 173 | MPN_COPY (PTR (rstate->_mp_seed), tp, tn); |

174 | SIZ (rstate->_mp_seed) = tn; | |

[15293] | 175 | |

[18190] | 176 | { |

177 | /* Discard the lower m2exp/2 bits of result. */ | |

178 | unsigned long int bits = m2exp / 2; | |

179 | mp_size_t xn = bits / GMP_NUMB_BITS; | |

[15293] | 180 | |

[18190] | 181 | tn -= xn; |

182 | if (tn > 0) | |

183 | { | |

184 | unsigned int cnt = bits % GMP_NUMB_BITS; | |

185 | if (cnt != 0) | |

186 | { | |

187 | mpn_rshift (tp, tp + xn, tn, cnt); | |

188 | MPN_COPY_INCR (rp, tp, xn + 1); | |

189 | } | |

190 | else /* Even limb boundary. */ | |

191 | MPN_COPY_INCR (rp, tp + xn, tn); | |

192 | } | |

193 | } | |

[15293] | 194 | |

195 | TMP_FREE (mark); | |

196 | ||

197 | /* Return number of valid bits in the result. */ | |

[18190] | 198 | return (m2exp + 1) / 2; |

[15293] | 199 | } |

200 | ||

201 | #ifdef RAWRANDEBUG | |

202 | /* Set even bits to EVENBITS and odd bits to ! EVENBITS in RP. | |

203 | Number of bits is m2exp in state. */ | |

204 | /* FIXME: Remove. */ | |

205 | unsigned long int | |

206 | lc_test (mp_ptr rp, gmp_randstate_t s, const int evenbits) | |

207 | { | |

208 | unsigned long int rn, nbits; | |

209 | int f; | |

210 | ||

[18190] | 211 | nbits = s->_mp_algdata._mp_lc->_mp_m2exp / 2; |

212 | rn = nbits / GMP_NUMB_BITS + (nbits % GMP_NUMB_BITS != 0); | |

[15293] | 213 | MPN_ZERO (rp, rn); |

214 | ||

215 | for (f = 0; f < nbits; f++) | |

216 | { | |

217 | mpn_lshift (rp, rp, rn, 1); | |

218 | if (f % 2 == ! evenbits) | |

219 | rp[0] += 1; | |

220 | } | |

221 | ||

222 | return nbits; | |

223 | } | |

224 | #endif /* RAWRANDEBUG */ | |

225 | ||

226 | void | |

227 | _gmp_rand (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits) | |

228 | { | |

229 | mp_size_t rn; /* Size of R. */ | |

230 | ||

[18190] | 231 | rn = (nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; |

[15293] | 232 | |

[18190] | 233 | switch (rstate->_mp_alg) |

[15293] | 234 | { |

235 | case GMP_RAND_ALG_LC: | |

236 | { | |

237 | unsigned long int rbitpos; | |

238 | int chunk_nbits; | |

239 | mp_ptr tp; | |

240 | mp_size_t tn; | |

241 | TMP_DECL (lcmark); | |

242 | ||

243 | TMP_MARK (lcmark); | |

244 | ||

[18190] | 245 | chunk_nbits = rstate->_mp_algdata._mp_lc->_mp_m2exp / 2; |

246 | tn = (chunk_nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; | |

[15293] | 247 | |

248 | tp = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB); | |

249 | ||

250 | rbitpos = 0; | |

251 | while (rbitpos + chunk_nbits <= nbits) | |

252 | { | |

[18190] | 253 | mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS; |

[15293] | 254 | |

[18190] | 255 | if (rbitpos % GMP_NUMB_BITS != 0) |

[15293] | 256 | { |

257 | mp_limb_t savelimb, rcy; | |

258 | /* Target of of new chunk is not bit aligned. Use temp space | |

259 | and align things by shifting it up. */ | |

260 | lc (tp, rstate); | |

261 | savelimb = r2p[0]; | |

[18190] | 262 | rcy = mpn_lshift (r2p, tp, tn, rbitpos % GMP_NUMB_BITS); |

[15293] | 263 | r2p[0] |= savelimb; |

[18190] | 264 | /* bogus */ if ((chunk_nbits % GMP_NUMB_BITS + rbitpos % GMP_NUMB_BITS) |

265 | > GMP_NUMB_BITS) | |

[15293] | 266 | r2p[tn] = rcy; |

267 | } | |

268 | else | |

269 | { | |

270 | /* Target of of new chunk is bit aligned. Let `lc' put bits | |

271 | directly into our target variable. */ | |

272 | lc (r2p, rstate); | |

273 | } | |

274 | rbitpos += chunk_nbits; | |

275 | } | |

276 | ||

277 | /* Handle last [0..chunk_nbits) bits. */ | |

278 | if (rbitpos != nbits) | |

279 | { | |

[18190] | 280 | mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS; |

[15293] | 281 | int last_nbits = nbits - rbitpos; |

[18190] | 282 | tn = (last_nbits + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS; |

[15293] | 283 | lc (tp, rstate); |

[18190] | 284 | if (rbitpos % GMP_NUMB_BITS != 0) |

[15293] | 285 | { |

286 | mp_limb_t savelimb, rcy; | |

287 | /* Target of of new chunk is not bit aligned. Use temp space | |

288 | and align things by shifting it up. */ | |

289 | savelimb = r2p[0]; | |

[18190] | 290 | rcy = mpn_lshift (r2p, tp, tn, rbitpos % GMP_NUMB_BITS); |

[15293] | 291 | r2p[0] |= savelimb; |

[18190] | 292 | if (rbitpos + tn * GMP_NUMB_BITS - rbitpos % GMP_NUMB_BITS < nbits) |

[15293] | 293 | r2p[tn] = rcy; |

294 | } | |

295 | else | |

296 | { | |

297 | MPN_COPY (r2p, tp, tn); | |

298 | } | |

299 | /* Mask off top bits if needed. */ | |

[18190] | 300 | if (nbits % GMP_NUMB_BITS != 0) |

301 | rp[nbits / GMP_NUMB_BITS] | |

302 | &= ~ ((~(mp_limb_t) 0) << nbits % GMP_NUMB_BITS); | |

[15293] | 303 | } |

304 | ||

305 | TMP_FREE (lcmark); | |

306 | break; | |

307 | } | |

308 | ||

309 | default: | |

[18190] | 310 | ASSERT (0); |

[15293] | 311 | break; |

312 | } | |

313 | } |

**Note:**See TracBrowser for help on using the repository browser.