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