From 411bc0b59f46a46c779d6ee6705710c664ca3109 Mon Sep 17 00:00:00 2001 From: Viral Shah Date: Fri, 8 Jul 2011 19:12:33 +0530 Subject: [PATCH] Use DSFMT to implement randn. This is done in jl_random.c. About 50% faster now. Use libm sqrt() for now, since it seems well-behaved, until we move to libm64. Switch random.j to use DSFMT for all random number generation. --- external/random/jl_random.c | 23 +++++++++++++ external/random/jl_random.h | 1 + j/math_libm.j | 2 +- j/random.j | 66 ++++++++++++++++++++----------------- 4 files changed, 60 insertions(+), 32 deletions(-) diff --git a/external/random/jl_random.c b/external/random/jl_random.c index 3efb72d9f1a00..4c5f8d32922df 100644 --- a/external/random/jl_random.c +++ b/external/random/jl_random.c @@ -1,9 +1,11 @@ /* random numbers */ + #include "jl_random.h" #include "mt19937ar.c" +#include "dsfmt-2.1/dSFMT.c" double rand_double() { @@ -70,3 +72,24 @@ void randomize() randomseed64(a); } + +double dsfmt_randn() +{ + double s, vre, vim, ure, uim; + + if (randn_next != -42) { + s = randn_next; + randn_next = -42; + return s; + } + do { + ure = dsfmt_gv_genrand_close1_open2(); + uim = dsfmt_gv_genrand_close1_open2(); + vre = 2*ure - 3; + vim = 2*uim - 3; + s = vre*vre + vim*vim; + } while (s >= 1); + s = sqrt(-2*log(s)/s); + randn_next = s * vre; + return s * vim; +} diff --git a/external/random/jl_random.h b/external/random/jl_random.h index 01aa4813f14e3..f40aaa26d01ea 100644 --- a/external/random/jl_random.h +++ b/external/random/jl_random.h @@ -15,5 +15,6 @@ void randomize(); uint32_t genrand_int32(); void randomseed32(uint32_t s); void randomseed64(uint64_t s); +double dsfmt_randn(); #endif diff --git a/j/math_libm.j b/j/math_libm.j index 163f1342e87b5..cea183bfcdc1b 100644 --- a/j/math_libm.j +++ b/j/math_libm.j @@ -50,7 +50,6 @@ macro libfdmfunc_2arg(f) end end -@libfdmfunc_1arg_float sqrt @libfdmfunc_1arg_float cbrt @libfdmfunc_1arg_float sin @libfdmfunc_1arg_float cos @@ -75,6 +74,7 @@ end @libfdmfunc_1arg_float rint @libfdmfunc_1arg_float lgamma +@libmfunc_1arg_float sqrt @libmfunc_1arg_float exp2 @libmfunc_1arg_float nearbyint @libmfunc_1arg_float trunc diff --git a/j/random.j b/j/random.j index 084f89a06545d..cf619f07d94db 100644 --- a/j/random.j +++ b/j/random.j @@ -3,48 +3,52 @@ libmt = dlopen("libMT") randomize() = ccall(dlsym(libmt, :randomize), Void, ()) function mt_init() - randomize() - dsfmt_init_gen_rand(0) + #randomize() + srand(0) end ### DSFMT ### -dsfmt_init_gen_rand(seed::Union(Int32,Uint32)) = ccall(dlsym(libmt, :dsfmt_gv_init_gen_rand), - Void, (Uint32, ), uint32(seed)) - dsfmt_get_min_array_size() = ccall(dlsym(libmt, :dsfmt_get_min_array_size), Int32, ()) -dsfmt_genrand_open_open() = ccall(dlsym(libmt, :dsfmt_gv_genrand_open_open), Float64, ()) +srand(seed::Union(Int32,Uint32)) = ccall(dlsym(libmt, :dsfmt_gv_init_gen_rand), + Void, (Uint32, ), uint32(seed)) + +randf() = float32(rand()) + +rand() = ccall(dlsym(libmt, :dsfmt_gv_genrand_open_open), Float64, ()) + +randui32() = ccall(dlsym(libmt, :dsfmt_gv_genrand_uint32), Uint32, ()) -dsfmt_genrand_uint32() = ccall(dlsym(libmt, :dsfmt_gv_genrand_uint32), Uint32, ()) +randn() = ccall(dlsym(libmt, :dsfmt_randn), Float64, ()) function dsfmt_fill_array_open_open(A::Array{Float64}) n = numel(A) if (n <= dsfmt_get_min_array_size()) for i=1:numel(A) - A[i] = dsfmt_genrand_open_open() + A[i] = rand() end else - ccall(dlsym(libmt, :dsfmt_gv_fill_array_open_open), Void, (Ptr{Void}, Int32), A, n) + if isodd(n) + ccall(dlsym(libmt, :dsfmt_gv_fill_array_open_open), Void, (Ptr{Void}, Int32), A, n-1) + A[n] = rand() + else + ccall(dlsym(libmt, :dsfmt_gv_fill_array_open_open), Void, (Ptr{Void}, Int32), A, n) + end end return A end -# rand() = dsfmt_genrand_open_open() -# randf() = float32(rand()) -# randui32() = dsfmt_genrand_uint32() -# srand(s) = dsfmt_init_gen_rand(s) - # jl_randn_next = -42.0 # function randn() # global jl_randn_next - +# # if (jl_randn_next != -42.0) # s = jl_randn_next # jl_randn_next = -42.0 # return s # end - +# # s = 1.0 # vre = 0.0 # vim = 0.0 @@ -55,7 +59,7 @@ end # vim = 2.0*uim - 1.0 # s = vre*vre + vim*vim # end - +# # s = sqrt(-2.0*log(s)/s) # jl_randn_next = s * vre # return s * vim @@ -63,13 +67,14 @@ end ### MT ### +# This is the old code based on the original Mersenne Twister -rand() = ccall(dlsym(libmt, :rand_double), Float64, ()) -randf() = ccall(dlsym(libmt, :rand_float), Float32, ()) -randui32() = ccall(dlsym(libmt, :genrand_int32), Uint32, ()) -randn() = ccall(dlsym(libmt, :randn), Float64, ()) -srand(s::Union(Int32,Uint32)) = ccall(dlsym(libmt, :randomseed32), Void, (Uint32,), uint32(s)) -srand(s::Union(Int64,Uint64)) = ccall(dlsym(libmt, :randomseed64), Void, (Uint64,), uint64(s)) +#rand() = ccall(dlsym(libmt, :rand_double), Float64, ()) +#randf() = ccall(dlsym(libmt, :rand_float), Float32, ()) +#randui32() = ccall(dlsym(libmt, :genrand_int32), Uint32, ()) +#randn() = ccall(dlsym(libmt, :randn), Float64, ()) +#srand(s::Union(Int32,Uint32)) = ccall(dlsym(libmt, :randomseed32), Void, (Uint32,), uint32(s)) +#srand(s::Union(Int64,Uint64)) = ccall(dlsym(libmt, :randomseed64), Void, (Uint64,), uint64(s)) ## Random integers @@ -107,13 +112,13 @@ randint(n::Int) = randint(one(n), n) ## Arrays of random numbers -# function rand(dims::Dims) -# A = Array(Float64, dims) -# dsfmt_fill_array_open_open(A) -# return A -# end +function rand(dims::Dims) + A = Array(Float64, dims) + dsfmt_fill_array_open_open(A) + return A +end -# rand(dims::Size...) = rand(dims) +rand(dims::Size...) = rand(dims) macro rand_matrix_builder(t, f) quote @@ -131,7 +136,6 @@ macro rand_matrix_builder(t, f) end # quote end # macro -@rand_matrix_builder Float64 rand @rand_matrix_builder Float64 randn -@rand_matrix_builder Uint32 randui32 @rand_matrix_builder Float32 randf +@rand_matrix_builder Uint32 randui32