Problems with Octave random number generation, part 1

classic Classic list List threaded Threaded
7 messages Options
Reply | Threaded
Open this post in threaded view
|

Problems with Octave random number generation, part 1

Rik-4
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





Reply | Threaded
Open this post in threaded view
|

Re: Problems with Octave random number generation, part 1

Juan Pablo Carbajal-2
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 <[hidden email]> 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
>
>
>
>
>

Reply | Threaded
Open this post in threaded view
|

Re: Problems with Octave random number generation, part 1

michaell
In reply to this post by Rik-4
First of all, let me note that whenever I do serious simulation programming
in C (to finally wind up in peer-reviewed publications, that is, as serious
as it gets), I always take directly the integer output from a PRNG and
transform it to the distribution I need, even if it is only a uniform
floating-point valued one. The reason is that I do not want to rely on such
issues being implemented correctly. Octave's rand() is in this sense just a
convenience function to be used when I quickly try something out. Thus, for
me non-uniformities on the level of 2^-32 are not essential, but of course,
we should try to be correct as long as it is reasonably doable.


Rik-4 wrote
>  static float randu32 (void)
>   {
>     return (static_cast
> <float>
>  (randi32 ()) + 0.5) * (1.0/4294967296.0);
>     /* divided by 2^32 */
>   }

Yes, that's a bug.


Rik-4 wrote
>   do
>     {
>       i = randi32 ();
>     }
>   while (i == 0 || i > (2^32 - 129);
>
>   return i * (1.0 / 4294967296.0);

No, if that is really what the patch boils down to, that's not a problem in
my eyes.

Both Philipp's and your proposals give uniformly distributed random numbers
as long as one looks at them on a scale that is coarser than eps. The
difference between your and Philipp's proposal is that his uses all the
representable floating-point numbers (at least all above 2^-8, probably --
if you used a randi32(), you would be even finer at small values), and as
the discretization becomes finer at small values, a given number becomes
less probable. But that is okay, because if you look into how often you hit
a given interval, it comes out correctly. Actually I would prefer his
proposal, _specifically_ if you want to pass it through log() to generate
exponentially-distributed numbers (but see below).

By the way, you definitely should drop the condition (i == 0xFFFFFF) in your
code, it will make you miss the last float value below 1.

But I am still not quite convinced by either of your proposals, because both
of you introduce some bias (and with 24-bit mantissas, this will be easily
visible): let's assume randi32() is quite dumb and has a cycle length of
2^32, during which time it returns every uint32 exactly once. Let's consider
the histogram of floats after one such cycle resulting from your approaches.
In the region around, say, 0.6, both of your proposals give the correct
density of 2^32 per unit interval. However, in Philipp's proposal 2^32-128
values have been distributed over the unit interval, while in your case it
were only 2^32-256 (after you have dropped the (i == 0xFFFFFF) condition),
because you skip on the zeros. So somewhere these numbers have to be
missing. In your case it is clear, they are missing from the endpoints of
the interval (symmetrically). In Philipp's case, they are missing at 0.5,
0.25... and so on, every time the exponent decreases by one (because any
number above 0.5 gets a range of plus and minus eps/2 around itself, but the
interval that is assigned to exact 0.5 extends only eps/4 to the left).

I think that the best solution would be to a simple

return (randi32 () + 0.5) * (1.0 / 4294967296.0);

but use one of the non-default IEEE-754 rounding-modes "round towards minus
infinity" or "round towards zero" (and verify that casting doubles to floats
honours the rounding mode, which I do believe). The vast majority of
computers that today run octave will have no problem with that, but is there
a relevant architecture where this does not work?


Rik-4 wrote
> #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);

If this is really true, this is a massive problem, but fortunately with an
easy solution: isn't there a randu64()? If not, one should construct a
double random number with all mantissa bits filled from two calls to
randu32(), and at the very least pass that through log. The point is that
2^32 is not so large, you can have computers where you can fill the RAM in a
few seconds by 2^32 exponentially distributed random numbers. With your
proposal, Rik, the largest number you could ever get is 16.636, with
Philipp's it is 22.181, but actually you would expect that the maximum of N
independent exponentially-distributed random numbers goes with log(N). Thus,
already at 2^24 or 2^32 numbers you would start to see significant
deviations in the statistics of the maximum. If you use the full width of
double, you need 2^53 numbers for that, which is not within the reach of
computers on which octave runs today.




