RFC: Adds random number generation for BigInts and BigFloats

base/exports.jl

@@ -21,6 +21,7 @@ export
     Bidiagonal,
     BigFloat,
     BigInt,
+    BigRNG,
     BitArray,
     BitMatrix,
     BitVector,

base/gmp.jl

@@ -1,13 +1,13 @@
 module GMP
 
-export BigInt
+export BigInt, BigRNG
 
 import Base: *, +, -, /, <, <<, >>, >>>, <=, ==, >, >=, ^, (~), (&), (|), ($),
              binomial, cmp, convert, div, divrem, factorial, fld, gcd, gcdx, lcm, mod,
              ndigits, promote_rule, rem, show, isqrt, string, isprime, powermod,
              widemul, sum, trailing_zeros, trailing_ones, count_ones, base, parseint,
              serialize, deserialize, bin, oct, dec, hex, isequal, invmod,
-             prevpow2, nextpow2
+             prevpow2, nextpow2, rand, rand!, srand
 
 type BigInt <: Integer
     alloc::Cint
@@ -421,4 +421,87 @@ widemul{T<:Integer}(x::T, y::T) = BigInt(x)*BigInt(y)
 prevpow2(x::BigInt) = x < 0 ? -prevpow2(-x) : (x <= 2 ? x : one(BigInt) << (ndigits(x, 2)-1))
 nextpow2(x::BigInt) = x < 0 ? -nextpow2(-x) : (x <= 2 ? x : one(BigInt) << ndigits(x-1, 2))
 
+# Random number generation
+type BigRNG <: AbstractRNG
+    seed_alloc::Cint
+    seed_size::Cint
+    seed_d::Ptr{Void}
+    alg::Cint
+    alg_data::Ptr{Void}
+
+    function BigRNG()
+        randstate = new(zero(Cint), zero(Cint), C_NULL, zero(Cint), C_NULL)
+        # The default algorithm in GMP is currently MersenneTwister
+        ccall((:__gmp_randinit_default, :libgmp), Void, (Ptr{BigRNG},),
+          &randstate)
+        finalizer(randstate, BigRNG_clear)
+        randstate
+    end
+end
+
+# Initialize with seed
+function BigRNG(seed::Uint64)
+    randstate = BigRNG()
+    ccall((:__gmp_randseed_ui, :libgmp), Void, (Ptr{BigRNG}, Culong),
+      &randstate, seed)
+    randstate
+end
+function BigRNG(seed::BigInt)
+    randstate = BigRNG()
+    ccall((:__gmp_randseed, :libgmp), Void, (Ptr{BigRNG}, Ptr{BigInt}),
+      &randstate, &seed)
+    randstate
+end
+
+function BigRNG_clear(x::BigRNG)
+    ccall((:__gmp_randclear, :libgmp), Void, (Ptr{BigRNG},), &x)
+end
+
+function rand(::Type{BigInt}, randstate::BigRNG, n::Integer)
+    m = convert(BigInt, n)
+    randu(randstate, m)
+end
+
+function rand(r::Range1{BigInt})
+    ulen = convert(BigInt, length(r))
+    convert(BigInt, first(r) + randu(ulen))
+end
+
+function rand!(r::Range1{BigInt}, A::Array{BigInt})
+    for i = 1:length(A)
+        A[i] = rand(r)
+    end
+    A
+end
+
+rand(r::Range1{BigInt}, dims::Dims) = rand!(r, Array(BigInt, dims))
+rand(r::Range1{BigInt}, dims::Int...) = rand!(r, Array(BigInt, dims...))
+
+function randu(randstate::BigRNG, k::BigInt)
+    z = BigInt()
+    ccall((:__gmpz_urandomm, :libgmp), Void,
+      (Ptr{BigInt}, Ptr{BigRNG}, Ptr{BigInt}), &z, &randstate, &k)
+    z
+end
+randu(k::BigInt) = randu(Base.Random.DEFAULT_BIGRNG, k)
+
+srand(r::BigRNG, seed) = srand(r, convert(BigInt, seed))
+function srand(randstate::BigRNG, seed::Vector{Uint32})
+    s = big(0)
+    for i in seed
+        s = s << 32 + i
+    end
+    srand(randstate, s)
+end
+function srand(randstate::BigRNG, seed::Uint64)
+    ccall((:__gmp_randseed_ui, :libgmp), Void, (Ptr{BigRNG}, Culong),
+      &randstate, seed)
+    return
+end
+function srand(randstate::BigRNG, seed::BigInt)
+    ccall((:__gmp_randseed, :libgmp), Void, (Ptr{BigRNG}, Ptr{BigInt}),
+      &randstate, &seed)
+    return
+end
+
 end # module

