[Octave-bug-tracker] [bug #54619] randi() is biased

[Michael Leitner]

Follow-up Comment #28, bug #54619 (project octave):

I cannot remember why I thought this flintmax()-1 was needed, and now I am
actually quite convinced that it gives incorrect results, as the line 


r_prim = floor ( rand (M,1) * (flintmax ()-1) );


would lead to 0 appearing in r_prim with double probability. No, using just
flintmax() everywhere instead of flintmax()+1 should be perfectly okay -- the
results should stay perfectly the same, but the reader of the code would not
be confused. 

By the way, I will document what I found about octave's rand() to answer my
questions in comment #19: according to https://hal.science/hal-02427338,
octave's rand() works by first generating a uniformly distributed integer
between 0 and 2^53-1, adds 0.4 to it (looks like an ad hoc-decision, as no
other popular software package does it that way) and divides by 2^53. Thus, it
obviates (partly?) the problem that the less-significant bits of the returned
doubles below 2^-11 are all zero, at the price of a very small bias towards
lower values. 

In the context of randi(), everything should be perfect: multiplication by
flintmax() undoes the division by 2^53, floor undoes the addition of 0.4, so
that r_prim should really be a uniformly distributed integer. Thus, apart from
the potential issue of inexact floating-point division as mentioned in comment
#19, as far as I can tell randi() should work as advertised.


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?54619>

_______________________________________________
Message sent via Savannah
https://savannah.gnu.org/