--- /dev/null
+/* Linear Congruential pseudo-random number generator functions.
+
+Copyright 1999, 2000, 2001, 2002, 2003, 2005 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 3 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. If not, see http://www.gnu.org/licenses/. */
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
+
+/* State structure for LC, the RNG_STATE() pointer in a gmp_randstate_t.
+
+ _mp_seed holds the current seed value, in the range 0 to 2^m2exp-1.
+ SIZ(_mp_seed) is fixed at BITS_TO_LIMBS(_mp_m2exp) and the value is
+ padded with high zero limbs if necessary. ALLOC(_mp_seed) is the current
+ size of PTR(_mp_seed) in the usual way. There only needs to be
+ BITS_TO_LIMBS(_mp_m2exp) allocated, but the mpz functions in the
+ initialization and seeding end up making it a bit more than this.
+
+ _mp_a is the "a" multiplier, in the range 0 to 2^m2exp-1. SIZ(_mp_a) is
+ the size of the value in the normal way for an mpz_t, except that a value
+ of zero is held with SIZ(_mp_a)==1 and PTR(_mp_a)[0]==0. This makes it
+ easy to call mpn_mul, and the case of a==0 is highly un-random and not
+ worth any trouble to optimize.
+
+ {_cp,_cn} is the "c" addend. Normally _cn is 1, but when nails are in
+ use a ulong can be bigger than one limb, and in this case _cn is 2 if
+ necessary. c==0 is stored as _cp[0]==0 and _cn==1, which makes it easy
+ to call __GMPN_ADD. c==0 is fairly un-random so isn't worth optimizing.
+
+ _mp_m2exp gives the modulus, namely 2^m2exp. We demand m2exp>=1, since
+ m2exp==0 would mean no bits at all out of each iteration, which makes no
+ sense. */
+
+typedef struct {
+ mpz_t _mp_seed;
+ mpz_t _mp_a;
+ mp_size_t _cn;
+ mp_limb_t _cp[LIMBS_PER_ULONG];
+ unsigned long _mp_m2exp;
+} gmp_rand_lc_struct;
+
+
+/* lc (rp, state) -- Generate next number in LC sequence. Return the
+ number of valid bits in the result. Discards the lower half of the
+ result. */
+
+static unsigned long int
+lc (mp_ptr rp, gmp_randstate_t rstate)
+{
+ mp_ptr tp, seedp, ap;
+ mp_size_t ta;
+ mp_size_t tn, seedn, an;
+ unsigned long int m2exp;
+ unsigned long int bits;
+ int cy;
+ mp_size_t xn;
+ gmp_rand_lc_struct *p;
+ TMP_DECL;
+
+ p = (gmp_rand_lc_struct *) RNG_STATE (rstate);
+
+ m2exp = p->_mp_m2exp;
+
+ seedp = PTR (p->_mp_seed);
+ seedn = SIZ (p->_mp_seed);
+
+ ap = PTR (p->_mp_a);
+ an = SIZ (p->_mp_a);
+
+ /* Allocate temporary storage. Let there be room for calculation of
+ (A * seed + C) % M, or M if bigger than that. */
+
+ TMP_MARK;
+
+ ta = an + seedn + 1;
+ tn = BITS_TO_LIMBS (m2exp);
+ if (ta <= tn) /* that is, if (ta < tn + 1) */
+ {
+ mp_size_t tmp = an + seedn;
+ ta = tn + 1;
+ tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
+ MPN_ZERO (&tp[tmp], ta - tmp); /* mpn_mul won't zero it out. */
+ }
+ else
+ tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
+
+ /* t = a * seed. NOTE: an is always > 0; see initialization. */
+ ASSERT (seedn >= an && an > 0);
+ mpn_mul (tp, seedp, seedn, ap, an);
+
+ /* t = t + c. NOTE: tn is always >= p->_cn (precondition for __GMPN_ADD);
+ see initialization. */
+ ASSERT (tn >= p->_cn);
+ __GMPN_ADD (cy, tp, tp, tn, p->_cp, p->_cn);
+
+ /* t = t % m */
+ tp[m2exp / GMP_NUMB_BITS] &= (CNST_LIMB (1) << m2exp % GMP_NUMB_BITS) - 1;
+
+ /* Save result as next seed. */
+ MPN_COPY (PTR (p->_mp_seed), tp, tn);
+
+ /* Discard the lower m2exp/2 of the result. */
+ bits = m2exp / 2;
+ xn = bits / GMP_NUMB_BITS;
+
+ tn -= xn;
+ if (tn > 0)
+ {
+ unsigned int cnt = bits % GMP_NUMB_BITS;
+ if (cnt != 0)
+ {
+ mpn_rshift (tp, tp + xn, tn, cnt);
+ MPN_COPY_INCR (rp, tp, xn + 1);
+ }
+ else /* Even limb boundary. */
+ MPN_COPY_INCR (rp, tp + xn, tn);
+ }
+
+ TMP_FREE;
+
+ /* Return number of valid bits in the result. */
+ return (m2exp + 1) / 2;
+}
+
+
+/* Obtain a sequence of random numbers. */
+static void
+randget_lc (gmp_randstate_t rstate, mp_ptr rp, unsigned long int nbits)
+{
+ unsigned long int rbitpos;
+ int chunk_nbits;
+ mp_ptr tp;
+ mp_size_t tn;
+ gmp_rand_lc_struct *p;
+ TMP_DECL;
+
+ p = (gmp_rand_lc_struct *) RNG_STATE (rstate);
+
+ TMP_MARK;
+
+ chunk_nbits = p->_mp_m2exp / 2;
+ tn = BITS_TO_LIMBS (chunk_nbits);
+
+ tp = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB);
+
+ rbitpos = 0;
+ while (rbitpos + chunk_nbits <= nbits)
+ {
+ mp_ptr r2p = rp + rbitpos / GMP_NUMB_BITS;
+
+ if (rbitpos % GMP_NUMB_BITS != 0)
+ {
+ mp_limb_t savelimb, rcy;
+ /* Target 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 % GMP_NUMB_BITS);
+ r2p[0] |= savelimb;
+ /* bogus */
+ if ((chunk_nbits % GMP_NUMB_BITS + rbitpos % GMP_NUMB_BITS)
+ > GMP_NUMB_BITS)
+ r2p[tn] = rcy;
+ }
+ else
+ {
+ /* Target 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 / GMP_NUMB_BITS;
+ int last_nbits = nbits - rbitpos;
+ tn = BITS_TO_LIMBS (last_nbits);
+ lc (tp, rstate);
+ if (rbitpos % GMP_NUMB_BITS != 0)
+ {
+ mp_limb_t savelimb, rcy;
+ /* Target 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 % GMP_NUMB_BITS);
+ r2p[0] |= savelimb;
+ if (rbitpos + tn * GMP_NUMB_BITS - rbitpos % GMP_NUMB_BITS < nbits)
+ r2p[tn] = rcy;
+ }
+ else
+ {
+ MPN_COPY (r2p, tp, tn);
+ }
+ /* Mask off top bits if needed. */
+ if (nbits % GMP_NUMB_BITS != 0)
+ rp[nbits / GMP_NUMB_BITS]
+ &= ~(~CNST_LIMB (0) << nbits % GMP_NUMB_BITS);
+ }
+
+ TMP_FREE;
+}
+
+
+static void
+randseed_lc (gmp_randstate_t rstate, mpz_srcptr seed)
+{
+ gmp_rand_lc_struct *p = (gmp_rand_lc_struct *) RNG_STATE (rstate);
+ mpz_ptr seedz = p->_mp_seed;
+ mp_size_t seedn = BITS_TO_LIMBS (p->_mp_m2exp);
+
+ /* Store p->_mp_seed as an unnormalized integer with size enough
+ for numbers up to 2^m2exp-1. That size can't be zero. */
+ mpz_fdiv_r_2exp (seedz, seed, p->_mp_m2exp);
+ MPN_ZERO (&PTR (seedz)[SIZ (seedz)], seedn - SIZ (seedz));
+ SIZ (seedz) = seedn;
+}
+
+
+static void
+randclear_lc (gmp_randstate_t rstate)
+{
+ gmp_rand_lc_struct *p = (gmp_rand_lc_struct *) RNG_STATE (rstate);
+
+ mpz_clear (p->_mp_seed);
+ mpz_clear (p->_mp_a);
+ (*__gmp_free_func) (p, sizeof (gmp_rand_lc_struct));
+}
+
+static void randiset_lc __GMP_PROTO ((gmp_randstate_ptr dst, gmp_randstate_srcptr src));
+
+static const gmp_randfnptr_t Linear_Congruential_Generator = {
+ randseed_lc,
+ randget_lc,
+ randclear_lc,
+ randiset_lc
+};
+
+static void
+randiset_lc (gmp_randstate_ptr dst, gmp_randstate_srcptr src)
+{
+ gmp_rand_lc_struct *dstp, *srcp;
+
+ srcp = (gmp_rand_lc_struct *) RNG_STATE (src);
+ dstp = (*__gmp_allocate_func) (sizeof (gmp_rand_lc_struct));
+
+ RNG_STATE (dst) = (void *) dstp;
+ RNG_FNPTR (dst) = (void *) &Linear_Congruential_Generator;
+
+ /* _mp_seed and _mp_a might be unnormalized (high zero limbs), but
+ mpz_init_set won't worry about that */
+ mpz_init_set (dstp->_mp_seed, srcp->_mp_seed);
+ mpz_init_set (dstp->_mp_a, srcp->_mp_a);
+
+ dstp->_cn = srcp->_cn;
+
+ dstp->_cp[0] = srcp->_cp[0];
+ if (LIMBS_PER_ULONG > 1)
+ dstp->_cp[1] = srcp->_cp[1];
+ if (LIMBS_PER_ULONG > 2) /* usually there's only 1 or 2 */
+ MPN_COPY (dstp->_cp + 2, srcp->_cp + 2, LIMBS_PER_ULONG - 2);
+
+ dstp->_mp_m2exp = srcp->_mp_m2exp;
+}
+
+
+void
+gmp_randinit_lc_2exp (gmp_randstate_t rstate,
+ mpz_srcptr a,
+ unsigned long int c,
+ unsigned long int m2exp)
+{
+ gmp_rand_lc_struct *p;
+ mp_size_t seedn = BITS_TO_LIMBS (m2exp);
+
+ ASSERT_ALWAYS (m2exp != 0);
+
+ p = __GMP_ALLOCATE_FUNC_TYPE (1, gmp_rand_lc_struct);
+ RNG_STATE (rstate) = (void *) p;
+ RNG_FNPTR (rstate) = (void *) &Linear_Congruential_Generator;
+
+ /* allocate m2exp bits of space for p->_mp_seed, and initial seed "1" */
+ mpz_init2 (p->_mp_seed, m2exp);
+ MPN_ZERO (PTR (p->_mp_seed), seedn);
+ SIZ (p->_mp_seed) = seedn;
+ PTR (p->_mp_seed)[0] = 1;
+
+ /* "a", forced to 0 to 2^m2exp-1 */
+ mpz_init (p->_mp_a);
+ mpz_fdiv_r_2exp (p->_mp_a, a, m2exp);
+
+ /* Avoid SIZ(a) == 0 to avoid checking for special case in lc(). */
+ if (SIZ (p->_mp_a) == 0)
+ {
+ SIZ (p->_mp_a) = 1;
+ PTR (p->_mp_a)[0] = CNST_LIMB (0);
+ }
+
+ MPN_SET_UI (p->_cp, p->_cn, c);
+
+ /* Internally we may discard any bits of c above m2exp. The following
+ code ensures that __GMPN_ADD in lc() will always work. */
+ if (seedn < p->_cn)
+ p->_cn = (p->_cp[0] != 0);
+
+ p->_mp_m2exp = m2exp;
+}
projects / msp430-gcc.git / blobdiff