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.
JuliaLang · Jul 8, 2011 · 411bc0b · 411bc0b
1 parent 5dd15dd
commit 411bc0b
Show file tree

Hide file tree

Showing 4 changed files with 60 additions and 32 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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,21 +59,22 @@ 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
 # 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