summaryrefslogtreecommitdiff
path: root/randlc2x.c
blob: f1994187d9141b6b63bce35b89622a32544bc38e (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
/* Linear Congruential pseudo-random number generator functions.

Copyright 1999, 2000, 2001, 2002, 2003, 2005 Free Software Foundation, Inc.

This file is part of the GNU MP Library.

The GNU MP Library is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 3 of the License, or (at your
option) any later version.

The GNU MP Library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
License for more details.

You should have received a copy of the GNU Lesser General Public License
along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */

#include "gmp.h"
#include "gmp-impl.h"


/* State structure for LC, the RNG_STATE() pointer in a gmp_randstate_t.

   _mp_seed holds the current seed value, in the range 0 to 2^m2exp-1.
   SIZ(_mp_seed) is fixed at BITS_TO_LIMBS(_mp_m2exp) and the value is
   padded with high zero limbs if necessary.  ALLOC(_mp_seed) is the current
   size of PTR(_mp_seed) in the usual way.  There only needs to be
   BITS_TO_LIMBS(_mp_m2exp) allocated, but the mpz functions in the
   initialization and seeding end up making it a bit more than this.

   _mp_a is the "a" multiplier, in the range 0 to 2^m2exp-1.  SIZ(_mp_a) is
   the size of the value in the normal way for an mpz_t, except that a value
   of zero is held with SIZ(_mp_a)==1 and PTR(_mp_a)[0]==0.  This makes it
   easy to call mpn_mul, and the case of a==0 is highly un-random and not
   worth any trouble to optimize.

   {_cp,_cn} is the "c" addend.  Normally _cn is 1, but when nails are in
   use a ulong can be bigger than one limb, and in this case _cn is 2 if
   necessary.  c==0 is stored as _cp[0]==0 and _cn==1, which makes it easy
   to call __GMPN_ADD.  c==0 is fairly un-random so isn't worth optimizing.

   _mp_m2exp gives the modulus, namely 2^m2exp.  We demand m2exp>=1, since
   m2exp==0 would mean no bits at all out of each iteration, which makes no
   sense.  */

typedef struct {
  mpz_t          _mp_seed;
  mpz_t          _mp_a;
  mp_size_t      _cn;
  mp_limb_t      _cp[LIMBS_PER_ULONG];
  unsigned long  _mp_m2exp;
} gmp_rand_lc_struct;


/* lc (rp, state) -- Generate next number in LC sequence.  Return the
   number of valid bits in the result.  Discards the lower half of the
   result.  */

static unsigned long int
lc (mp_ptr rp, gmp_randstate_t rstate)
{
  mp_ptr tp, seedp, ap;
  mp_size_t ta;
  mp_size_t tn, seedn, an;
  unsigned long int m2exp;
  unsigned long int bits;
  int cy;
  mp_size_t xn;
  gmp_rand_lc_struct *p;
  TMP_DECL;

  p = (gmp_rand_lc_struct *) RNG_STATE (rstate);

  m2exp = p->_mp_m2exp;

  seedp = PTR (p->_mp_seed);
  seedn = SIZ (p->_mp_seed);

  ap = PTR (p->_mp_a);
  an = SIZ (p->_mp_a);

  /* Allocate temporary storage.  Let there be room for calculation of
     (A * seed + C) % M, or M if bigger than that.  */

  TMP_MARK;

  ta = an + seedn + 1;
  tn = BITS_TO_LIMBS (m2exp);
  if (ta <= tn) /* that is, if (ta < tn + 1) */
    {
      mp_size_t tmp = an + seedn;
      ta = tn + 1;
      tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
      MPN_ZERO (&tp[tmp], ta - tmp); /* mpn_mul won't zero it out.  */
    }
  else
    tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);

  /* t = a * seed.  NOTE: an is always > 0; see initialization.  */
  ASSERT (seedn >= an && an > 0);
  mpn_mul (tp, seedp, seedn, ap, an);

  /* t = t + c.  NOTE: tn is always >= p->_cn (precondition for __GMPN_ADD);
     see initialization.  */
  ASSERT (tn >= p->_cn);
  __GMPN_ADD (cy, tp, tp, tn, p->_cp, p->_cn);

  /* t = t % m */
  tp[m2exp / GMP_NUMB_BITS] &= (CNST_LIMB (1) << m2exp % GMP_NUMB_BITS) - 1;

  /* Save result as next seed.  */
  MPN_COPY (PTR (p->_mp_seed), tp, tn);

  /* Discard the lower m2exp/2 of the result.  */
  bits = m2exp / 2;
  xn = bits / GMP_NUMB_BITS;

  tn -= xn;
  if (tn > 0)
    {
      unsigned int cnt = bits % GMP_NUMB_BITS;
      if (cnt != 0)
	{
	  mpn_rshift (tp, tp + xn, tn, cnt);
	  MPN_COPY_INCR (rp, tp, xn + 1);
	}
      else			/* Even limb boundary.  */
	MPN_COPY_INCR (rp, tp + xn, tn);
    }

  TMP_FREE;

  /* Return number of valid bits in the result.  */
  return (m2exp + 1) / 2;
}


