fix rand(Float16|Float32) equal to 1

[rfourquet]

Cf. bug reported by @tkelman in #9301.
As noted by @ivarne, converting a Float64 different from 1 to Float16
or Float32 can produce a value equal to one, which doesn't agree with
the documented API.
So instead of converting, let's produce a Float16/Float32 in [1,2) at
the bit level, and then substract 1 (as was done in the array versions).
As a side effect, this also makes rand(Float16) about twice as fast.

[ivarne]

Why not use

rand(r::AbstractRNG, ::Type{Float16}) =
    Float16(reinterpret(Float32, rand(r, UInt16) % UInt32 & 0x007fe000 | 0x3f800000) - 1)

rand(r::AbstractRNG, ::Type{Float32}) =
    reinterpret(Float32, rand(r, UInt32) % UInt32 & 0x007fffff | 0x3f800000) - 1

so that all other AbstracRNG implementations becomes correct for Float32 and Float64 automatically?

Be aware that I don't understand what rand_ui52_raw means, so maybe we should keep the speical version for MersenneTwister as well.

[rfourquet]

The reason was that if a particular RNG defines e.g. rand{T<:Union(Float16,Float32,Float64)}(r::ThisRNG, ::Type{T}) = ..., there will be ambiguity warning with your solution, forcing the implementor of ThisRNG to loop over the types and use @eval or some other workaround (that was my mistake in #9274 with rand(r::AbstractRNG, ::Type{Close1Open2}) = rand(r) + 1.0, which I removed later). But you are totally right that the problem to make this implementation available to other RNGs is open. I would lean towards having something like

stagedfunction{T} rand(r::AbstractRNG, ::Type{T})
    if T === Float16 and implements_UInt32(r):
        :(reinterpret(Float32, rand($r, UInt32) % UInt32 & 0x007fffff | 0x3f800000) - 1))
    elseif ...
...
end

That way, the general fallback will be chosen only if not other implementation at all is available. Another solution would be for an RNG to call a like macro @generate_rand_Float32_from_UInt32.

Concerning rand_ui52_raw: it returns an UInt64 with at least 52 random bits, raw meaning that the 12 higher bits are unspecified (may not be zero). So here I define rand_ui23_raw because only 23 random bits are needed, and each RNG will have its own best way to produce an unsigned with such 23 bits.

[ivarne]

If we don't figure out a way to implement

rand{T<:Union(Float16,Float32,Float64)}(r::AbstractRNG, ::Type{T})

correctly, I can't see why other implementations would want to do that. Either there is a way to write that generic implementation, or there isn't.

For rand_ui**_raw; could we at least relax the definition for r::RandomDevice to r::AbstractRNG. That definition should be correct, if not optimal, for all AbstractRNG, right?

PS: I hope you enjoy my comments, and don't get frustrated with my lack of understanding. I have big interest in random numbers and extensible generic interfaces, and believe that any code review is typically better than none.

[rfourquet]

That was an improbable example, but which still could happen: let's say I'm writing a wrapper RNG: for floating points, I want it to wrap FloatRNG, and wrap IntRNG for integers, then

rand{T<:Union(Float16,Float32,Float64)}(r::ThisRNG, ::Type{T}) = rand(r.floatrng, T)

is the normal way of writing the wrapping, but the Base generic implementation would get in the way.
But you are totally right that the definition for rand_ui**_raw can be relaxed to AbstractRNG, as there is a unique parameter (and no risk of ambiguity). I will update accordingly.
It will be nice to have more RNG to play with (an RNG package?), which will help to come up eventually with a good solution to provide default implementations for certain types in terms of others.

PS: I enjoy a lot your comments! to the point and lead to better code most of the time. Having my code reviewed is new for me and appreciate how useful/efficient it can be.

[ivarne]

My vision for the fallback system is that we should only need a single method to plug a new RNG into the system (eg rand(r::MyRNG, ::Type{UInt64}), or something similar). As far as I can tell, all other types and ranges, can be derived from this simple definition. Naturally you'd often want to implement special array filling methods for speed, and maybe a 52 bit version if you can earn some efficiency that way, but those should be very specific methods that exploits the internal structure of the RNG to be faster (or conserve random bits from a limited randomness pool).

