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

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

Line | |
---|---|

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 1999, 2000, 2001, 2002, 2004 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 | |

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 | |

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 | |

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

113 | |

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); |

117 | |

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 | |

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) |

127 | { |

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

129 | with C smaller than M*/ |

130 | *rp = c; |

131 | |

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 | |

142 | *seedp = c; |

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

144 | return m2exp; |

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) |

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

157 | else |

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

159 | tn = an + seedn; |

160 | |

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

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

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

164 | |

165 | if (tn > m2exp / GMP_NUMB_BITS) |

166 | { |

167 | /* t = t % m */ |

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; |

170 | } |

171 | |

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

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

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

175 | |

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; |

180 | |

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 | } |

194 | |

195 | TMP_FREE (mark); |

196 | |

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

198 | return (m2exp + 1) / 2; |

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 | |

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

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

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 | |

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

232 | |

233 | switch (rstate->_mp_alg) |

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 | |

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

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

247 | |

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

249 | |

250 | rbitpos = 0; |

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

252 | { |

253 | mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS; |

254 | |

255 | if (rbitpos % GMP_NUMB_BITS != 0) |

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]; |

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

263 | r2p[0] |= savelimb; |

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

265 | > GMP_NUMB_BITS) |

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 | { |

280 | mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS; |

281 | int last_nbits = nbits - rbitpos; |

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

283 | lc (tp, rstate); |

284 | if (rbitpos % GMP_NUMB_BITS != 0) |

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]; |

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

291 | r2p[0] |= savelimb; |

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

293 | r2p[tn] = rcy; |

294 | } |

295 | else |

296 | { |

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

298 | } |

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

300 | if (nbits % GMP_NUMB_BITS != 0) |

301 | rp[nbits / GMP_NUMB_BITS] |

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

303 | } |

304 | |

305 | TMP_FREE (lcmark); |

306 | break; |

307 | } |

308 | |

309 | default: |

310 | ASSERT (0); |

311 | break; |

312 | } |

313 | } |

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