/* Obtain a sequence of random numbers.  */
static void
randget_lc (gmp_randstate_t rstate, mp_ptr rp, unsigned long int nbits)
{
  unsigned long int rbitpos;
  int chunk_nbits;
  mp_ptr tp;
  mp_size_t tn;
  gmp_rand_lc_struct *p;
  TMP_DECL;

  p = (gmp_rand_lc_struct *) RNG_STATE (rstate);

  TMP_MARK;

  chunk_nbits = p->_mp_m2exp / 2;
  tn = BITS_TO_LIMBS (chunk_nbits);

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

  rbitpos = 0;
  while (rbitpos + chunk_nbits <= nbits)
    {
      mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS;

      if (rbitpos % GMP_NUMB_BITS != 0)
	{
	  mp_limb_t savelimb, rcy;
	  /* Target of new chunk is not bit aligned.  Use temp space
	     and align things by shifting it up.  */
	  lc (tp, rstate);
	  savelimb = r2p[0];
	  rcy = mpn_lshift (r2p, tp, tn, rbitpos % GMP_NUMB_BITS);
	  r2p[0] |= savelimb;
	  /* bogus */
	  if ((chunk_nbits % GMP_NUMB_BITS + rbitpos % GMP_NUMB_BITS)
	      > GMP_NUMB_BITS)
	    r2p[tn] = rcy;
	}
      else
	{
	  /* Target of new chunk is bit aligned.  Let `lc' put bits
	     directly into our target variable.  */
	  lc (r2p, rstate);
	}
      rbitpos += chunk_nbits;
    }

  /* Handle last [0..chunk_nbits) bits.  */
  if (rbitpos != nbits)
    {
      mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS;
      int last_nbits = nbits - rbitpos;
      tn = BITS_TO_LIMBS (last_nbits);
      lc (tp, rstate);
      if (rbitpos % GMP_NUMB_BITS != 0)
	{
	  mp_limb_t savelimb, rcy;
	  /* Target of new chunk is not bit aligned.  Use temp space
	     and align things by shifting it up.  */
	  savelimb = r2p[0];
	  rcy = mpn_lshift (r2p, tp, tn, rbitpos % GMP_NUMB_BITS);
	  r2p[0] |= savelimb;
	  if (rbitpos + tn * GMP_NUMB_BITS - rbitpos % GMP_NUMB_BITS < nbits)
	    r2p[tn] = rcy;
	}
      else
	{
	  MPN_COPY (r2p, tp, tn);
	}
      /* Mask off top bits if needed.  */
      if (nbits % GMP_NUMB_BITS != 0)
	rp[nbits / GMP_NUMB_BITS]
	  &= ~(~CNST_LIMB (0) << nbits % GMP_NUMB_BITS);
    }

  TMP_FREE;
}


