rts/gmp/randraw.c

   1 /* _gmp_rand (rp, state, nbits) -- Generate a random bitstream of
   2    length NBITS in RP.  RP must have enough space allocated to hold
   3    NBITS.
   4
   5 Copyright (C) 1999, 2000 Free Software Foundation, Inc.
   6
   7 This file is part of the GNU MP Library.
   8
   9 The GNU MP Library is free software; you can redistribute it and/or modify
  10 it under the terms of the GNU Lesser General Public License as published by
  11 the Free Software Foundation; either version 2.1 of the License, or (at your
  12 option) any later version.
  13
  14 The GNU MP Library is distributed in the hope that it will be useful, but
  15 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  16 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  17 License for more details.
  18
  19 You should have received a copy of the GNU Lesser General Public License
  20 along with the GNU MP Library; see the file COPYING.LIB.  If not, write to
  21 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  22 MA 02111-1307, USA. */
  23
  24 #include "gmp.h"
  25 #include "gmp-impl.h"
  26 #include "longlong.h"
  27
  28 /* For linear congruential (LC), we use one of algorithms (1) or (2).
  29    (gmp-3.0 uses algorithm (1) with 'm' as a power of 2.)
  30
  31 LC algorithm (1).
  32
  33         X = (aX + c) mod m
  34
  35 [D. Knuth, "The Art of Computer Programming: Volume 2, Seminumerical Algorithms",
  36 Third Edition, Addison Wesley, 1998, pp. 184-185.]
  37
  38         X is the seed and the result
  39         a is chosen so that
  40             a mod 8 = 5 [3.2.1.2] and [3.2.1.3]
  41             .01m < a < .99m
  42             its binary or decimal digits is not a simple, regular pattern
  43             it has no large quotients when Euclid's algorithm is used to find
  44               gcd(a, m) [3.3.3]
  45             it passes the spectral test [3.3.4]
  46             it passes several tests of [3.3.2]
  47         c has no factor in common with m (c=1 or c=a can be good)
  48         m is large (2^30)
  49           is a power of 2 [3.2.1.1]
  50
  51 The least significant digits of the generated number are not very
  52 random.  It should be regarded as a random fraction X/m.  To get a
  53 random integer between 0 and n-1, multiply X/m by n and truncate.
  54 (Don't use X/n [ex 3.4.1-3])
  55
  56 The ``accuracy'' in t dimensions is one part in ``the t'th root of m'' [3.3.4].
  57
  58 Don't generate more than about m/1000 numbers without changing a, c, or m.
  59
  60 The sequence length depends on chosen a,c,m.
  61
  62
  63 LC algorithm (2).
  64
  65         X = a * (X mod q) - r * (long) (X/q)
  66         if X<0 then X+=m
  67
  68 [Knuth, pp. 185-186.]
  69
  70         X is the seed and the result
  71           as a seed is nonzero and less than m
  72         a is a primitive root of m (which means that a^2 <= m)
  73         q is (long) m / a
  74         r is m mod a
  75         m is a prime number near the largest easily computed integer
  76
  77 which gives
  78
  79         X = a * (X % ((long) m / a)) -
  80             (M % a) * ((long) (X / ((long) m / a)))
  81
  82 Since m is prime, the least-significant bits of X are just as random as
  83 the most-significant bits. */
  84
  85 /* Blum, Blum, and Shub.
  86
  87    [Bruce Schneier, "Applied Cryptography", Second Edition, John Wiley
  88    & Sons, Inc., 1996, pp. 417-418.]
  89
  90    "Find two large prime numbers, p and q, which are congruent to 3
  91    modulo 4.  The product of those numbers, n, is a blum integer.
  92    Choose another random integer, x, which is relatively prime to n.
  93    Compute
  94         x[0] = x^2 mod n
  95    That's the seed for the generator."
  96
  97    To generate a random bit, compute
  98         x[i] = x[i-1]^2 mod n
  99    The least significant bit of x[i] is the one we want.
 100
 101    We can use more than one bit from x[i], namely the
 102         log2(bitlength of x[i])
 103    least significant bits of x[i].
 104
 105    So, for a 32-bit seed we get 5 bits per computation.
 106
 107    The non-predictability of this generator is based on the difficulty
 108    of factoring n.
 109  */
 110
 111 /* -------------------------------------------------- */
 112
 113 /* lc (rp, state) -- Generate next number in LC sequence.  Return the
 114    number of valid bits in the result.  NOTE: If 'm' is a power of 2
 115    (m2exp != 0), discard the lower half of the result.  */
 116
 117 static
 118 unsigned long int
 119 #if __STDC__
 120 lc (mp_ptr rp, gmp_randstate_t rstate)
 121 #else
 122 lc (rp, rstate)
 123      mp_ptr rp;
 124      gmp_randstate_t rstate;
 125 #endif
 126 {
 127   mp_ptr tp, seedp, ap;
 128   mp_size_t ta;
 129   mp_size_t tn, seedn, an;
 130   mp_size_t retval;
 131   int shiftcount = 0;
 132   unsigned long int m2exp;
 133   mp_limb_t c;
 134   TMP_DECL (mark);
 135
 136   m2exp = rstate->algdata.lc->m2exp;
 137   c = (mp_limb_t) rstate->algdata.lc->c;
 138
 139   seedp = PTR (rstate->seed);
 140   seedn = SIZ (rstate->seed);
 141
 142   if (seedn == 0)
 143     {
 144       /* Seed is 0.  Result is C % M.  */
 145       *rp = c;
 146
 147       if (m2exp != 0)
 148         {
 149           /* M is a power of 2.  */
 150           if (m2exp < BITS_PER_MP_LIMB)
 151             {
 152               /* Only necessary when M may be smaller than C.  */
 153               *rp &= (((mp_limb_t) 1 << m2exp) - 1);
 154             }
 155         }
 156       else
 157         {
 158           /* M is not a power of 2.  */
 159           abort ();             /* FIXME.  */
 160         }
 161
 162       /* Save result as next seed.  */
 163       *seedp = *rp;
 164       SIZ (rstate->seed) = 1;
 165       return BITS_PER_MP_LIMB;
 166     }
 167
 168   ap = PTR (rstate->algdata.lc->a);
 169   an = SIZ (rstate->algdata.lc->a);
 170
 171   /* Allocate temporary storage.  Let there be room for calculation of
 172      (A * seed + C) % M, or M if bigger than that.  */
 173
 174   ASSERT_ALWAYS (m2exp != 0);   /* FIXME.  */
 175
 176   TMP_MARK (mark);
 177   ta = an + seedn + 1;
 178   tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
 179   MPN_ZERO (tp, ta);
 180
 181   /* t = a * seed */
 182   if (seedn >= an)
 183     mpn_mul_basecase (tp, seedp, seedn, ap, an);
 184   else
 185     mpn_mul_basecase (tp, ap, an, seedp, seedn);
 186   tn = an + seedn;
 187
 188   /* t = t + c */
 189   mpn_incr_u (tp, c);
 190
 191   /* t = t % m */
 192   if (m2exp != 0)
 193     {
 194       /* M is a power of 2.  The mod operation is trivial.  */
 195
 196       tp[m2exp / BITS_PER_MP_LIMB] &= ((mp_limb_t) 1 << m2exp % BITS_PER_MP_LIMB) - 1;
 197       tn = (m2exp + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
 198     }
 199   else
 200     {
 201       abort ();                 /* FIXME.  */
 202     }
 203
 204   /* Save result as next seed.  */
 205   MPN_COPY (PTR (rstate->seed), tp, tn);
 206   SIZ (rstate->seed) = tn;
 207
 208   if (m2exp != 0)
 209     {
 210       /* Discard the lower half of the result.  */
 211       unsigned long int discardb = m2exp / 2;
 212       mp_size_t discardl = discardb / BITS_PER_MP_LIMB;
 213
 214       tn -= discardl;
 215       if (tn > 0)
 216         {
 217           if (discardb % BITS_PER_MP_LIMB != 0)
 218             {
 219               mpn_rshift (tp, tp + discardl, tn, discardb % BITS_PER_MP_LIMB);
 220               MPN_COPY (rp, tp, (discardb + BITS_PER_MP_LIMB -1) / BITS_PER_MP_LIMB);
 221             }
 222           else                  /* Even limb boundary.  */
 223             MPN_COPY_INCR (rp, tp + discardl, tn);
 224         }
 225     }
 226   else
 227     {
 228       MPN_COPY (rp, tp, tn);
 229     }
 230
 231   TMP_FREE (mark);
 232
 233   /* Return number of valid bits in the result.  */
 234   if (m2exp != 0)
 235     retval = (m2exp + 1) / 2;
 236   else
 237     retval = SIZ (rstate->algdata.lc->m) * BITS_PER_MP_LIMB - shiftcount;
 238   return retval;
 239 }
 240
 241 #ifdef RAWRANDEBUG
 242 /* Set even bits to EVENBITS and odd bits to ! EVENBITS in RP.
 243    Number of bits is m2exp in state.  */
 244 /* FIXME: Remove.  */
 245 unsigned long int
 246 lc_test (mp_ptr rp, gmp_randstate_t s, const int evenbits)
 247 {
 248   unsigned long int rn, nbits;
 249   int f;
 250
 251   nbits = s->algdata.lc->m2exp / 2;
 252   rn = nbits / BITS_PER_MP_LIMB + (nbits % BITS_PER_MP_LIMB != 0);
 253   MPN_ZERO (rp, rn);
 254
 255   for (f = 0; f < nbits; f++)
 256     {
 257       mpn_lshift (rp, rp, rn, 1);
 258       if (f % 2 == ! evenbits)
 259         rp[0] += 1;
 260     }
 261
 262   return nbits;
 263 }
 264 #endif /* RAWRANDEBUG */
 265
 266 void
 267 #if __STDC__
 268 _gmp_rand (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
 269 #else
 270 _gmp_rand (rp, rstate, nbits)
 271      mp_ptr rp;
 272      gmp_randstate_t rstate;
 273      unsigned long int nbits;
 274 #endif
 275 {
 276   mp_size_t rn;                 /* Size of R.  */
 277
 278   rn = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
 279
 280   switch (rstate->alg)
 281     {
 282     case GMP_RAND_ALG_LC:
 283       {
 284         unsigned long int rbitpos;
 285         int chunk_nbits;
 286         mp_ptr tp;
 287         mp_size_t tn;
 288         TMP_DECL (lcmark);
 289
 290         TMP_MARK (lcmark);
 291
 292         chunk_nbits = rstate->algdata.lc->m2exp / 2;
 293         tn = (chunk_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
 294
 295         tp = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB);
 296
 297         rbitpos = 0;
 298         while (rbitpos + chunk_nbits <= nbits)
 299           {
 300             mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
 301
 302             if (rbitpos % BITS_PER_MP_LIMB != 0)
 303               {
 304                 mp_limb_t savelimb, rcy;
 305                 /* Target of of new chunk is not bit aligned.  Use temp space
 306                    and align things by shifting it up.  */
 307                 lc (tp, rstate);
 308                 savelimb = r2p[0];
 309                 rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
 310                 r2p[0] |= savelimb;
 311 /* bogus */     if ((chunk_nbits % BITS_PER_MP_LIMB + rbitpos % BITS_PER_MP_LIMB)
 312                     > BITS_PER_MP_LIMB)
 313                   r2p[tn] = rcy;
 314               }
 315             else
 316               {
 317                 /* Target of of new chunk is bit aligned.  Let `lc' put bits
 318                    directly into our target variable.  */
 319                 lc (r2p, rstate);
 320               }
 321             rbitpos += chunk_nbits;
 322           }
 323
 324         /* Handle last [0..chunk_nbits) bits.  */
 325         if (rbitpos != nbits)
 326           {
 327             mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
 328             int last_nbits = nbits - rbitpos;
 329             tn = (last_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
 330             lc (tp, rstate);
 331             if (rbitpos % BITS_PER_MP_LIMB != 0)
 332               {
 333                 mp_limb_t savelimb, rcy;
 334                 /* Target of of new chunk is not bit aligned.  Use temp space
 335                    and align things by shifting it up.  */
 336                 savelimb = r2p[0];
 337                 rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
 338                 r2p[0] |= savelimb;
 339                 if (rbitpos + tn * BITS_PER_MP_LIMB - rbitpos % BITS_PER_MP_LIMB < nbits)
 340                   r2p[tn] = rcy;
 341               }
 342             else
 343               {
 344                 MPN_COPY (r2p, tp, tn);
 345               }
 346             /* Mask off top bits if needed.  */
 347             if (nbits % BITS_PER_MP_LIMB != 0)
 348               rp[nbits / BITS_PER_MP_LIMB]
 349                 &= ~ ((~(mp_limb_t) 0) << nbits % BITS_PER_MP_LIMB);
 350           }
 351
 352         TMP_FREE (lcmark);
 353         break;
 354       }
 355
 356     default:
 357       gmp_errno |= GMP_ERROR_UNSUPPORTED_ARGUMENT;
 358       break;
 359     }
 360 }