summaryrefslogtreecommitdiff
path: root/rts/gmp/randraw.c
diff options
context:
space:
mode:
Diffstat (limited to 'rts/gmp/randraw.c')
-rw-r--r--rts/gmp/randraw.c360
1 files changed, 360 insertions, 0 deletions
diff --git a/rts/gmp/randraw.c b/rts/gmp/randraw.c
new file mode 100644
index 0000000000..c0c3889d33
--- /dev/null
+++ b/rts/gmp/randraw.c
@@ -0,0 +1,360 @@
+/* _gmp_rand (rp, state, nbits) -- Generate a random bitstream of
+ length NBITS in RP. RP must have enough space allocated to hold
+ NBITS.
+
+Copyright (C) 1999, 2000 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 2.1 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; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "longlong.h"
+
+/* For linear congruential (LC), we use one of algorithms (1) or (2).
+ (gmp-3.0 uses algorithm (1) with 'm' as a power of 2.)
+
+LC algorithm (1).
+
+ X = (aX + c) mod m
+
+[D. Knuth, "The Art of Computer Programming: Volume 2, Seminumerical Algorithms",
+Third Edition, Addison Wesley, 1998, pp. 184-185.]
+
+ X is the seed and the result
+ a is chosen so that
+ a mod 8 = 5 [3.2.1.2] and [3.2.1.3]
+ .01m < a < .99m
+ its binary or decimal digits is not a simple, regular pattern
+ it has no large quotients when Euclid's algorithm is used to find
+ gcd(a, m) [3.3.3]
+ it passes the spectral test [3.3.4]
+ it passes several tests of [3.3.2]
+ c has no factor in common with m (c=1 or c=a can be good)
+ m is large (2^30)
+ is a power of 2 [3.2.1.1]
+
+The least significant digits of the generated number are not very
+random. It should be regarded as a random fraction X/m. To get a
+random integer between 0 and n-1, multiply X/m by n and truncate.
+(Don't use X/n [ex 3.4.1-3])
+
+The ``accuracy'' in t dimensions is one part in ``the t'th root of m'' [3.3.4].
+
+Don't generate more than about m/1000 numbers without changing a, c, or m.
+
+The sequence length depends on chosen a,c,m.
+
+
+LC algorithm (2).
+
+ X = a * (X mod q) - r * (long) (X/q)
+ if X<0 then X+=m
+
+[Knuth, pp. 185-186.]
+
+ X is the seed and the result
+ as a seed is nonzero and less than m
+ a is a primitive root of m (which means that a^2 <= m)
+ q is (long) m / a
+ r is m mod a
+ m is a prime number near the largest easily computed integer
+
+which gives
+
+ X = a * (X % ((long) m / a)) -
+ (M % a) * ((long) (X / ((long) m / a)))
+
+Since m is prime, the least-significant bits of X are just as random as
+the most-significant bits. */
+
+/* Blum, Blum, and Shub.
+
+ [Bruce Schneier, "Applied Cryptography", Second Edition, John Wiley
+ & Sons, Inc., 1996, pp. 417-418.]
+
+ "Find two large prime numbers, p and q, which are congruent to 3
+ modulo 4. The product of those numbers, n, is a blum integer.
+ Choose another random integer, x, which is relatively prime to n.
+ Compute
+ x[0] = x^2 mod n
+ That's the seed for the generator."
+
+ To generate a random bit, compute
+ x[i] = x[i-1]^2 mod n
+ The least significant bit of x[i] is the one we want.
+
+ We can use more than one bit from x[i], namely the
+ log2(bitlength of x[i])
+ least significant bits of x[i].
+
+ So, for a 32-bit seed we get 5 bits per computation.
+
+ The non-predictability of this generator is based on the difficulty
+ of factoring n.
+ */
+
+/* -------------------------------------------------- */
+
+/* lc (rp, state) -- Generate next number in LC sequence. Return the
+ number of valid bits in the result. NOTE: If 'm' is a power of 2
+ (m2exp != 0), discard the lower half of the result. */
+
+static
+unsigned long int
+#if __STDC__
+lc (mp_ptr rp, gmp_randstate_t rstate)
+#else
+lc (rp, rstate)
+ mp_ptr rp;
+ gmp_randstate_t rstate;
+#endif
+{
+ mp_ptr tp, seedp, ap;
+ mp_size_t ta;
+ mp_size_t tn, seedn, an;
+ mp_size_t retval;
+ int shiftcount = 0;
+ unsigned long int m2exp;
+ mp_limb_t c;
+ TMP_DECL (mark);
+
+ m2exp = rstate->algdata.lc->m2exp;
+ c = (mp_limb_t) rstate->algdata.lc->c;
+
+ seedp = PTR (rstate->seed);
+ seedn = SIZ (rstate->seed);
+
+ if (seedn == 0)
+ {
+ /* Seed is 0. Result is C % M. */
+ *rp = c;
+
+ if (m2exp != 0)
+ {
+ /* M is a power of 2. */
+ if (m2exp < BITS_PER_MP_LIMB)
+ {
+ /* Only necessary when M may be smaller than C. */
+ *rp &= (((mp_limb_t) 1 << m2exp) - 1);
+ }
+ }
+ else
+ {
+ /* M is not a power of 2. */
+ abort (); /* FIXME. */
+ }
+
+ /* Save result as next seed. */
+ *seedp = *rp;
+ SIZ (rstate->seed) = 1;
+ return BITS_PER_MP_LIMB;
+ }
+
+ ap = PTR (rstate->algdata.lc->a);
+ an = SIZ (rstate->algdata.lc->a);
+
+ /* Allocate temporary storage. Let there be room for calculation of
+ (A * seed + C) % M, or M if bigger than that. */
+
+ ASSERT_ALWAYS (m2exp != 0); /* FIXME. */
+
+ TMP_MARK (mark);
+ ta = an + seedn + 1;
+ tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
+ MPN_ZERO (tp, ta);
+
+ /* t = a * seed */
+ if (seedn >= an)
+ mpn_mul_basecase (tp, seedp, seedn, ap, an);
+ else
+ mpn_mul_basecase (tp, ap, an, seedp, seedn);
+ tn = an + seedn;
+
+ /* t = t + c */
+ mpn_incr_u (tp, c);
+
+ /* t = t % m */
+ if (m2exp != 0)
+ {
+ /* M is a power of 2. The mod operation is trivial. */
+
+ tp[m2exp / BITS_PER_MP_LIMB] &= ((mp_limb_t) 1 << m2exp % BITS_PER_MP_LIMB) - 1;
+ tn = (m2exp + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
+ }
+ else
+ {
+ abort (); /* FIXME. */
+ }
+
+ /* Save result as next seed. */
+ MPN_COPY (PTR (rstate->seed), tp, tn);
+ SIZ (rstate->seed) = tn;
+
+ if (m2exp != 0)
+ {
+ /* Discard the lower half of the result. */
+ unsigned long int discardb = m2exp / 2;
+ mp_size_t discardl = discardb / BITS_PER_MP_LIMB;
+
+ tn -= discardl;
+ if (tn > 0)
+ {
+ if (discardb % BITS_PER_MP_LIMB != 0)
+ {
+ mpn_rshift (tp, tp + discardl, tn, discardb % BITS_PER_MP_LIMB);
+ MPN_COPY (rp, tp, (discardb + BITS_PER_MP_LIMB -1) / BITS_PER_MP_LIMB);
+ }
+ else /* Even limb boundary. */
+ MPN_COPY_INCR (rp, tp + discardl, tn);
+ }
+ }
+ else
+ {
+ MPN_COPY (rp, tp, tn);
+ }
+
+ TMP_FREE (mark);
+
+ /* Return number of valid bits in the result. */
+ if (m2exp != 0)
+ retval = (m2exp + 1) / 2;
+ else
+ retval = SIZ (rstate->algdata.lc->m) * BITS_PER_MP_LIMB - shiftcount;
+ return retval;
+}
+
+#ifdef RAWRANDEBUG
+/* Set even bits to EVENBITS and odd bits to ! EVENBITS in RP.
+ Number of bits is m2exp in state. */
+/* FIXME: Remove. */
+unsigned long int
+lc_test (mp_ptr rp, gmp_randstate_t s, const int evenbits)
+{
+ unsigned long int rn, nbits;
+ int f;
+
+ nbits = s->algdata.lc->m2exp / 2;
+ rn = nbits / BITS_PER_MP_LIMB + (nbits % BITS_PER_MP_LIMB != 0);
+ MPN_ZERO (rp, rn);
+
+ for (f = 0; f < nbits; f++)
+ {
+ mpn_lshift (rp, rp, rn, 1);
+ if (f % 2 == ! evenbits)
+ rp[0] += 1;
+ }
+
+ return nbits;
+}
+#endif /* RAWRANDEBUG */
+
+void
+#if __STDC__
+_gmp_rand (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
+#else
+_gmp_rand (rp, rstate, nbits)
+ mp_ptr rp;
+ gmp_randstate_t rstate;
+ unsigned long int nbits;
+#endif
+{
+ mp_size_t rn; /* Size of R. */
+
+ rn = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
+
+ switch (rstate->alg)
+ {
+ case GMP_RAND_ALG_LC:
+ {
+ unsigned long int rbitpos;
+ int chunk_nbits;
+ mp_ptr tp;
+ mp_size_t tn;
+ TMP_DECL (lcmark);
+
+ TMP_MARK (lcmark);
+
+ chunk_nbits = rstate->algdata.lc->m2exp / 2;
+ tn = (chunk_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
+
+ tp = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB);
+
+ rbitpos = 0;
+ while (rbitpos + chunk_nbits <= nbits)
+ {
+ mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
+
+ if (rbitpos % BITS_PER_MP_LIMB != 0)
+ {
+ mp_limb_t savelimb, rcy;
+ /* Target of 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 % BITS_PER_MP_LIMB);
+ r2p[0] |= savelimb;
+/* bogus */ if ((chunk_nbits % BITS_PER_MP_LIMB + rbitpos % BITS_PER_MP_LIMB)
+ > BITS_PER_MP_LIMB)
+ r2p[tn] = rcy;
+ }
+ else
+ {
+ /* Target of 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 / BITS_PER_MP_LIMB;
+ int last_nbits = nbits - rbitpos;
+ tn = (last_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
+ lc (tp, rstate);
+ if (rbitpos % BITS_PER_MP_LIMB != 0)
+ {
+ mp_limb_t savelimb, rcy;
+ /* Target of 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 % BITS_PER_MP_LIMB);
+ r2p[0] |= savelimb;
+ if (rbitpos + tn * BITS_PER_MP_LIMB - rbitpos % BITS_PER_MP_LIMB < nbits)
+ r2p[tn] = rcy;
+ }
+ else
+ {
+ MPN_COPY (r2p, tp, tn);
+ }
+ /* Mask off top bits if needed. */
+ if (nbits % BITS_PER_MP_LIMB != 0)
+ rp[nbits / BITS_PER_MP_LIMB]
+ &= ~ ((~(mp_limb_t) 0) << nbits % BITS_PER_MP_LIMB);
+ }
+
+ TMP_FREE (lcmark);
+ break;
+ }
+
+ default:
+ gmp_errno |= GMP_ERROR_UNSUPPORTED_ARGUMENT;
+ break;
+ }
+}
View Lorry Controller Status