static void
randseed_lc (gmp_randstate_t rstate, mpz_srcptr seed)
{
  gmp_rand_lc_struct *p = (gmp_rand_lc_struct *) RNG_STATE (rstate);
  mpz_ptr seedz = p->_mp_seed;
  mp_size_t seedn = BITS_TO_LIMBS (p->_mp_m2exp);

  /* Store p->_mp_seed as an unnormalized integer with size enough
     for numbers up to 2^m2exp-1.  That size can't be zero.  */
  mpz_fdiv_r_2exp (seedz, seed, p->_mp_m2exp);
  MPN_ZERO (&PTR (seedz)[SIZ (seedz)], seedn - SIZ (seedz));
  SIZ (seedz) = seedn;
}


static void
randclear_lc (gmp_randstate_t rstate)
{
  gmp_rand_lc_struct *p = (gmp_rand_lc_struct *) RNG_STATE (rstate);

  mpz_clear (p->_mp_seed);
  mpz_clear (p->_mp_a);
  (*__gmp_free_func) (p, sizeof (gmp_rand_lc_struct));
}

static void randiset_lc __GMP_PROTO ((gmp_randstate_ptr dst, gmp_randstate_srcptr src));

static const gmp_randfnptr_t Linear_Congruential_Generator = {
  randseed_lc,
  randget_lc,
  randclear_lc,
  randiset_lc
};

static void
randiset_lc (gmp_randstate_ptr dst, gmp_randstate_srcptr src)
{
  gmp_rand_lc_struct *dstp, *srcp;

  srcp = (gmp_rand_lc_struct *) RNG_STATE (src);
  dstp = (*__gmp_allocate_func) (sizeof (gmp_rand_lc_struct));

  RNG_STATE (dst) = (void *) dstp;
  RNG_FNPTR (dst) = (void *) &Linear_Congruential_Generator;

  /* _mp_seed and _mp_a might be unnormalized (high zero limbs), but
     mpz_init_set won't worry about that */
  mpz_init_set (dstp->_mp_seed, srcp->_mp_seed);
  mpz_init_set (dstp->_mp_a,    srcp->_mp_a);

  dstp->_cn = srcp->_cn;

  dstp->_cp[0] = srcp->_cp[0];
  if (LIMBS_PER_ULONG > 1)
    dstp->_cp[1] = srcp->_cp[1];
  if (LIMBS_PER_ULONG > 2)  /* usually there's only 1 or 2 */
    MPN_COPY (dstp->_cp + 2, srcp->_cp + 2, LIMBS_PER_ULONG - 2);

  dstp->_mp_m2exp = srcp->_mp_m2exp;
}


void
gmp_randinit_lc_2exp (gmp_randstate_t rstate,
		      mpz_srcptr a,
		      unsigned long int c,
		      unsigned long int m2exp)
{
  gmp_rand_lc_struct *p;
  mp_size_t seedn = BITS_TO_LIMBS (m2exp);

  ASSERT_ALWAYS (m2exp != 0);

  p = __GMP_ALLOCATE_FUNC_TYPE (1, gmp_rand_lc_struct);
  RNG_STATE (rstate) = (void *) p;
  RNG_FNPTR (rstate) = (void *) &Linear_Congruential_Generator;

  /* allocate m2exp bits of space for p->_mp_seed, and initial seed "1" */
  mpz_init2 (p->_mp_seed, m2exp);
  MPN_ZERO (PTR (p->_mp_seed), seedn);
  SIZ (p->_mp_seed) = seedn;
  PTR (p->_mp_seed)[0] = 1;

  /* "a", forced to 0 to 2^m2exp-1 */
  mpz_init (p->_mp_a);
  mpz_fdiv_r_2exp (p->_mp_a, a, m2exp);

  /* Avoid SIZ(a) == 0 to avoid checking for special case in lc().  */
  if (SIZ (p->_mp_a) == 0)
    {
      SIZ (p->_mp_a) = 1;
      PTR (p->_mp_a)[0] = CNST_LIMB (0);
    }

  MPN_SET_UI (p->_cp, p->_cn, c);

  /* Internally we may discard any bits of c above m2exp.  The following
     code ensures that __GMPN_ADD in lc() will always work.  */
  if (seedn < p->_cn)
    p->_cn = (p->_cp[0] != 0);

  p->_mp_m2exp = m2exp;
}