1 /* Linear Congruential pseudo-random number generator functions.
2
3 Copyright 1999, 2000, 2001, 2002, 2003, 2005 Free Software Foundation, Inc.
4
5 This file is part of the GNU MP Library.
6
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of the GNU Lesser General Public License as published by
9 the Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 The GNU MP Library is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
16
17 You should have received a copy of the GNU Lesser General Public License
18 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
19
20 #include "gmp.h"
21 #include "gmp-impl.h"
22
23
24 /* State structure for LC, the RNG_STATE() pointer in a gmp_randstate_t.
25
26 _mp_seed holds the current seed value, in the range 0 to 2^m2exp-1.
27 SIZ(_mp_seed) is fixed at BITS_TO_LIMBS(_mp_m2exp) and the value is
28 padded with high zero limbs if necessary. ALLOC(_mp_seed) is the current
29 size of PTR(_mp_seed) in the usual way. There only needs to be
30 BITS_TO_LIMBS(_mp_m2exp) allocated, but the mpz functions in the
31 initialization and seeding end up making it a bit more than this.
32
33 _mp_a is the "a" multiplier, in the range 0 to 2^m2exp-1. SIZ(_mp_a) is
34 the size of the value in the normal way for an mpz_t, except that a value
35 of zero is held with SIZ(_mp_a)==1 and PTR(_mp_a)[0]==0. This makes it
36 easy to call mpn_mul, and the case of a==0 is highly un-random and not
37 worth any trouble to optimize.
38
39 {_cp,_cn} is the "c" addend. Normally _cn is 1, but when nails are in
40 use a ulong can be bigger than one limb, and in this case _cn is 2 if
41 necessary. c==0 is stored as _cp[0]==0 and _cn==1, which makes it easy
42 to call __GMPN_ADD. c==0 is fairly un-random so isn't worth optimizing.
43
44 _mp_m2exp gives the modulus, namely 2^m2exp. We demand m2exp>=1, since
45 m2exp==0 would mean no bits at all out of each iteration, which makes no
46 sense. */
47
48 typedef struct {
49 mpz_t _mp_seed;
50 mpz_t _mp_a;
51 mp_size_t _cn;
52 mp_limb_t _cp[LIMBS_PER_ULONG];
53 unsigned long _mp_m2exp;
54 } gmp_rand_lc_struct;
55
56
57 /* lc (rp, state) -- Generate next number in LC sequence. Return the
58 number of valid bits in the result. Discards the lower half of the
59 result. */
60
61 static unsigned long int
lc(mp_ptr rp,gmp_randstate_t rstate)62 lc (mp_ptr rp, gmp_randstate_t rstate)
63 {
64 mp_ptr tp, seedp, ap;
65 mp_size_t ta;
66 mp_size_t tn, seedn, an;
67 unsigned long int m2exp;
68 unsigned long int bits;
69 int cy;
70 mp_size_t xn;
71 gmp_rand_lc_struct *p;
72 TMP_DECL;
73
74 p = (gmp_rand_lc_struct *) RNG_STATE (rstate);
75
76 m2exp = p->_mp_m2exp;
77
78 seedp = PTR (p->_mp_seed);
79 seedn = SIZ (p->_mp_seed);
80
81 ap = PTR (p->_mp_a);
82 an = SIZ (p->_mp_a);
83
84 /* Allocate temporary storage. Let there be room for calculation of
85 (A * seed + C) % M, or M if bigger than that. */
86
87 TMP_MARK;
88
89 ta = an + seedn + 1;
90 tn = BITS_TO_LIMBS (m2exp);
91 if (ta <= tn) /* that is, if (ta < tn + 1) */
92 {
93 mp_size_t tmp = an + seedn;
94 ta = tn + 1;
95 tp = TMP_ALLOC_LIMBS (ta);
96 MPN_ZERO (&tp[tmp], ta - tmp); /* mpn_mul won't zero it out. */
97 }
98 else
99 tp = TMP_ALLOC_LIMBS (ta);
100
101 /* t = a * seed. NOTE: an is always > 0; see initialization. */
102 ASSERT (seedn >= an && an > 0);
103 mpn_mul (tp, seedp, seedn, ap, an);
104
105 /* t = t + c. NOTE: tn is always >= p->_cn (precondition for __GMPN_ADD);
106 see initialization. */
107 ASSERT (tn >= p->_cn);
108 __GMPN_ADD (cy, tp, tp, tn, p->_cp, p->_cn);
109
110 /* t = t % m */
111 tp[m2exp / GMP_NUMB_BITS] &= (CNST_LIMB (1) << m2exp % GMP_NUMB_BITS) - 1;
112
113 /* Save result as next seed. */
114 MPN_COPY (PTR (p->_mp_seed), tp, tn);
115
116 /* Discard the lower m2exp/2 of the result. */
117 bits = m2exp / 2;
118 xn = bits / GMP_NUMB_BITS;
119
120 tn -= xn;
121 if (tn > 0)
122 {
123 unsigned int cnt = bits % GMP_NUMB_BITS;
124 if (cnt != 0)
125 {
126 mpn_rshift (tp, tp + xn, tn, cnt);
127 MPN_COPY_INCR (rp, tp, xn + 1);
128 }
129 else /* Even limb boundary. */
130 MPN_COPY_INCR (rp, tp + xn, tn);
131 }
132
133 TMP_FREE;
134
135 /* Return number of valid bits in the result. */
136 return (m2exp + 1) / 2;
137 }
138
139
140 /* Obtain a sequence of random numbers. */
141 static void
randget_lc(gmp_randstate_t rstate,mp_ptr rp,unsigned long int nbits)142 randget_lc (gmp_randstate_t rstate, mp_ptr rp, unsigned long int nbits)
143 {
144 unsigned long int rbitpos;
145 int chunk_nbits;
146 mp_ptr tp;
147 mp_size_t tn;
148 gmp_rand_lc_struct *p;
149 TMP_DECL;
150
151 p = (gmp_rand_lc_struct *) RNG_STATE (rstate);
152
153 TMP_MARK;
154
155 chunk_nbits = p->_mp_m2exp / 2;
156 tn = BITS_TO_LIMBS (chunk_nbits);
157
158 tp = TMP_ALLOC_LIMBS (tn);
159
160 rbitpos = 0;
161 while (rbitpos + chunk_nbits <= nbits)
162 {
163 mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS;
164
165 if (rbitpos % GMP_NUMB_BITS != 0)
166 {
167 mp_limb_t savelimb, rcy;
168 /* Target of new chunk is not bit aligned. Use temp space
169 and align things by shifting it up. */
170 lc (tp, rstate);
171 savelimb = r2p[0];
172 rcy = mpn_lshift (r2p, tp, tn, rbitpos % GMP_NUMB_BITS);
173 r2p[0] |= savelimb;
174 /* bogus */
175 if ((chunk_nbits % GMP_NUMB_BITS + rbitpos % GMP_NUMB_BITS)
176 > GMP_NUMB_BITS)
177 r2p[tn] = rcy;
178 }
179 else
180 {
181 /* Target of new chunk is bit aligned. Let `lc' put bits
182 directly into our target variable. */
183 lc (r2p, rstate);
184 }
185 rbitpos += chunk_nbits;
186 }
187
188 /* Handle last [0..chunk_nbits) bits. */
189 if (rbitpos != nbits)
190 {
191 mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS;
192 int last_nbits = nbits - rbitpos;
193 tn = BITS_TO_LIMBS (last_nbits);
194 lc (tp, rstate);
195 if (rbitpos % GMP_NUMB_BITS != 0)
196 {
197 mp_limb_t savelimb, rcy;
198 /* Target of new chunk is not bit aligned. Use temp space
199 and align things by shifting it up. */
200 savelimb = r2p[0];
201 rcy = mpn_lshift (r2p, tp, tn, rbitpos % GMP_NUMB_BITS);
202 r2p[0] |= savelimb;
203 if (rbitpos + tn * GMP_NUMB_BITS - rbitpos % GMP_NUMB_BITS < nbits)
204 r2p[tn] = rcy;
205 }
206 else
207 {
208 MPN_COPY (r2p, tp, tn);
209 }
210 /* Mask off top bits if needed. */
211 if (nbits % GMP_NUMB_BITS != 0)
212 rp[nbits / GMP_NUMB_BITS]
213 &= ~(~CNST_LIMB (0) << nbits % GMP_NUMB_BITS);
214 }
215
216 TMP_FREE;
217 }
218
219
220 static void
randseed_lc(gmp_randstate_t rstate,mpz_srcptr seed)221 randseed_lc (gmp_randstate_t rstate, mpz_srcptr seed)
222 {
223 gmp_rand_lc_struct *p = (gmp_rand_lc_struct *) RNG_STATE (rstate);
224 mpz_ptr seedz = p->_mp_seed;
225 mp_size_t seedn = BITS_TO_LIMBS (p->_mp_m2exp);
226
227 /* Store p->_mp_seed as an unnormalized integer with size enough
228 for numbers up to 2^m2exp-1. That size can't be zero. */
229 mpz_fdiv_r_2exp (seedz, seed, p->_mp_m2exp);
230 MPN_ZERO (&PTR (seedz)[SIZ (seedz)], seedn - SIZ (seedz));
231 SIZ (seedz) = seedn;
232 }
233
234
235 static void
randclear_lc(gmp_randstate_t rstate)236 randclear_lc (gmp_randstate_t rstate)
237 {
238 gmp_rand_lc_struct *p = (gmp_rand_lc_struct *) RNG_STATE (rstate);
239
240 mpz_clear (p->_mp_seed);
241 mpz_clear (p->_mp_a);
242 (*__gmp_free_func) (p, sizeof (gmp_rand_lc_struct));
243 }
244
245 static void randiset_lc __GMP_PROTO ((gmp_randstate_ptr dst, gmp_randstate_srcptr src));
246
247 static const gmp_randfnptr_t Linear_Congruential_Generator = {
248 randseed_lc,
249 randget_lc,
250 randclear_lc,
251 randiset_lc
252 };
253
254 static void
randiset_lc(gmp_randstate_ptr dst,gmp_randstate_srcptr src)255 randiset_lc (gmp_randstate_ptr dst, gmp_randstate_srcptr src)
256 {
257 gmp_rand_lc_struct *dstp, *srcp;
258
259 srcp = (gmp_rand_lc_struct *) RNG_STATE (src);
260 dstp = (*__gmp_allocate_func) (sizeof (gmp_rand_lc_struct));
261
262 RNG_STATE (dst) = (void *) dstp;
263 RNG_FNPTR (dst) = (void *) &Linear_Congruential_Generator;
264
265 /* _mp_seed and _mp_a might be unnormalized (high zero limbs), but
266 mpz_init_set won't worry about that */
267 mpz_init_set (dstp->_mp_seed, srcp->_mp_seed);
268 mpz_init_set (dstp->_mp_a, srcp->_mp_a);
269
270 dstp->_cn = srcp->_cn;
271
272 dstp->_cp[0] = srcp->_cp[0];
273 if (LIMBS_PER_ULONG > 1)
274 dstp->_cp[1] = srcp->_cp[1];
275 if (LIMBS_PER_ULONG > 2) /* usually there's only 1 or 2 */
276 MPN_COPY (dstp->_cp + 2, srcp->_cp + 2, LIMBS_PER_ULONG - 2);
277
278 dstp->_mp_m2exp = srcp->_mp_m2exp;
279 }
280
281
282 void
gmp_randinit_lc_2exp(gmp_randstate_t rstate,mpz_srcptr a,unsigned long int c,mp_bitcnt_t m2exp)283 gmp_randinit_lc_2exp (gmp_randstate_t rstate,
284 mpz_srcptr a,
285 unsigned long int c,
286 mp_bitcnt_t m2exp)
287 {
288 gmp_rand_lc_struct *p;
289 mp_size_t seedn = BITS_TO_LIMBS (m2exp);
290
291 ASSERT_ALWAYS (m2exp != 0);
292
293 p = __GMP_ALLOCATE_FUNC_TYPE (1, gmp_rand_lc_struct);
294 RNG_STATE (rstate) = (void *) p;
295 RNG_FNPTR (rstate) = (void *) &Linear_Congruential_Generator;
296
297 /* allocate m2exp bits of space for p->_mp_seed, and initial seed "1" */
298 mpz_init2 (p->_mp_seed, m2exp);
299 MPN_ZERO (PTR (p->_mp_seed), seedn);
300 SIZ (p->_mp_seed) = seedn;
301 PTR (p->_mp_seed)[0] = 1;
302
303 /* "a", forced to 0 to 2^m2exp-1 */
304 mpz_init (p->_mp_a);
305 mpz_fdiv_r_2exp (p->_mp_a, a, m2exp);
306
307 /* Avoid SIZ(a) == 0 to avoid checking for special case in lc(). */
308 if (SIZ (p->_mp_a) == 0)
309 {
310 SIZ (p->_mp_a) = 1;
311 PTR (p->_mp_a)[0] = CNST_LIMB (0);
312 }
313
314 MPN_SET_UI (p->_cp, p->_cn, c);
315
316 /* Internally we may discard any bits of c above m2exp. The following
317 code ensures that __GMPN_ADD in lc() will always work. */
318 if (seedn < p->_cn)
319 p->_cn = (p->_cp[0] != 0);
320
321 p->_mp_m2exp = m2exp;
322 }
323