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

[Michael Leitner]

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

See the attached patch (a diff against the current randi.m). I think it is
correct, at least it passes the test on the ratio of evens and odds. This
could be also added in the tests, but I do not know how to write a test on
stochastic functions. Of course, that it does pass this specific test does not
prove anything, so I invite also other people to have a look on it.
Essentially it just tries to create integers between 0 and flintmax-1 (2^53-1
for IEEE-754 doubles), divides them by another integer K so that the requested
range is covered (for exact integer arithmetic this should result in exactly
uniform probabilities within the requested range) and rejects all above the
range. It initially generates a buffer of random numbers that is in all
likelihood larger than what is requested, so that it is only insignificantly
slower than the previous implementation for cases where not many numbers have
to be dropped (for small ranges, that is, where the previous implementation
was already nearly correct). In these cases, the sequence of returned numbers
is also the same, provided not too many numbers are requested. 

A potential issue: is IEEE-754 double division predicted to be exact on
exactly representable integers, when the result is again an integer? If
division is multiplication by the inverse, then probably not, because the
inverse is in general not exactly representable. In this case, it would
probably be better to do these manipulations in uint64, which should be
required to be exact. 

As a side note: I would be very much in favour of exposing the fundamental
random number generator not via rand(), which returns doubles in the range
[0,1], but by a function returning uniformly distributed uint64, which they
internally indeed are. This would save operations in the present problem which
have only the purpose to undo this conversion to scaled doubles, further you
can get the scaled doubles by a two-line m-file, and finally this would then
document what rand() _really_ returns (is exact 0 and/or exact 1 included, if
so, with which probability, and is indeed every double between 0 and 1
returned or only those that are rounded uint64 divided by 2^64 and so on).

(file #44961)
    _______________________________________________________

Additional Item Attachment:

File name: bug54619_patch                 Size:0 KB


    _______________________________________________________

Reply to this item at:

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

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