1 /* Linear Congruential pseudo-random number generator functions. 2 3 Copyright 1999, 2000, 2001, 2002, 2003, 2005 Free Software Foundation, Inc. 4 5 This file is part of the GNU MP Library. 6 7 The GNU MP Library is free software; you can redistribute it and/or modify 8 it under the terms of the GNU Lesser General Public License as published by 9 the Free Software Foundation; either version 3 of the License, or (at your 10 option) any later version. 11 12 The GNU MP Library is distributed in the hope that it will be useful, but 13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 15 License for more details. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along 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 48 typedef 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 61 static unsigned long int 62 lc (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. */ 141 static void 142 randget_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 220 static void 221 randseed_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 235 static void 236 randclear_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 245 static void randiset_lc __GMP_PROTO ((gmp_randstate_ptr dst, gmp_randstate_srcptr src)); 246 247 static const gmp_randfnptr_t Linear_Congruential_Generator = { 248 randseed_lc, 249 randget_lc, 250 randclear_lc, 251 randiset_lc 252 }; 253 254 static void 255 randiset_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 282 void 283 gmp_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