--
Sent from: https://octave.1599824.n4.nabble.com/Octave-Maintainers-f1638794.html

Reply | Threaded
Open this post in threaded view
|

Re: Problems with Octave random number generation, part 1

michaell
In reply to this post by Juan Pablo Carbajal-2
That's in essence Rik's approach, but not rejecting 0. Thus, it would be
incompatible with Matlab, which specifies (0,1).



--
Sent from: https://octave.1599824.n4.nabble.com/Octave-Maintainers-f1638794.html

Reply | Threaded
Open this post in threaded view
|

Re: Problems with Octave random number generation, part 1

Daniel Sebald
In reply to this post by Rik-4
On 12/17/19 2:27 PM, Rik 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);
> is
>    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.(1.0/4294967296.0)

In the discussion between Jordi, Mike and Philipp is the link to the code:

http://hg.savannah.gnu.org/hgweb/octave/file/c579bd4e12c9/liboctave/numeric/randmtzig.c#l392

Yes, I don't like the notion of tossing out a range of numbers using the
while loop.  RNG's have been extensively analyzed for randomness, and
the conditional statement above has the *potential* to effect that.
Then again, one analysis might be that good pseudo-RNG should have
subintervals that look properly random.

Small technicality is that

   return ((float)randi32 () + 0.5) * (1.0/4294967296.0);
   /* divided by 2^32 */

is actually multiplying by the inverse, not dividing.  Fine for real
numbers and a CPU-saving technique in a general sense, but there will
always be a bit of a guess as to exactly what a hardware floating point
processor is doing.  It probably doesn't make a difference in this case,
though.

But along that line, is there a slightly smarter way of doing the
arithmetic closer to integer math?  For example, if there is a way of
avoiding the rounding problems by generating a range that is as

delta + [0,M-1] * alpha

it might avoid the rare value going out of range (which is kind of what
you did in your proposed alternative later).


> 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  return ((float)randi32 () + 0.5) * (1.0/4294967296.0);

   /* divided by 2^32 */

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

Creative, but mapping out that whole upper level of pseudo-random
numbers doesn't feel good.  It may be as good as anything though
because, yeah, we could change the name of the function to randu24() and
it would have the same meaning because of the imposed range of (0,1).
In some way, that range must implicitly toss away the exponent and the
floating point bit.  It's effectively now fractional point, not floating
point.

Are we doing a "horse-is-out-of-the-barn" type of thing by generating
this (0:1) floating point range and then expecting to scale it?  Say I
take a vector of numbers from the RNG with range [1/2^24 : 2^24-1/2^24]
(or think of mantissa [1:0xFFFFFE] if you like) and multiply by floating
point value greater than 1?  How about less than 1?  Probably some
strange dithering effect on the floating point numbers in either case,
unless the scaling is a power of 2.


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

This is probably where the significance of [0:1) versus (0:1) makes a
difference.  Several math functions throw an error or require special
treatment with an input of 0.  It would be a mistake to think of
floating point in terms of real number (0:1).  Instead it means
[delta:1-delta] where delta is the resolution of the floating point number.

Dan

Reply | Threaded
Open this post in threaded view
|

Re: Problems with Octave random number generation, part 1