base/mpfr.jl

@@ -19,7 +19,7 @@ import
         gamma, lgamma, digamma, erf, erfc, zeta, log1p, airyai, iceil, ifloor,
         itrunc, eps, signbit, sin, cos, tan, sec, csc, cot, acos, asin, atan,
         cosh, sinh, tanh, sech, csch, coth, acosh, asinh, atanh, atan2,
-        serialize, deserialize, inf, nan, hash, cbrt
+        serialize, deserialize, inf, nan, hash, cbrt, rand, rand!, randn, randn!
 
 import Base.Math.lgamma_r
 
@@ -739,4 +739,41 @@ function hash(x::BigFloat)
     h
 end
 
+# RNG-related functions
+function rand(::Type{BigFloat}, randstate::BigRNG)
+    z = BigFloat()
+    ccall((:mpfr_urandom,:libmpfr), Int32,
+          (Ptr{BigFloat}, Ptr{BigRNG}, Int32),
+           &z, &randstate, ROUNDING_MODE[end])
+    z
+end
+rand(::Type{BigFloat}) = rand(BigFloat, Base.Random.DEFAULT_BIGRNG)
+
+function rand!(r::BigRNG, A::Array{BigFloat})
+    for i = 1:length(A)
+        A[i] = rand(BigFloat, r)
+    end
+    A
+end
+rand!(A::Array{BigFloat}) = rand!(Base.Random.DEFAULT_BIGRNG, A)
+
+function randn(::Type{BigFloat}, randstate::BigRNG)
+    z = BigFloat()
+    ccall((:mpfr_grandom,:libmpfr), Int32,
+          (Ptr{BigFloat}, Ptr{BigFloat}, Ptr{BigRNG}, Int32),
+           &z, C_NULL, &randstate, ROUNDING_MODE[end])
+    z
+end
+randn(::Type{BigFloat}) = randn(BigFloat, Base.Random.DEFAULT_BIGRNG)
+randn(::Type{BigFloat}, dims::Dims) = randn!(Array(BigFloat, dims))
+randn(::Type{BigFloat}, dims::Int...) = randn!(Array(BigFloat, dims...))
+
+function randn!(r::BigRNG, A::Array{BigFloat})
+    for i = 1:length(A)
+        A[i] = randn(BigFloat, r)
+    end
+    A
+end
+randn!(A::Array{BigFloat}) = randn!(Base.Random.DEFAULT_BIGRNG, A)
+
 end #module

base/random.jl

@@ -54,7 +54,8 @@ function librandom_init()
         seed = reinterpret(Uint64, time())
         seed = bitmix(seed, uint64(getpid()))
         try
-            seed = bitmix(seed, parseint(Uint64, readall(`ifconfig` |> `sha1sum`)[1:40], 16))
+            @linux_only seed = bitmix(seed, parseint(Uint64, readall(`ifconfig` |> `sha1sum`)[1:40], 16))
+            @osx_only seed = bitmix(seed, parseint(Uint64, readall(`ifconfig` |> `shasum`)[1:40], 16))
         catch
             # ignore
         end
@@ -75,6 +76,12 @@ end
 function srand(seed::Vector{Uint32})
     global RANDOM_SEED = seed
     dsfmt_gv_init_by_array(seed)
+    s = big(0)
+    for i in Base.Random.RANDOM_SEED
+        s = s << 32 + i
+    end
+    global DEFAULT_BIGRNG = BigRNG(s)
+    return
 end
 srand(n::Integer) = srand(make_seed(n))
 

doc/helpdb.jl

