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

Revision 18191, 7.6 KB checked in by ghudson, 22 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 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 | m2exp = rstate->_mp_algdata._mp_lc->_mp_m2exp; |

102 | |

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

111 | |

112 | if (seedn == 0) |

113 | { |

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

115 | with C smaller than M*/ |

116 | *rp = c; |

117 | |

118 | *seedp = c; |

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

120 | return m2exp; |

121 | } |

122 | |

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

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

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

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

136 | else |

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

138 | tn = an + seedn; |

139 | |

140 | /* t = t + c */ |

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

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

143 | |

144 | ASSERT_ALWAYS (m2exp / GMP_NUMB_BITS < ta); |

145 | |

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

149 | |

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

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

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

153 | |

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

172 | |

173 | TMP_FREE (mark); |

174 | |

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

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

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

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

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

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

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

210 | |

211 | switch (rstate->_mp_alg) |

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

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

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

225 | |

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

227 | |

228 | rbitpos = 0; |

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

230 | { |

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

232 | |

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

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

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

241 | r2p[0] |= savelimb; |

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

243 | > GMP_NUMB_BITS) |

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

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

259 | int last_nbits = nbits - rbitpos; |

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

261 | lc (tp, rstate); |

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

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

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

269 | r2p[0] |= savelimb; |

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

271 | r2p[tn] = rcy; |

272 | } |

273 | else |

274 | { |

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

276 | } |

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

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

279 | rp[nbits / GMP_NUMB_BITS] |

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

281 | } |

282 | |

283 | TMP_FREE (lcmark); |

284 | break; |

285 | } |

286 | |

287 | default: |

288 | ASSERT (0); |

289 | break; |

290 | } |

291 | } |

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