Sat 01 Mar 2014 12:09:31 PM UTC, original submission:
According to the documentation rand returns a number between (0, 1) and it is conventional practice to be able to map the returned number to an integer between 1 and a using an expression such as:
x = round (rand * a + 0.5)
However, on rare occasions rand actually returns exactly 1.0 which causes the expression to return a+1 instead. I've reviewed the code with Jordi Gutierrez Hermoso, and I think the problem may be in lines 396 and 404-406 of the file http://hg.savannah.gnu.org/hgweb/octave/file/c579bd4e12c9/liboctave/numeric/randmtzig.c#l400
I am suspicious of the float cast and the addition of a constant. Jordi has other ideas which I will quote below, since I don't fully understand his interpretation on the loss of precision issue.
- Mike
---
So I think I see the bug:
http://hg.savannah.gnu.org/hgweb/octave/file/c579bd4e12c9/liboctave/numeric/randmtzig.c#l392
Here the cast to float is pointless. The constants are doubles, not
floats, so the computation is done in doubles, with 32 bits. However,
a float mantissa only has 23 bits, so casting to float loses precision
bits, and I think this casting also performs rounding.
I think it might be better to do bit manipulations to light the
appropriate mantissa bits in the float representation than to be
relying on implicit casts, explicit casts, and the arithmetic between
them.
- Jordi G. H.
|