@@ -3885,7 +3885,8 @@ popdisplay(d::Display)
 
 ("Mathematical Functions","Base","gcd","gcd(x, y)
 
-   Greatest common (positive) divisor (or zero if x and y are both zero).
+   Greatest common (positive) divisor (or zero if x and y are both
+   zero).
 
 "),
 
@@ -3897,8 +3898,8 @@ popdisplay(d::Display)
 
 ("Mathematical Functions","Base","gcdx","gcdx(x, y)
 
-   Greatest common (positive) divisor, also returning integer coefficients \"u\"
-   and \"v\" that solve \"ux+vy == gcd(x,y)\"
+   Greatest common (positive) divisor, also returning integer
+   coefficients \"u\" and \"v\" that solve \"ux+vy == gcd(x,y)\"
 
 "),
 
@@ -4769,7 +4770,8 @@ popdisplay(d::Display)
    Seed the RNG with a \"seed\", which may be an unsigned integer or a
    vector of unsigned integers. \"seed\" can even be a filename, in
    which case the seed is read from a file. If the argument \"rng\" is
-   not provided, the default global RNG is seeded.
+   not provided, the default global RNG and the default BigRNG are
+   seeded.
 
 "),
 
@@ -4781,9 +4783,24 @@ popdisplay(d::Display)
 
 "),
 
+("Random Numbers","Base","BigRNG","BigRNG([seed])
+
+   Create a \"BigRNG\" RNG object, used exclusively to generate random
+   BigInts and BigFloats. Different RNG objects can have their own
+   seeds, which may be useful for generating different streams of
+   random numbers.
+
+"),
+
 ("Random Numbers","Base","rand","rand()
 
-   Generate a \"Float64\" random number uniformly in [0,1)
+   Generate a \"Float64\" random number uniformly in [0,1).
+
+"),
+
+("Random Numbers","Base","rand","rand(BigFloat)
+
+   Generate a \"BigFloat\" random number uniformly in [0,1).
 
 "),
 
@@ -4803,9 +4820,24 @@ popdisplay(d::Display)
 
 "),
 
+("Random Numbers","Base","rand","rand(BigFloat, rng::AbstractRNG[, dims...])
+
+   Generate a random \"BigFloat\" number or array of the size
+   specified by dims, using the specified RNG object. Currently,
+   \"BigRNG\" is the only available Random Number Generator (RNG) for
+   BigFloats, which may be seeded using srand.
+
+"),
+
 ("Random Numbers","Base","rand","rand(dims or [dims...])
 
-   Generate a random \"Float64\" array of the size specified by dims
+   Generate a random \"Float64\" array of the size specified by dims.
+
+"),
+
+("Random Numbers","Base","rand","rand(BigFloat, dims or [dims...])
+
+   Generate a random \"BigFloat\" array of the size specified by dims.
 
 "),
 
@@ -4846,13 +4878,28 @@ popdisplay(d::Display)
 
 "),
 
+("Random Numbers","Base","randn","randn(BigFloat, dims or [dims...])
+
+   Generate a normally-distributed random BigFloat with mean 0 and
+   standard deviation 1. Optionally generate an array of normally-
+   distributed random BigFloats.
+
+"),
+
 ("Random Numbers","Base","randn!","randn!(A::Array{Float64, N})
 
    Fill the array A with normally-distributed (mean 0, standard
    deviation 1) random numbers. Also see the rand function.
 
 "),
 
+("Random Numbers","Base","randn!","randn!(A::Array{BigFloat, N})
+
+   Fill the array A with normally-distributed (mean 0, standard
+   deviation 1) random BigFloats. Also see the rand function.
+
+"),
+
 ("Random Numbers","Base","randsym","randsym(n)
 
    Generate a \"nxn\" symmetric array of normally-distributed random

doc/stdlib/base.rst

@@ -3206,19 +3206,27 @@ The `BigFloat` type implements arbitrary-precision floating-point aritmetic usin
 Random Numbers
 --------------
 
-Random number generateion in Julia uses the `Mersenne Twister library <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/#dSFMT>`_. Julia has a global RNG, which is used by default. Multiple RNGs can be plugged in using the ``AbstractRNG`` object, which can then be used to have multiple streams of random numbers. Currently, only ``MersenneTwister`` is supported.
+Random number generation in Julia uses the `Mersenne Twister library <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/#dSFMT>`_. Julia has a global RNG, which is used by default. Multiple RNGs can be plugged in using the ``AbstractRNG`` object, which can then be used to have multiple streams of random numbers. Currently, only ``MersenneTwister`` is supported. For the generation of BigInts and BigFloats, GMP's own Mersenne Twister implementation is used instead.
 
 .. function:: srand([rng], seed)
 
-   Seed the RNG with a ``seed``, which may be an unsigned integer or a vector of unsigned integers. ``seed`` can even be a filename, in which case the seed is read from a file. If the argument ``rng`` is not provided, the default global RNG is seeded.
+   Seed the RNG with a ``seed``, which may be an unsigned integer or a vector of unsigned integers. ``seed`` can even be a filename, in which case the seed is read from a file. If the argument ``rng`` is not provided, the default global RNG and the default BigRNG are seeded.
 
 .. function:: MersenneTwister([seed])
 
    Create a ``MersenneTwister`` RNG object. Different RNG objects can have their own seeds, which may be useful for generating different streams of random numbers.
 
+.. function:: BigRNG([seed])
+
+   Create a ``BigRNG`` RNG object, used exclusively to generate random BigInts and BigFloats. Different RNG objects can have their own seeds, which may be useful for generating different streams of random numbers.
+
 .. function:: rand()
 
-   Generate a ``Float64`` random number uniformly in [0,1)
+   Generate a ``Float64`` random number uniformly in [0,1).
+
+.. function:: rand(BigFloat)
+
+   Generate a ``BigFloat`` random number uniformly in [0,1).
 
 .. function:: rand!([rng], A)
 
@@ -3228,9 +3236,17 @@ Random number generateion in Julia uses the `Mersenne Twister library <http://ww
 
    Generate a random ``Float64`` number or array of the size specified by dims, using the specified RNG object. Currently, ``MersenneTwister`` is the only available Random Number Generator (RNG), which may be seeded using srand.
 
+.. function:: rand(BigFloat, rng::AbstractRNG, [dims...])
+
+   Generate a random ``BigFloat`` number or array of the size specified by dims, using the specified RNG object. Currently, ``BigRNG`` is the only available Random Number Generator (RNG) for BigFloats, which may be seeded using srand.
+
 .. function:: rand(dims or [dims...])
 
-   Generate a random ``Float64`` array of the size specified by dims
+   Generate a random ``Float64`` array of the size specified by dims.
+
+.. function:: rand(BigFloat, dims or [dims...])
+
+   Generate a random ``BigFloat`` array of the size specified by dims.
 
 .. function:: rand(Int32|Uint32|Int64|Uint64|Int128|Uint128, [dims...])
 
@@ -3252,10 +3268,18 @@ Random number generateion in Julia uses the `Mersenne Twister library <http://ww
 
    Generate a normally-distributed random number with mean 0 and standard deviation 1. Optionally generate an array of normally-distributed random numbers.
 
+.. function:: randn(BigFloat, dims or [dims...])
+
+   Generate a normally-distributed random BigFloat with mean 0 and standard deviation 1. Optionally generate an array of normally-distributed random BigFloats.
+
 .. function:: randn!(A::Array{Float64,N})
 
    Fill the array A with normally-distributed (mean 0, standard deviation 1) random numbers. Also see the rand function.
 
+.. function:: randn!(A::Array{BigFloat,N})
+
+   Fill the array A with normally-distributed (mean 0, standard deviation 1) random BigFloats. Also see the rand function.
+
 .. function:: randsym(n)
 
    Generate a ``nxn`` symmetric array of normally-distributed random numbers with mean 0 and standard deviation 1.

test/random.jl

@@ -17,3 +17,25 @@ if sizeof(Int32) < sizeof(Int)
     @test typeof(r) == Int32
     @test -1 <= r <= typemax(Int32)
 end
+
+range = big(typemax(Int128)):(big(typemax(Int128)) + 10000)
+@test rand(range) != rand(range)
+@test big(typemax(Int128)) <= rand(range) <= big(typemax(Int128)) + 10000
+r = rand(range, 1, 2)
+@test size(r) == (1, 2)
+@test typeof(r) == Matrix{BigInt}
+
+srand(0)
+r = rand(range)
+f = rand(BigFloat)
+srand(0)
+@test rand(range) == r
+@test rand(BigFloat) == f
+
+@test rand(BigFloat) != rand(BigFloat)
+@test 0.0 <= rand(BigFloat) < 1.0
+@test typeof(rand(BigFloat)) == BigFloat
+
+r = rand(BigFloat, 1, 3)
+@test size(r) == (1, 3)
+@test typeof(r) == Matrix{BigFloat}