It's sad that we have to resort to @eval loops to specialize for groups of types in cases like this. It seems to me like a general problem with using the first argument to mimic Object Oriented programming in Julia though. It will be an equally annoying problem for the IO system when someone defines a special method for how to write a FloatingPoint to their MyIO < :IO. It would be nice if one could somehow further configure the method dispatch system, but we should be very careful about exposing more user visible complexity.

[ViralBShah]

It would certainly be nice if the implementation for Float16 and Float32 here could be for AbstractRNG. Is there something we can do here to fix this particular issue, and a more general solution that we adopt in the near future?

We certainly could use an RNG.jl package that adds more RNGs to julia, and helps us refine the interfaces in base.

[rfourquet]

I updated to make rand_ui**_raw applicable to AbstractRNG.

I agree that it would be better to provide rand{T<:Union(Float16, FLoat32)}(::AbstractRNG, ::Type{T}), but I'm not sure yet what is the best solution. Also the problem is not specific to floating types, e.g.: rand(r::AbstractRNG, ::Type{UInt128}) = rand(r, UInt64) % UInt128 << 64 $ rand(r, UInt64) could be provided, but not like this (it would e.g. be ambiguous with RandomDevice methods).

Within a stagedfunction rand{T}(r::AbstractRNG, t::Type{T}), which I believe would not trigger ambiguities, it's probably possible via methods to know what are the "native" types that are implemented for a given RNG and then select the best implementation for type T.

[ViralBShah]

I am inclined to merge this to fix the reported issue, and figure out the APIs in future work. @ivarne would that be ok?

[ivarne]

I don't like postponing API discussions, but as this actually fixes a real problem in a correct way for the existing RNGs, I think we can merge this now.

I see the potential ambiguity problem @rfourquet is concerned about, but I don't think he makes the right tradeoff. I might be wrong, and it would be great if others could comment on the issue.

[ViralBShah]

I too would prefer that if you have a new RNG, all you have to do is provide a method for UInt64 and everything else is generic and automatically works out. Producing Float16 from RandomDevice should ideally just work out. I suspect we are awfully close to it, and some refactoring is in order.

dSFMT is a bit different in that it produces Float64, but that's a small detail that can be accommodated in the refactor.

[rfourquet]

I'm all for having everything work out automatically, but my reserve comes from having had more that once ambiguity warning caused by a method in base which is not enough constrained. As long as there is no way of telling the compiler that a method should take precedence over another one (i.e. explicitly playing with the "hierarchy" of methods), I think a method f(A::AbstractSmthg, B::Concrete) must be defined only if we know that f(A::ConcreteSmthg, Union(Concrete, AnotherConcrete) should not (have to) be defined in user code. If we don't know, then another way of defining f must be found, but I lack experience with Julia to know the good way.
In this particular case, I think we could define something like

stagedfunction rand{T<:Union(Float16,Float32)}(r::AbstractRNG, ::Type{T})
    T===Float16 ? :(rand_Float16(r)) :(rand_Float32(r)))
end

but I wish we did the same for all types.
One example of non-trivial fall back is: rand_ui52 is used in randn, so it is useful to define it for any RNG in terms of rand(r, Close1Open2). But then rand(r, Close1Open2) can be defined either in terms of rand(r, CloseOpen) or rand(r, UInt64). It can't be decided without knowing more about the possibilities, "traits", of the used RNG, and circular definitions must be avoided. We may also want rand(r::AbstractRNG, CloseOpen) = rand(r, Close1Open2)-1 or rand(r, CloseOpen) = rand(r, Float64)...

Asking the user to define traits for her RNG seems workable but not friendly/robust enough given Julia's introspection capabilities. There is probably a solution with stagedfunction and exploring the method tables (the problem with stagedfunction rand{T<:Number}(r::AbstractRNG, ::Type{T}) is that applicable(rand, r, ::Type{T}) will always be true, so another mean (than applicable) of introspection is needed).

I would see a package RNG.jl as an easier playgroud to devise a satisfying automating mecanism. I guess I could start by properly wrapping rdrand, hopefully within few weeks.

[rfourquet]

I didn't understand @ViralBShah 's "Producing Float16 from RandomDevice should ideally just work out" as I thought it worked, and then saw I forgot to add a shift when I changed the implementation to use Float32 internally (before submitting the PR), very sorry for that! Here is a fix, where I also added a test for the uniformity of distribution, does that look good for merging?

[rfourquet]

As the build looks ok, I'll merge the fix.

[ivarne]

Yes, please do