# Random integer generation

Kevin Ryde user42@zip.com.au
Fri, 21 Feb 2003 09:44:28 +1000

--=-=-=

"Sisyphus" <kalinabears@hdc.com.au> writes:
>
> It was suggested to me (off list) that I read up on on linear congruential
> generators in Knuth.

Pedro Gimeno points out that it's not necessarily silly to have a
smallish "a", a poor choice of words on my part.

I had a go at improving the code, it should fix the case of small or
zero "a", with a bit of luck.

--=-=-=
Content-Disposition: attachment; filename=randraw.c.small-a.diff

--- randraw.c.old	2003-02-21 08:50:32.000000000 +1000
+++ randraw.c	2003-02-21 09:41:04.000000000 +1000
@@ -109,7 +109,9 @@
seedp = PTR (rstate->_mp_seed);
seedn = SIZ (rstate->_mp_seed);

-  if (seedn == 0)
+  an = SIZ (rstate->_mp_algdata._mp_lc->_mp_a);
+
+  if (seedn == 0 || an == 0)
{
/* Seed is 0.  Result is C % M.  Assume table is sensibly stored,
with C smaller than M*/
@@ -121,13 +123,12 @@
}

ap = PTR (rstate->_mp_algdata._mp_lc->_mp_a);
-  an = SIZ (rstate->_mp_algdata._mp_lc->_mp_a);

/* Allocate temporary storage.  Let there be room for calculation of
(A * seed + C) % M, or M if bigger than that.  */

TMP_MARK (mark);
-  ta = an + seedn + 1;
+  ta = MAX (m2exp/GMP_NUMB_BITS, an + seedn) + 1;
tp = (mp_ptr) TMP_ALLOC (ta * BYTES_PER_MP_LIMB);

/* t = a * seed */
@@ -137,12 +138,13 @@
mpn_mul (tp, ap, an, seedp, seedn);
tn = an + seedn;

+  if (tn < m2exp/GMP_NUMB_BITS)
+    MPN_ZERO (tp+tn, m2exp/GMP_NUMB_BITS - tn);
+
/* t = t + c */
tp[tn] = 0;			/* sentinel, stops MPN_INCR_U */
MPN_INCR_U (tp, tn, c);

-  ASSERT_ALWAYS (m2exp / GMP_NUMB_BITS < ta);
-
/* t = t % m */
tp[m2exp / GMP_NUMB_BITS] &= ((mp_limb_t) 1 << m2exp % GMP_NUMB_BITS) - 1;
tn = (m2exp + GMP_NUMB_BITS - 1) / GMP_NUMB_BITS;

--=-=-=--