1/* Linear Congruential pseudo-random number generator functions.
2
3Copyright 1999, 2000, 2001, 2002, 2003, 2005 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library.
6
7The GNU MP Library is free software; you can redistribute it and/or modify
8it under the terms of the GNU Lesser General Public License as published by
9the Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12The GNU MP Library is distributed in the hope that it will be useful, but
13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
15License for more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along 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
48typedef 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
61static unsigned long int
62lc (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.  */
141static void
142randget_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
220static void
221randseed_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
235static void
236randclear_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
245static void randiset_lc __GMP_PROTO ((gmp_randstate_ptr dst, gmp_randstate_srcptr src));
246
247static const gmp_randfnptr_t Linear_Congruential_Generator = {
248  randseed_lc,
249  randget_lc,
250  randclear_lc,
251  randiset_lc
252};
253
254static void
255randiset_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
282void
283gmp_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