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