Re: Random number generation quirk

[Jordi Gutiérrez Hermoso]

Moving this to the dev list...

On Fri, 2014-02-28 at 12:14 -0500, Jordi Gutiérrez Hermoso wrote:
> Ah, you're using single. I can reproduce now on 3.8.0:
> 
>     octave:1> x = rand(3000,"single"); max(round(x*10+0.5)(:))
>     ans =  10
>     octave:2> x = rand(3000,"single"); max(round(x*10+0.5)(:))
>     ans =  11
>     octave:3> max(x(:))
>     ans =  1
> 
> This a bug one way or another, since the documentation clearly
> states that this should not happen. Can you please file a bug
> report?

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.