From ccb11668f19b3f7211fc16983d3c4aa124fbd366 Mon Sep 17 00:00:00 2001
From: Rafael Fourquet <fourquet.rafael@gmail.com>
Date: Tue, 25 Nov 2014 14:06:40 +0530
Subject: [PATCH] rename randmtzig_exprnd -> randexp faster randexp (separate
 unlikely branch + @inline) add array versions of randexp and export these
 functions

---
 base/exports.jl     |  2 ++
 base/random.jl      | 48 +++++++++++++++++++++++++++++++--------------
 doc/stdlib/base.rst |  8 ++++++++
 test/random.jl      |  9 +++++++++
 4 files changed, 52 insertions(+), 15 deletions(-)

diff --git a/base/exports.jl b/base/exports.jl
index 7e2f9340909f9..8d614c7e81224 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -905,6 +905,8 @@ export
     randbool,
     randn!,
     randn,
+    randexp!,
+    randexp,
     srand,
 
 # bigfloat & precision
diff --git a/base/random.jl b/base/random.jl
index 2353f2bd50e98..620ef46adef0d 100644
--- a/base/random.jl
+++ b/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
diff --git a/doc/stdlib/base.rst b/doc/stdlib/base.rst
index 0771bb78db31a..3ab5822e46ca1 100644
--- a/doc/stdlib/base.rst
+++ b/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
 ------
 
diff --git a/test/random.jl b/test/random.jl
index 3bf0d87fec087..c42f9860f0b42 100644
--- a/test/random.jl
+++ b/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