contrib/gmp/randlc2x.c

   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
  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 = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
  96       MPN_ZERO (&tp[tmp], ta - tmp); /* mpn_mul won't zero it out.  */
  97     }
  98   else
  99     tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
 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
 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 = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB);
 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
 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
 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
 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
 283 gmp_randinit_lc_2exp (gmp_randstate_t rstate,
 284                       mpz_srcptr a,
 285                       unsigned long int c,
 286                       unsigned long int 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 }