randlc2x.c
上传用户:qaz666999
上传日期:2022-08-06
资源大小:2570k
文件大小:9k
- /* 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 = TMP_ALLOC_LIMBS (ta);
- MPN_ZERO (&tp[tmp], ta - tmp); /* mpn_mul won't zero it out. */
- }
- else
- tp = TMP_ALLOC_LIMBS (ta);
- /* 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 = TMP_ALLOC_LIMBS (tn);
- 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,
- mp_bitcnt_t 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;
- }