Re: Random numbers (again) broken on 64-bit platforms

[Andreas Rottmann]

Andy Wingo <address@hidden> writes:

> Heya,
>
> On Mon 26 Jul 2010 23:01, Andreas Rottmann <address@hidden> writes:
>
>> scheme@(guile-user)> (random (ash 1 32))
>> ERROR: In procedure random:
>> ERROR: Argument 1 out of range: 4294967296
>
> Well, it's not a segfault at least :)
>
> The point of that change was to indicate that the RNG did not return
> sizeof(unsigned long) bits of randomness; it always returned 32
> bits. See the note in "BSD Random Number Functions" in libc's manual:
>
>      *NB:* Temporarily this function was defined to return a `int32_t'
>      value to indicate that the return value always contains 32 bits
>      even if `long int' is wider.  The standard demands it differently.
>      Users must always be aware of the 32-bit limitation, though.
>
> I'll fix this now to avoid the error, but there is still work to do on
> the RNG -- we really need to update the RNG, I think. Brian Gough, the
> GSL maintainer, says MT19937 is the one to use, and specifically the new
> SIMD version at
> http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/index.html. We
> should replace our MWC RNG with that one.
>
I've spotted an issue with your fix: the way 64-bit random numbers are
generated in scm_random() is bogus; to test that, I used the following
code:

(define (random-max n iterations)
  (let loop ((i 0) (result 0))
    (if (< i iterations)
        (loop (+ i 1) (max (random n) result))
        result)))

(random-max #x1ffffffff 10000000)

If you try this, you'll notice that `random-max' returns a much lower
number than it should. The attached patch hopefully fixes this issue

From: Andreas Rottmann <address@hidden>
Subject: Fix the range of `random' on 64-bit platforms

For > 32 bit integers still in the fixnum range, scm_random() would
return random numbers with a lower range than specified.

* libguile/random.c (scm_i_mask32): New static inline function.
  (scm_c_random): Use `scm_i_mask32'.
  (scm_c_random64): New function, 64-bit variant of scm_c_random.
  (scm_random): Use `scm_c_random64' instead of forming the 64-bit random 
  number in a bogus way.
* libguile/random.h: Added `scm_c_random64'.

---
 libguile/random.c |   42 +++++++++++++++++++++++++++---------------
 libguile/random.h |    1 +
 2 files changed, 28 insertions(+), 15 deletions(-)

diff --git a/libguile/random.c b/libguile/random.c
index 83870f6..c0a3e04 100644
--- a/libguile/random.c
+++ b/libguile/random.c
@@ -241,21 +241,42 @@ scm_c_exp1 (scm_t_rstate *state)
 
 unsigned char scm_masktab[256];
 
-scm_t_uint32
-scm_c_random (scm_t_rstate *state, scm_t_uint32 m)
+static inline scm_t_uint32
+scm_i_mask32 (scm_t_uint32 m)
 {
-  scm_t_uint32 r, mask;
-  mask = (m < 0x100
+  return (m < 0x100
          ? scm_masktab[m]
          : (m < 0x10000
             ? scm_masktab[m >> 8] << 8 | 0xff
             : (m < 0x1000000
                ? scm_masktab[m >> 16] << 16 | 0xffff
                : scm_masktab[m >> 24] << 24 | 0xffffff)));
+}
+
+scm_t_uint32
+scm_c_random (scm_t_rstate *state, scm_t_uint32 m)
+{
+  scm_t_uint32 r, mask = scm_i_mask32 (m);
   while ((r = state->rng->random_bits (state) & mask) >= m);
   return r;
 }
 
+scm_t_uint64
+scm_c_random64 (scm_t_rstate *state, scm_t_uint64 m)
+{
+  scm_t_uint64 r;
+  scm_t_uint32 mask;
+
+  if (m <= SCM_T_UINT32_MAX)
+    return scm_c_random (state, (scm_t_uint32) m);
+  
+  mask = scm_i_mask32 (m >> 32);
+  while ((r = ((scm_t_uint64) (state->rng->random_bits (state) & mask) << 32)
+          | state->rng->random_bits (state)) >= m)
+    ;
+  return r;
+}
+
 /*
   SCM scm_c_random_bignum (scm_t_rstate *state, SCM m)
 
@@ -377,17 +398,8 @@ SCM_DEFINE (scm_random, "random", 1, 1, 0,
       return scm_from_uint32 (scm_c_random (SCM_RSTATE (state),
                                             (scm_t_uint32) m));
 #elif SCM_SIZEOF_UNSIGNED_LONG <= 8
-      if (m <= SCM_T_UINT32_MAX)
-        return scm_from_uint32 (scm_c_random (SCM_RSTATE (state),
-                                              (scm_t_uint32) m));
-      else
-        {
-          scm_t_uint64 upper, lower;
-          
-          upper = scm_c_random (SCM_RSTATE (state), (scm_t_uint32) (m >> 32));
-          lower = scm_c_random (SCM_RSTATE (state), SCM_T_UINT32_MAX);
-          return scm_from_uint64 ((upper << 32) | lower);
-        }
+      return scm_from_uint64 (scm_c_random64 (SCM_RSTATE (state),
+                                              (scm_t_uint64) m));
 #else
 #error "Cannot deal with this platform's unsigned long size"
 #endif
diff --git a/libguile/random.h b/libguile/random.h
index 512b7d2..2f1f0f6 100644
--- a/libguile/random.h
+++ b/libguile/random.h
@@ -67,6 +67,7 @@ SCM_API double scm_c_uniform01 (scm_t_rstate *);
 SCM_API double scm_c_normal01 (scm_t_rstate *);
 SCM_API double scm_c_exp1 (scm_t_rstate *);
 SCM_API scm_t_uint32 scm_c_random (scm_t_rstate *, scm_t_uint32 m);
+SCM_API scm_t_uint64 scm_c_random64 (scm_t_rstate *state, scm_t_uint64 m);
 SCM_API SCM scm_c_random_bignum (scm_t_rstate *, SCM m);
 
 
-- 
tg: (7387c23..) t/random64-fix-2 (depends on: master)

> Would you be interested in doing this? We would need some test
> suites too, I think, and possibly changes to the scm_t_rng structure.
>
Sorry, I don't have the inclination to work on this ATM. Also, random
number tests are kinda hard to write -- it's random stuff, after all
:-).

Regards, Rotty
-- 
Andreas Rottmann -- <http://rotty.yi.org/>