source/libs/gmp/gmp-src/rand/randlc2x.c

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