Random numbers (again) broken on 64-bit platforms

Random numbers (again) broken on 64-bit platforms

[Andreas Rottmann]

Hi!

After pulling from Git today, I noticed that `random' is still/again
broken on 64-bit platforms:

scheme@(guile-user)> (random (ash 1 32))
ERROR: In procedure random:
ERROR: Argument 1 out of range: 4294967296

This looks like being caused by
b606ff6af9c9ec7fc3c4473c09ce1e95c18f024a, which doesn't take into
account (in scm_random) that an INUM may well exceed 32 bits on 64-bit
platforms. I'll try to come up with a patch for that, unless someone
beats me to it.

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

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

[Andy Wingo]

Heya,

On Mon 26 Jul 2010 23:01, Andreas Rottmann <a.rottmann@gmx.at> 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.

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.

Cheers,

Andy
-- 
http://wingolog.org/

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

[Andreas Rottmann]

[-- Attachment #1: Type: text/plain, Size: 1622 bytes --]

Andy Wingo <wingo@pobox.com> writes:

> Heya,
>
> On Mon 26 Jul 2010 23:01, Andreas Rottmann <a.rottmann@gmx.at> 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


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #2: random64-fix-2.diff --]
[-- Type: text/x-diff, Size: 3413 bytes --]

From: Andreas Rottmann <a.rottmann@gmx.at>
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);
 
 \f
-- 
tg: (7387c23..) t/random64-fix-2 (depends on: master)

[-- Attachment #3: Type: text/plain, Size: 347 bytes --]


> 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/>

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

[Andy Wingo]

Hi,

On Sun 01 Aug 2010 20:44, Andreas Rottmann <a.rottmann@gmx.at> writes:

> I've spotted an issue with your fix: the way 64-bit random numbers are
> generated in scm_random() is bogus;

Grr, indeed. Fix applied, thanks.

>> Would you be interested in [importing another rng]? 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
> :-).

Too bad, we could really use it, and now is the time for any needed ABI
breaks. Thanks anyway for the fix!

Andy
-- 
http://wingolog.org/

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

[No Itisnt]

On Sun, Aug 1, 2010 at 2:57 PM, Andy Wingo <wingo@pobox.com> wrote:
> Hi,
>
> On Sun 01 Aug 2010 20:44, Andreas Rottmann <a.rottmann@gmx.at> writes:
>
>> I've spotted an issue with your fix: the way 64-bit random numbers are
>> generated in scm_random() is bogus;
>
> Grr, indeed. Fix applied, thanks.
>
>>> Would you be interested in [importing another rng]? 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
>> :-).
>
> Too bad, we could really use it, and now is the time for any needed ABI
> breaks. Thanks anyway for the fix!

I may have time for doing this later on but definitely not this month;
do you guys have any sort of time frame for a release?

And also: is there any way to set the lower bound of a random call?

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

[Andy Wingo]

Hi,

On Mon 02 Aug 2010 00:06, No Itisnt <theseaisinhere@gmail.com> writes:

> I may have time for doing [rng work] later on but definitely not this
> month; do you guys have any sort of time frame for a release?

Very cool!

We plan to release a 1.9.12 within a week or two, then 2.0 should be the
next one. I know we've been slipping with our self-imposed deadlines,
but it does feel like we're getting closer.

> And also: is there any way to set the lower bound of a random call?

I don't think so, no. One can of course call random with the range, and
add an offset.

Cheers,

Andy
-- 
http://wingolog.org/