Re: Random number generation quirk

[Michael Pender]

I have had the problem occur with double rands, but it is much easier to reproduce with singles.

I think loss of precision could account for the problem.  I was under the impression that the division was taking place as a single precision calculation, not a double precision calculation?

I would be more confident in the correctness of the result if all operands were doubles, and the division took place as a double precision division, and then the final result was cast as a single.

That said, it would probably be much faster to do bit manipulations place the bits in the appropriate locations.  However, I am not sure whether assumptions about x86 versus other architectures will require specialized coding.  I don't completely understand the comments from lines 104-108 on performance impacts of the i686 64-bit architecture.

- Mike

On Fri, Feb 28, 2014 at 1:39 PM, Jordi Gutiérrez Hermoso <address@hidden> wrote:

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.