Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 1 commit into from
Nov 30, 2014
Merged

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

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions base/exports.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -905,6 +905,8 @@ export
randbool,
randn!,
randn,
randexp!,
randexp,
srand,

# bigfloat & precision
Expand Down
48 changes: 33 additions & 15 deletions base/random.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
8 changes: 8 additions & 0 deletions doc/stdlib/base.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
------

Expand Down
9 changes: 9 additions & 0 deletions test/random.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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))

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You should probably @test something too. This code will not verify anything other than that it doesn't throw exceptions. I know the results are random, but the size, type and range is not.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here I only duplicated what already existed for randn, but I agree that testing values should be added.
I personally find that it is good to test somewhere only the API (smaller sizes than 100000 would do it though). What we could do here is test that the global version and the version with a passed rng give the same results. I am planning to rewrite #8339 differently, with all API checks grouped together. Then it could be enough to do all the other tests on the global RNG.

# Test ziggurat tables
ziggurat_table_size = 256
nmantissa = int64(2)^51 # one bit for the sign
Expand Down