randn() and BigFloat

[AshtonSBradley]

In the docs:

https://docs.julialang.org/en/latest/stdlib/Random/

A MersenneTwister or RandomDevice RNG can generate random numbers of the following types: Float16, Float32, Float64, BigFloat, Bool, Int8, UInt8, Int16, UInt16, Int32, UInt32, Int64, UInt64, Int128, UInt128, BigInt (or complex numbers of those types). Random floating point numbers are generated uniformly in [0,1). As BigInt represents unbounded integers, the interval must be specified (e.g. rand(big.(1:6))).

However,

julia> randn(BigFloat)
ERROR: MethodError: no method matching randn(::MersenneTwister, ::Type{BigFloat})
Closest candidates are:
  randn(::AbstractRNG, ::Type{T}, ::Tuple{Vararg{Int64,N}} where N) where T at random.jl:1387
  randn(::AbstractRNG, ::Type{T}, ::Integer, ::Integer...) where T at random.jl:1389
  randn(::AbstractRNG) at random.jl:1257
  ...
Stacktrace:
 [1] randn(::Type{BigFloat}) at ./random.jl:1374

Is there some special syntax to call randn() for BigFloat? Is it supported?

[ChrisRackauckas]

@AshtonSBradley thanks for opening this. One thing to note though is that we use a lot of RandomNumbers.jl for the application here, so if this gets into v0.7 we can modify the default for the BigFloat case but I hope that whatever solution is found can be downstreamed to some of RandomNumbers.jl as well.

[stevengj]

I don't think we have ever supported randn for BigFloat? See #17629. So this seems like a documentation bug.

[rfourquet]

Yes, the introduction of the Random documentation is quite old and should be updated, the list of supported types implicitly assumes uniform distribution. On the other hand, the doc of randn is quite clear on what is available.

I'm in the process of re-rewriting this doc, so will clarify this point while I'm at it, thanks.

Chris, I don't really understand what you said, but if this ever gets implemented in Random, it should support out of the box other RNGs, e.g. from RandomNumbers.jl.

[ChrisRackauckas]

This isn't just a documentation bug. This was noticed by Ashton when trying to solve stochastic differential equations with BigFloats, and instead we're trying to find a way to get random normal BigFloats. While there does seem to be a documentation bug here (it saying it exists when it doesn't), the real underlying thing is a feature request for random BigFloats.

Chris, I don't really understand what you said, but if this ever gets implemented in Random, it should support out of the box other RNGs, e.g. from RandomNumbers.jl.

Part of that response was continuing the conversation which led to this issue, sorry.

I'm not sure you can map a single BigFloat handling to all PRNG algorithms? If you take the current RNGs can you just upscale them to larger digit sizes? For example, for a 256 BigFloat, just take 8 64-bit Floats from the same RNG to build a uniform distribution RNG and use that in the normal distribution? I am not very familiar with all of the RNG stuff, so maybe it is that simple and I just don't know.

[stevengj]

the real underlying thing is a feature request for random BigFloats.

That's a dup of #17629.

[rfourquet]

I'm not sure you can map a single BigFloat handling to all PRNG algorithms?

Maybe not so easily, but this is on the radar to have an API allowing RNGs implementors to plug into already implemented random generation algorithms. For example it's already quite easy to define an RNG and have it generate BigFloat uniformly. Implementing normal generation should have to be done only once.

[ViralBShah]

I thought that the details of the stuff inside the ziggurat are hardcoded, and that it would be different for different precisions.

[stevengj]

There are various non-Ziggurat algorithms, some of which are cited in the other issue.

[andreasnoack]

The documentation is pretty clear about which number types are supported for randn so this is a dup of #17629

[MilesCranmer]

For anybody like me who landed here while Googling for a workaround: the following hack works if you just want your codebase to be generic to BigFloat types without needing a high-precision RNG.

using Random: AbstractRNG

try
    # In case randn(BigFloat) gets defined in the future
    randn(BigFloat)
catch e
    if isa(e, MethodError)
        Base.randn(::Type{BigFloat}, args::Integer...) = BigFloat.(randn(Float64, args...))
        Base.randn(rng::AbstractRNG, ::Type{BigFloat}, args::Integer...) = BigFloat.(randn(rng, Float64, args...))
    else
        throw(e)
    end
end

If you do want a high-precision RNG and don't need it to be fast you could use erfinv and the existing rand(BigFloat):

using Random: AbstractRNG
using SpecialFunctions: erfinv

try
    # In case randn(BigFloat) gets defined in the future
    randn(BigFloat)
catch e
    if isa(e, MethodError)
        gauss_q(p::T) where {T} = sqrt(T(2))*erfinv(2*p-1)
        Base.randn(::Type{BigFloat}, args::Integer...) = gauss_q.(rand(BigFloat, args...))
        Base.randn(rng::AbstractRNG, ::Type{BigFloat}, args::Integer...) = gauss_q.(rand(rng, BigFloat, args...))
    else
        throw(e)
    end
end

[oscardssmith]

Note that in 1.9, randn(BigFloat) already works. (added in #44714)

[MilesCranmer]

Nice! Thanks 😅