Random integer generation

[Kevin Ryde]

--=-=-=

"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;

--=-=-=--