Annotation of OpenXM_contrib/gmp/randraw.c, Revision 1.1
1.1 ! maekawa 1: /* _gmp_rand (rp, state, nbits) -- Generate a random bitstream of
! 2: length NBITS in RP. RP must have enough space allocated to hold
! 3: NBITS.
! 4:
! 5: Copyright (C) 1999, 2000 Free Software Foundation, Inc.
! 6:
! 7: This file is part of the GNU MP Library.
! 8:
! 9: The GNU MP Library is free software; you can redistribute it and/or modify
! 10: it under the terms of the GNU Lesser General Public License as published by
! 11: the Free Software Foundation; either version 2.1 of the License, or (at your
! 12: option) any later version.
! 13:
! 14: The GNU MP Library is distributed in the hope that it will be useful, but
! 15: WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
! 16: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
! 17: License for more details.
! 18:
! 19: You should have received a copy of the GNU Lesser General Public License
! 20: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
! 21: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
! 22: MA 02111-1307, USA. */
! 23:
! 24: #include "gmp.h"
! 25: #include "gmp-impl.h"
! 26: #include "longlong.h"
! 27:
! 28: /* For linear congruential (LC), we use one of algorithms (1) or (2).
! 29: (gmp-3.0 uses algorithm (1) with 'm' as a power of 2.)
! 30:
! 31: LC algorithm (1).
! 32:
! 33: X = (aX + c) mod m
! 34:
! 35: [D. Knuth, "The Art of Computer Programming: Volume 2, Seminumerical Algorithms",
! 36: Third Edition, Addison Wesley, 1998, pp. 184-185.]
! 37:
! 38: X is the seed and the result
! 39: a is chosen so that
! 40: a mod 8 = 5 [3.2.1.2] and [3.2.1.3]
! 41: .01m < a < .99m
! 42: its binary or decimal digits is not a simple, regular pattern
! 43: it has no large quotients when Euclid's algorithm is used to find
! 44: gcd(a, m) [3.3.3]
! 45: it passes the spectral test [3.3.4]
! 46: it passes several tests of [3.3.2]
! 47: c has no factor in common with m (c=1 or c=a can be good)
! 48: m is large (2^30)
! 49: is a power of 2 [3.2.1.1]
! 50:
! 51: The least significant digits of the generated number are not very
! 52: random. It should be regarded as a random fraction X/m. To get a
! 53: random integer between 0 and n-1, multiply X/m by n and truncate.
! 54: (Don't use X/n [ex 3.4.1-3])
! 55:
! 56: The ``accuracy'' in t dimensions is one part in ``the t'th root of m'' [3.3.4].
! 57:
! 58: Don't generate more than about m/1000 numbers without changing a, c, or m.
! 59:
! 60: The sequence length depends on chosen a,c,m.
! 61:
! 62:
! 63: LC algorithm (2).
! 64:
! 65: X = a * (X mod q) - r * (long) (X/q)
! 66: if X<0 then X+=m
! 67:
! 68: [Knuth, pp. 185-186.]
! 69:
! 70: X is the seed and the result
! 71: as a seed is nonzero and less than m
! 72: a is a primitive root of m (which means that a^2 <= m)
! 73: q is (long) m / a
! 74: r is m mod a
! 75: m is a prime number near the largest easily computed integer
! 76:
! 77: which gives
! 78:
! 79: X = a * (X % ((long) m / a)) -
! 80: (M % a) * ((long) (X / ((long) m / a)))
! 81:
! 82: Since m is prime, the least-significant bits of X are just as random as
! 83: the most-significant bits. */
! 84:
! 85: /* Blum, Blum, and Shub.
! 86:
! 87: [Bruce Schneier, "Applied Cryptography", Second Edition, John Wiley
! 88: & Sons, Inc., 1996, pp. 417-418.]
! 89:
! 90: "Find two large prime numbers, p and q, which are congruent to 3
! 91: modulo 4. The product of those numbers, n, is a blum integer.
! 92: Choose another random integer, x, which is relatively prime to n.
! 93: Compute
! 94: x[0] = x^2 mod n
! 95: That's the seed for the generator."
! 96:
! 97: To generate a random bit, compute
! 98: x[i] = x[i-1]^2 mod n
! 99: The least significant bit of x[i] is the one we want.
! 100:
! 101: We can use more than one bit from x[i], namely the
! 102: log2(bitlength of x[i])
! 103: least significant bits of x[i].
! 104:
! 105: So, for a 32-bit seed we get 5 bits per computation.
! 106:
! 107: The non-predictability of this generator is based on the difficulty
! 108: of factoring n.
! 109: */
! 110:
! 111: /* -------------------------------------------------- */
! 112:
! 113: /* lc (rp, state) -- Generate next number in LC sequence. Return the
! 114: number of valid bits in the result. NOTE: If 'm' is a power of 2
! 115: (m2exp != 0), discard the lower half of the result. */
! 116:
! 117: static
! 118: unsigned long int
! 119: #if __STDC__
! 120: lc (mp_ptr rp, gmp_randstate_t rstate)
! 121: #else
! 122: lc (rp, rstate)
! 123: mp_ptr rp;
! 124: gmp_randstate_t rstate;
! 125: #endif
! 126: {
! 127: mp_ptr tp, seedp, ap;
! 128: mp_size_t ta;
! 129: mp_size_t tn, seedn, an;
! 130: mp_size_t retval;
! 131: int shiftcount = 0;
! 132: unsigned long int m2exp;
! 133: mp_limb_t c;
! 134: TMP_DECL (mark);
! 135:
! 136: m2exp = rstate->algdata.lc->m2exp;
! 137: c = (mp_limb_t) rstate->algdata.lc->c;
! 138:
! 139: seedp = PTR (rstate->seed);
! 140: seedn = SIZ (rstate->seed);
! 141:
! 142: if (seedn == 0)
! 143: {
! 144: /* Seed is 0. Result is C % M. */
! 145: *rp = c;
! 146:
! 147: if (m2exp != 0)
! 148: {
! 149: /* M is a power of 2. */
! 150: if (m2exp < BITS_PER_MP_LIMB)
! 151: {
! 152: /* Only necessary when M may be smaller than C. */
! 153: *rp &= (((mp_limb_t) 1 << m2exp) - 1);
! 154: }
! 155: }
! 156: else
! 157: {
! 158: /* M is not a power of 2. */
! 159: abort (); /* FIXME. */
! 160: }
! 161:
! 162: /* Save result as next seed. */
! 163: *seedp = *rp;
! 164: SIZ (rstate->seed) = 1;
! 165: return BITS_PER_MP_LIMB;
! 166: }
! 167:
! 168: ap = PTR (rstate->algdata.lc->a);
! 169: an = SIZ (rstate->algdata.lc->a);
! 170:
! 171: /* Allocate temporary storage. Let there be room for calculation of
! 172: (A * seed + C) % M, or M if bigger than that. */
! 173:
! 174: ASSERT_ALWAYS (m2exp != 0); /* FIXME. */
! 175:
! 176: TMP_MARK (mark);
! 177: ta = an + seedn + 1;
! 178: tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);
! 179: MPN_ZERO (tp, ta);
! 180:
! 181: /* t = a * seed */
! 182: if (seedn >= an)
! 183: mpn_mul_basecase (tp, seedp, seedn, ap, an);
! 184: else
! 185: mpn_mul_basecase (tp, ap, an, seedp, seedn);
! 186: tn = an + seedn;
! 187:
! 188: /* t = t + c */
! 189: mpn_incr_u (tp, c);
! 190:
! 191: /* t = t % m */
! 192: if (m2exp != 0)
! 193: {
! 194: /* M is a power of 2. The mod operation is trivial. */
! 195:
! 196: tp[m2exp / BITS_PER_MP_LIMB] &= ((mp_limb_t) 1 << m2exp % BITS_PER_MP_LIMB) - 1;
! 197: tn = (m2exp + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
! 198: }
! 199: else
! 200: {
! 201: abort (); /* FIXME. */
! 202: }
! 203:
! 204: /* Save result as next seed. */
! 205: MPN_COPY (PTR (rstate->seed), tp, tn);
! 206: SIZ (rstate->seed) = tn;
! 207:
! 208: if (m2exp != 0)
! 209: {
! 210: /* Discard the lower half of the result. */
! 211: unsigned long int discardb = m2exp / 2;
! 212: mp_size_t discardl = discardb / BITS_PER_MP_LIMB;
! 213:
! 214: tn -= discardl;
! 215: if (tn > 0)
! 216: {
! 217: if (discardb % BITS_PER_MP_LIMB != 0)
! 218: {
! 219: mpn_rshift (tp, tp + discardl, tn, discardb % BITS_PER_MP_LIMB);
! 220: MPN_COPY (rp, tp, (discardb + BITS_PER_MP_LIMB -1) / BITS_PER_MP_LIMB);
! 221: }
! 222: else /* Even limb boundary. */
! 223: MPN_COPY_INCR (rp, tp + discardl, tn);
! 224: }
! 225: }
! 226: else
! 227: {
! 228: MPN_COPY (rp, tp, tn);
! 229: }
! 230:
! 231: TMP_FREE (mark);
! 232:
! 233: /* Return number of valid bits in the result. */
! 234: if (m2exp != 0)
! 235: retval = (m2exp + 1) / 2;
! 236: else
! 237: retval = SIZ (rstate->algdata.lc->m) * BITS_PER_MP_LIMB - shiftcount;
! 238: return retval;
! 239: }
! 240:
! 241: #ifdef RAWRANDEBUG
! 242: /* Set even bits to EVENBITS and odd bits to ! EVENBITS in RP.
! 243: Number of bits is m2exp in state. */
! 244: /* FIXME: Remove. */
! 245: unsigned long int
! 246: lc_test (mp_ptr rp, gmp_randstate_t s, const int evenbits)
! 247: {
! 248: unsigned long int rn, nbits;
! 249: int f;
! 250:
! 251: nbits = s->algdata.lc->m2exp / 2;
! 252: rn = nbits / BITS_PER_MP_LIMB + (nbits % BITS_PER_MP_LIMB != 0);
! 253: MPN_ZERO (rp, rn);
! 254:
! 255: for (f = 0; f < nbits; f++)
! 256: {
! 257: mpn_lshift (rp, rp, rn, 1);
! 258: if (f % 2 == ! evenbits)
! 259: rp[0] += 1;
! 260: }
! 261:
! 262: return nbits;
! 263: }
! 264: #endif /* RAWRANDEBUG */
! 265:
! 266: void
! 267: #if __STDC__
! 268: _gmp_rand (mp_ptr rp, gmp_randstate_t rstate, unsigned long int nbits)
! 269: #else
! 270: _gmp_rand (rp, rstate, nbits)
! 271: mp_ptr rp;
! 272: gmp_randstate_t rstate;
! 273: unsigned long int nbits;
! 274: #endif
! 275: {
! 276: mp_size_t rn; /* Size of R. */
! 277:
! 278: rn = (nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
! 279:
! 280: switch (rstate->alg)
! 281: {
! 282: case GMP_RAND_ALG_LC:
! 283: {
! 284: unsigned long int rbitpos;
! 285: int chunk_nbits;
! 286: mp_ptr tp;
! 287: mp_size_t tn;
! 288: TMP_DECL (lcmark);
! 289:
! 290: TMP_MARK (lcmark);
! 291:
! 292: chunk_nbits = rstate->algdata.lc->m2exp / 2;
! 293: tn = (chunk_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
! 294:
! 295: tp = (mp_ptr) TMP_ALLOC (tn * BYTES_PER_MP_LIMB);
! 296:
! 297: rbitpos = 0;
! 298: while (rbitpos + chunk_nbits <= nbits)
! 299: {
! 300: mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
! 301:
! 302: if (rbitpos % BITS_PER_MP_LIMB != 0)
! 303: {
! 304: mp_limb_t savelimb, rcy;
! 305: /* Target of of new chunk is not bit aligned. Use temp space
! 306: and align things by shifting it up. */
! 307: lc (tp, rstate);
! 308: savelimb = r2p[0];
! 309: rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
! 310: r2p[0] |= savelimb;
! 311: /* bogus */ if ((chunk_nbits % BITS_PER_MP_LIMB + rbitpos % BITS_PER_MP_LIMB)
! 312: > BITS_PER_MP_LIMB)
! 313: r2p[tn] = rcy;
! 314: }
! 315: else
! 316: {
! 317: /* Target of of new chunk is bit aligned. Let `lc' put bits
! 318: directly into our target variable. */
! 319: lc (r2p, rstate);
! 320: }
! 321: rbitpos += chunk_nbits;
! 322: }
! 323:
! 324: /* Handle last [0..chunk_nbits) bits. */
! 325: if (rbitpos != nbits)
! 326: {
! 327: mp_ptr r2p = rp + rbitpos / BITS_PER_MP_LIMB;
! 328: int last_nbits = nbits - rbitpos;
! 329: tn = (last_nbits + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB;
! 330: lc (tp, rstate);
! 331: if (rbitpos % BITS_PER_MP_LIMB != 0)
! 332: {
! 333: mp_limb_t savelimb, rcy;
! 334: /* Target of of new chunk is not bit aligned. Use temp space
! 335: and align things by shifting it up. */
! 336: savelimb = r2p[0];
! 337: rcy = mpn_lshift (r2p, tp, tn, rbitpos % BITS_PER_MP_LIMB);
! 338: r2p[0] |= savelimb;
! 339: if (rbitpos + tn * BITS_PER_MP_LIMB - rbitpos % BITS_PER_MP_LIMB < nbits)
! 340: r2p[tn] = rcy;
! 341: }
! 342: else
! 343: {
! 344: MPN_COPY (r2p, tp, tn);
! 345: }
! 346: /* Mask off top bits if needed. */
! 347: if (nbits % BITS_PER_MP_LIMB != 0)
! 348: rp[nbits / BITS_PER_MP_LIMB]
! 349: &= ~ ((~(mp_limb_t) 0) << nbits % BITS_PER_MP_LIMB);
! 350: }
! 351:
! 352: TMP_FREE (lcmark);
! 353: break;
! 354: }
! 355:
! 356: default:
! 357: gmp_errno |= GMP_ERROR_UNSUPPORTED_ARGUMENT;
! 358: break;
! 359: }
! 360: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to randraw.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp