[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: Random numbers (again) broken on 64-bit platforms
From: |
Andreas Rottmann |
Subject: |
Re: Random numbers (again) broken on 64-bit platforms |
Date: |
Sun, 01 Aug 2010 20:44:42 +0200 |
User-agent: |
Gnus/5.13 (Gnus v5.13) Emacs/23.2 (gnu/linux) |
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/>
- Re: Random numbers (again) broken on 64-bit platforms,
Andreas Rottmann <=