Gene Harvey
I think we can improve on Rik's solution and lose the loop. Consider
====================================================
static float randu32 (void)
{
  return ((randi32 () & uint32_t (0x00FFFFFE)) + 1.f)
           / std::exp2 (std::numeric_limits<float>::digits);
}
====================================================
which has a range of uniform (0, 1). However, this still only allows for
2^24 - 1 unique values in the interval, but it seems quite simple. I also
tried reworking Philipp's idea and ended up with
====================================================
static float randu32 (void)
{
  return std::nextafterf (randi32 () + 0.5f, 0)
           / std::numeric_limits<uint32_t>::max ();
}
====================================================
which ends up with 79,691,776 unique values. (that's 0x04C00000 if
you'd like to figure out the math behind it; I just ran it through a loop).

Hope these help.

On Tue, Dec 17, 2019 at 11:40 PM Daniel J Sebald <[hidden email]> wrote:

>
> On 12/17/19 2:27 PM, Rik 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);
> > is
> >    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.(1.0/4294967296.0)
>
> In the discussion between Jordi, Mike and Philipp is the link to the code:
>
> http://hg.savannah.gnu.org/hgweb/octave/file/c579bd4e12c9/liboctave/numeric/randmtzig.c#l392
>
> Yes, I don't like the notion of tossing out a range of numbers using the
> while loop.  RNG's have been extensively analyzed for randomness, and
> the conditional statement above has the *potential* to effect that.
> Then again, one analysis might be that good pseudo-RNG should have
> subintervals that look properly random.
>
> Small technicality is that
>
>    return ((float)randi32 () + 0.5) * (1.0/4294967296.0);
>    /* divided by 2^32 */
>
> is actually multiplying by the inverse, not dividing.  Fine for real
> numbers and a CPU-saving technique in a general sense, but there will
> always be a bit of a guess as to exactly what a hardware floating point
> processor is doing.  It probably doesn't make a difference in this case,
> though.
>
> But along that line, is there a slightly smarter way of doing the
> arithmetic closer to integer math?  For example, if there is a way of
> avoiding the rounding problems by generating a range that is as
>
> delta + [0,M-1] * alpha
>
> it might avoid the rare value going out of range (which is kind of what
> you did in your proposed alternative later).
>
>
> > 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  return ((float)randi32 () + 0.5) * (1.0/4294967296.0);
>
>    /* divided by 2^32 */
> > 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?
>
> Creative, but mapping out that whole upper level of pseudo-random
> numbers doesn't feel good.  It may be as good as anything though
> because, yeah, we could change the name of the function to randu24() and
> it would have the same meaning because of the imposed range of (0,1).
> In some way, that range must implicitly toss away the exponent and the
> floating point bit.  It's effectively now fractional point, not floating
> point.
>
> Are we doing a "horse-is-out-of-the-barn" type of thing by generating
> this (0:1) floating point range and then expecting to scale it?  Say I
> take a vector of numbers from the RNG with range [1/2^24 : 2^24-1/2^24]
> (or think of mantissa [1:0xFFFFFE] if you like) and multiply by floating
> point value greater than 1?  How about less than 1?  Probably some
> strange dithering effect on the floating point numbers in either case,
> unless the scaling is a power of 2.
>
>
> > 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.
>
> This is probably where the significance of [0:1) versus (0:1) makes a
> difference.  Several math functions throw an error or require special
> treatment with an input of 0.  It would be a mistake to think of
> floating point in terms of real number (0:1).  Instead it means
> [delta:1-delta] where delta is the resolution of the floating point number.
>
> Dan
>

Reply | Threaded
Open this post in threaded view
|

Re: Problems with Octave random number generation, part 1

Daniel Sebald
In reply to this post by Daniel Sebald
On 12/18/19 12:39 AM, Daniel J Sebald wrote:
> On 12/17/19 2:27 PM, Rik wrote:
>> Maintainers,
[snip]

>> 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?
>
> Creative, but mapping out that whole upper level of pseudo-random
> numbers doesn't feel good.  It may be as good as anything though
> because, yeah, we could change the name of the function to randu24() and
> it would have the same meaning because of the imposed range of (0,1). In
> some way, that range must implicitly toss away the exponent and the
> floating point bit.  It's effectively now fractional point, not floating
> point.

In fact, the range [1:0xFFFFFE] could go directly into the mantissa
rather than multiply or divide.  I.e., think of 1/16777216 as 2^-24 or
right bit shift by 24.  Is there a C macro for loading a mantissa directly?

Dan