Re: Problems with Octave random number generation, part 1

[Juan Pablo Carbajal]

Just for the sake of comparison here is the code used in numpy, the
convention (and seems to be widespread is to sample [0,1) )

https://github.com/numpy/numpy/blob/master/numpy/core/include/numpy/random/distributions.h
https://github.com/numpy/numpy/blob/master/numpy/random/src/distributions/distributions.c

https://github.com/numpy/numpy/blob/master/numpy/core/include/numpy/random/bitgen.h

There we see

static NPY_INLINE float next_float(bitgen_t *bitgen_state) {
return (next_uint32(bitgen_state) >> 9) * (1.0f / 8388608.0f);
}

and in the mt library we see

https://github.com/numpy/numpy/blob/master/numpy/random/src/mt19937/mt19937.h

static inline double mt19937_next_double(mt19937_state *state) {
int32_t a = mt19937_next(state) >> 5, b = mt19937_next(state) >> 6;
return (a * 67108864.0 + b) / 9007199254740992.0;
}


On Tue, Dec 17, 2019 at 8:27 PM Rik <address@hidden> wrote:
>
> Maintainers,
>
> I believe there are serious, but subtle, problems with Octave's random
> number generation.  Resolving them, without introducing further issues,
> requires careful thought and review beyond just a single coder (myself).
> Hence, I'm bringing this issue to the Maintainer's list.  I've also CC'ed
> Jordi and Philipp Kutin who were active in initially discovering the bugs,
> and Michael Leitner who was helpful in removing the bias from Octave's
> randi function.
>
> The tip of the iceberg was bug report #41742
> (https://savannah.gnu.org/bugs/?41742).  In this report, using rand to
> generate an integer in the range [1, a] using the straightforward code
>
> x = round (rand * a + 0.5)
>
> occasionally fails and returns a result of a+1 outside the desired range.
> This is only possible if rand should return exactly 1.  But, according to
> the documentation, rand produces a floating point number in the range (0,
> 1), i.e., endpoints are excluded.  Unfortunately, it is quite simple to
> demonstrate that the rand function does include the endpoint.  Test code
> that I found that always reproduces the problem is
>
> octave:1> rand ("state", 1);
> octave:2> x = max (rand ([20e6, 1], "single"))
> x = 1
>
> This result is contrary to the documentation, and contrary to the behavior
> of the rand function in Matlab.  Discussion between Jordi and Philipp
> ensued at this thread
> https://octave.1599824.n4.nabble.com/Re-Random-number-generation-quirk-td4662514.html.
>
> The Octave code which generates a "single" (32-bit IEEE floating point)
> random number is found in liboctave/numeric/randmtzig.cc
>
>   /* generates a random number on (0,1)-real-interval */
>   static float randu32 (void)
>   {
>     return (static_cast<float> (randi32 ()) + 0.5) * (1.0/4294967296.0);
>     /* divided by 2^32 */
>   }
>
> The function randi32 (which I have no reason to doubt) generates an
> unsigned 32-bit integer in the range [0, 2^32 - 1].  Unfortunately, I think
> there was some confusion by the original programmer between Octave and C
> implicit conversion rules.  In Octave, arguments are demoted to lower
> precision before the operand is applied:
>
> single + double = single
>
> This is exactly reversed in C where arguments are promoted to the more
> precise class.  In C,
>
> single + double = double
>
> So this code is actually performed completely with double precision values,
> and then truncated to single precision float at the end of the function
> because the return type for the function is float.  In Octave pseudo-code,
>
> single ( (double (randi32 ()) + 0.5) * ( 1.0 / 4294967296.0) )
>
> A useful Octave function for mimicking the C code is
>
> fs = @(x) single ((floor (x) + 0.5) * (1.0 / 4294967296.0))
>
> Using fs() one can determine the maximum valuable allowable in order for
> the result not to become exactly equal to 1.
>
> octave:28> fs (2^32-129)
> ans = 0.9999999
> octave:29> fs (2^32-128)
> ans = 1
>
> I assume Philipp did something similar.  He has provided a patch that,
> after unpacking various macros, does the following
>
> static float randu32 (void)
> {
>   uint32_t i;
>
>   do
>     {
>       i = randi32 ();
>     }
>   while (i == 0 || i > (2^32 - 129);
>
>   return i * (1.0 / 4294967296.0);
> }
>
> My fear, and more mathematical help in this specific would be nice, is that
> this will introduce a bias.  In Octave, the new code can be modeled as
>
> fs2 = @(x) single (floor (x) * (1.0 / 4294967296.0 ))
>
> Trying some sample values
>
> octave:78> fs2(1)
> ans = 2.3283064e-10
> octave:79> fs2(2)
> ans = 4.6566129e-10
> octave:80> fs2(3)
> ans = 6.9849193e-10
> octave:82> fs2(2^32-128)
> ans = 1
> octave:83> fs2(2^32-129)
> ans = 0.9999999
> octave:84> fs2(2^32-130)
> ans = 0.9999999
> octave:100> fs2(2^32-129) - fs2(2^32-383)
> ans = 0
>
> What one notices is that the small integer values all get mapped to
> different floating point values, but at the upper end, large ranges of
> different integer values get mapped to the same floating point value.
>
> Given that IEEE single-precision floating point has only 24 bits for the
> mantissa, I think something like this would work better
>
>   static float randu32 (void)
>   {
>     uint32_t i;
>
>     do
>       {
>         i = randi32 () & static_cast<uint32_t> (0xFFFFFF);
>       }
>     while (i == 0 || i == 0xFFFFFF);
>
>     return i * (1.0f / 16777216.0f);
>   }
>
> In my testing, this does not produce exact values of 1, nor does it map
> multiple integer values to the same floating point value.  Is anyone aware
> of problems with this algorithm?
>
> Resolving this must happen, but is only part of the issue.  There is a
> #define in randmtzig.cc
>
> #define RANDU randu32()
>
> and then RANDU is used in the algorithms for rand_normal and
> rnd_exponential.  In particular, the use occurs as the argument to the log
> function.  For example,
>
> yy = - std::log (RANDU);
>
> The old RANDU = randu32() occasionally put out values of 1 which will be
> mapped by the log function to 0.  Whereas any new function will always be
> slightly less than one which means the log function output will be small,
> but not exactly 0.  Do rand_exponential and rand_normal rely on the old
> behavior in order to function correctly?  I have no idea.
>
> --Rik
>
>
>
>
>