Merge pull request #9144 from JuliaLang/rf/randexp

faster randmtzig_exprnd (separate unlikely branch + inline) and rename it "randexp"

base/exports.jl

@@ -905,6 +905,8 @@ export
     randbool,
     randn!,
     randn,
+    randexp!,
+    randexp,
     srand,
 
 # bigfloat & precision

base/random.jl

@@ -5,9 +5,10 @@ using Base.dSFMT
 export srand,
        rand, rand!,
        randn, randn!,
+       randexp, randexp!,
        randbool,
-       AbstractRNG, RNG, MersenneTwister,
-       randmtzig_exprnd
+       AbstractRNG, RNG, MersenneTwister
+
 
 abstract AbstractRNG
 
@@ -984,23 +985,40 @@ randn(dims::Int...) = randn!(Array(Float64, dims...))
 randn(rng::MersenneTwister, dims::Dims) = randn!(rng, Array(Float64, dims))
 randn(rng::MersenneTwister, dims::Int...) = randn!(rng, Array(Float64, dims...))
 
-function randmtzig_exprnd(rng::MersenneTwister=GLOBAL_RNG)
+@inline function randexp(rng::MersenneTwister=GLOBAL_RNG)
     @inbounds begin
-        while true
-            ri = rand_ui52(rng)
-            idx = ri & 0xFF
-            x = ri*we[idx+1]
-            if ri < ke[idx+1]
-                return x # 98.9% of the time we return here 1st try
-            elseif idx == 0
-                return ziggurat_exp_r - log(rand(rng))
-            elseif (fe[idx] - fe[idx+1])*rand(rng) + fe[idx+1] < exp(-x)
-                return x # return from the triangular area
-            end
-        end
+        ri = rand_ui52(rng)
+        idx = ri & 0xFF
+        x = ri*we[idx+1]
+        ri < ke[idx+1] && return x # 98.9% of the time we return here 1st try
+        return randexp_unlikely(rng, idx, x)
+    end
+end
+
+function randexp_unlikely(rng, idx, x)
+    @inbounds if idx == 0
+        return ziggurat_exp_r - log(rand(rng))
+    elseif (fe[idx] - fe[idx+1])*rand(rng) + fe[idx+1] < exp(-x)
+        return x # return from the triangular area
+    else
+        return randexp(rng)
+    end
+end
+
+function randexp!(rng::MersenneTwister, A::Array{Float64})
+    for i = 1:length(A)
+        @inbounds A[i] = randexp(rng)
     end
+    A
 end
 
+randexp!(A::Array{Float64}) = randexp!(GLOBAL_RNG, A)
+randexp(dims::Dims) = randexp!(Array(Float64, dims))
+randexp(dims::Int...) = randexp!(Array(Float64, dims))
+randexp(rng::MersenneTwister, dims::Dims) = randexp!(rng, Array(Float64, dims))
+randexp(rng::MersenneTwister, dims::Int...) = randexp!(rng, Array(Float64, dims))
+
+
 ## random UUID generation
 
 immutable UUID

doc/stdlib/base.rst

@@ -4119,6 +4119,14 @@ A ``MersenneTwister`` RNG can generate random numbers of the following types: ``
 
    Fill the array A with normally-distributed (mean 0, standard deviation 1) random numbers. Also see the rand function.
 
+.. function:: randexp([rng], [dims...])
+
+   Generate a random number according to the exponential distribution with scale 1. Optionally generate an array of such random numbers.
+
+.. function:: randexp!([rng], A::Array{Float64,N})
+
+   Fill the array A with random numbers following the exponential distribution (with scale 1).
+
 Arrays
 ------
 

test/random.jl

@@ -80,11 +80,20 @@ if sizeof(Int32) < sizeof(Int)
 
 end
 
+randn()
 randn(100000)
 randn!(Array(Float64, 100000))
+randn(MersenneTwister(10))
 randn(MersenneTwister(10), 100000)
 randn!(MersenneTwister(10), Array(Float64, 100000))
 
+randexp()
+randexp(100000)
+randexp!(Array(Float64, 100000))
+randexp(MersenneTwister(10))
+randexp(MersenneTwister(10), 100000)
+randexp!(MersenneTwister(10), Array(Float64, 100000))
+
 # Test ziggurat tables
 ziggurat_table_size = 256
 nmantissa           = int64(2)^51 # one bit for the sign