Ziggurat stuff compiled and callable now.

external/Makefile

@@ -132,8 +132,8 @@ distclean-fdlibm: clean-fdlibm
 
 ## DSFMT ##
 
-DSFMT_OBJ_TARGET = $(EXTROOTLIB)/libMT.$(SHLIB_EXT)
-DSFMT_OBJ_SOURCE = random/libMT.$(SHLIB_EXT)
+DSFMT_OBJ_TARGET = $(EXTROOTLIB)/librandom.$(SHLIB_EXT)
+DSFMT_OBJ_SOURCE = random/librandom.$(SHLIB_EXT)
 
 compile-dsfmt: $(DSFMT_OBJ_SOURCE)
 install-dsfmt: $(DSFMT_OBJ_TARGET)
@@ -149,14 +149,14 @@ random/jl_random.c: random/dsfmt-$(DSFMT_VER).tar.gz
 	touch $@
 $(DSFMT_OBJ_SOURCE): random/jl_random.c
 	cd random && \
-	$(CC) -O3 -finline-functions -fomit-frame-pointer -DNDEBUG -fno-strict-aliasing --param max-inline-insns-single=1800  -Wmissing-prototypes -Wall  -std=c99 -DDSFMT_MEXP=19937 -fPIC -shared -DDSFMT_DO_NOT_USE_OLD_NAMES jl_random.c -o libMT.$(SHLIB_EXT)
+	$(CC) -O3 -finline-functions -fomit-frame-pointer -DNDEBUG -fno-strict-aliasing --param max-inline-insns-single=1800  -Wmissing-prototypes -Wall  -std=c99 -DDSFMT_MEXP=19937 -fPIC -shared -DDSFMT_DO_NOT_USE_OLD_NAMES -DUSE_X86_32=0 jl_random.c randmtzig.c -o librandom.$(SHLIB_EXT)
 
 $(DSFMT_OBJ_TARGET): $(DSFMT_OBJ_SOURCE)
 	mkdir -p $(EXTROOTLIB)
 	cp $< $@
 
 clean-dsfmt:
-	rm -f random/libMT.$(SHLIB_EXT)
+	rm -f random/librandom.$(SHLIB_EXT)
 distclean-dsfmt: clean-dsfmt
 	cd random && rm -rf *.tar.gz dsfmt-$(DSFMT_VER)
 

external/random/jl_random.h

@@ -15,5 +15,6 @@ uint32_t genrand_int32();
 void randomseed32(uint32_t s);
 void randomseed64(uint64_t s);
 double dsfmt_randn();
+uint32_t dsfmt_gv_genrand_uint32(void);
 
 #endif

external/random/randmtzig.c

@@ -136,6 +136,8 @@
 #include <time.h>
 #include <sys/time.h>
 
+#include "jl_random.h"
+
 typedef int randmtzig_idx_type;
 typedef signed char randmtzig_int8_t;
 typedef unsigned char randmtzig_uint8_t;
@@ -332,12 +334,12 @@ static randmtzig_uint32_t randmt (randmtzig_uint32_t *state)
 /* ===== Uniform generators ===== */
 
 /* Select which 32 bit generator to use */
-#define randi32 randmt
+//#define randi32 randmt
 
 static randmtzig_uint64_t randi53 (randmtzig_uint32_t *state)
 {
-  const randmtzig_uint32_t lo = randi32(state);
-  const randmtzig_uint32_t hi = randi32(state)&0x1FFFFF;
+  const randmtzig_uint32_t lo = dsfmt_gv_genrand_uint32();
+  const randmtzig_uint32_t hi = dsfmt_gv_genrand_uint32()&0x1FFFFF;
 #if HAVE_X86_32
   randmtzig_uint64_t u;
   randmtzig_uint32_t *p = (randmtzig_uint32_t *)&u;
@@ -351,8 +353,8 @@ static randmtzig_uint64_t randi53 (randmtzig_uint32_t *state)
 
 static randmtzig_uint64_t randi54 (randmtzig_uint32_t *state)
 {
-  const randmtzig_uint32_t lo = randi32(state);
-  const randmtzig_uint32_t hi = randi32(state)&0x3FFFFF;
+  const randmtzig_uint32_t lo = dsfmt_gv_genrand_uint32();
+  const randmtzig_uint32_t hi = dsfmt_gv_genrand_uint32()&0x3FFFFF;
 #if HAVE_X86_32
   randmtzig_uint64_t u;
   randmtzig_uint32_t *p = (randmtzig_uint32_t *)&u;
@@ -366,8 +368,8 @@ static randmtzig_uint64_t randi54 (randmtzig_uint32_t *state)
 
 static randmtzig_uint64_t randi64 (randmtzig_uint32_t *state)
 {
-  const randmtzig_uint32_t lo = randi32(state);
-  const randmtzig_uint32_t hi = randi32(state);
+  const randmtzig_uint32_t lo = dsfmt_gv_genrand_uint32();
+  const randmtzig_uint32_t hi = dsfmt_gv_genrand_uint32();
 #if HAVE_X86_32
   randmtzig_uint64_t u;
   randmtzig_uint32_t *p = (randmtzig_uint32_t *)&u;
@@ -382,15 +384,15 @@ static randmtzig_uint64_t randi64 (randmtzig_uint32_t *state)
 /* generates a random number on (0,1)-real-interval */
 static double randu32 (randmtzig_uint32_t *state)
 {
-  return ((double)randi32(state) + 0.5) * (1.0/4294967296.0); 
+  return ((double)dsfmt_gv_genrand_uint32() + 0.5) * (1.0/4294967296.0); 
   /* divided by 2^32 */
 }
 
 /* generates a random number on (0,1) with 53-bit resolution */
 static double randu53 (randmtzig_uint32_t *state) 
 { 
-  const randmtzig_uint32_t a=randi32(state)>>5;
-  const randmtzig_uint32_t b=randi32(state)>>6; 
+  const randmtzig_uint32_t a=dsfmt_gv_genrand_uint32()>>5;
+  const randmtzig_uint32_t b=dsfmt_gv_genrand_uint32()>>6; 
   return(a*67108864.0+b+0.4) * (1.0/9007199254740992.0);
 } 
 
@@ -402,7 +404,7 @@ double randmtzig_randu (randmtzig_uint32_t *state)
 
 randmtzig_uint32_t randmtzig_randi32 (randmtzig_uint32_t *state)
 {
-     return randi32(state);
+     return dsfmt_gv_genrand_uint32();
 }
 
 /* ===== Ziggurat normal and exponential generators ===== */
@@ -574,9 +576,9 @@ randmtzig_randn (randmtzig_uint32_t *state, ZIGINT *ki, ZIGINT *ke, double *wi,
       register randmtzig_uint32_t lo, hi;
       randmtzig_int64_t rabs;
       randmtzig_uint32_t *p = (randmtzig_uint32_t *)&rabs;
-      lo = randi32(state);
+      lo = dsfmt_gv_genrand_uint32();
       idx = lo&0xFF;
-      hi = randi32(state);
+      hi = dsfmt_gv_genrand_uint32();
       si = hi&UMASK;
       p[0] = lo;
       p[1] = hi&0x1FFFFF;
@@ -591,7 +593,7 @@ randmtzig_randn (randmtzig_uint32_t *state, ZIGINT *ki, ZIGINT *ke, double *wi,
       if (rabs < (randmtzig_int64_t)ki[idx])
 #else /* ALLBITS */
       /* 32-bit mantissa */
-      const randmtzig_uint32_t r = randi32(state);
+      const randmtzig_uint32_t r = dsfmt_gv_genrand_uint32();
       const randmtzig_uint32_t rabs = r&LMASK;
       const int idx = (int)(r&0xFF);
       const double x = ((randmtzig_int32_t)r) * wi[idx];

j/random.j

@@ -1,36 +1,38 @@
-libmt = dlopen("libMT")
+librandom = dlopen("librandom")
 
-function mt_init()
+function librandom_init()
     try
         srand("/dev/urandom", 4)
     catch
         srand(uint64(clock()*2.0^32))
     end
+
+    zig_randn_init_by_int(uint32(clock()))
 end
 
-dsfmt_get_min_array_size() = ccall(dlsym(libmt, :dsfmt_get_min_array_size), Int32, ())
+dsfmt_get_min_array_size() = ccall(dlsym(librandom, :dsfmt_get_min_array_size), Int32, ())
 
-dsfmt_randn_reset() = ccall(dlsym(libmt, :dsfmt_randn_reset), Void, ())
+dsfmt_randn_reset() = ccall(dlsym(librandom, :dsfmt_randn_reset), Void, ())
 
-srand(seed::Uint32) = (ccall(dlsym(libmt, :dsfmt_gv_init_gen_rand), Void, (Uint32, ), seed);
+srand(seed::Uint32) = (ccall(dlsym(librandom, :dsfmt_gv_init_gen_rand), Void, (Uint32, ), seed);
                        dsfmt_randn_reset())
 
 srand(seed::Uint64) = srand([uint32(seed),uint32(seed>>32)])
 
 function srand(seed::Vector{Uint32})
-    ccall(dlsym(libmt, :dsfmt_gv_init_by_array),
+    ccall(dlsym(librandom, :dsfmt_gv_init_by_array),
           Void, (Ptr{Uint32}, Int32),
           seed, int32(length(seed)))
     dsfmt_randn_reset()
 end
 
 randf() = float32(rand())
 
-rand() = ccall(dlsym(libmt, :dsfmt_gv_genrand_open_open), Float64, ())
+rand() = ccall(dlsym(librandom, :dsfmt_gv_genrand_open_open), Float64, ())
 
-randui32() = ccall(dlsym(libmt, :dsfmt_gv_genrand_uint32), Uint32, ())
+randui32() = ccall(dlsym(librandom, :dsfmt_gv_genrand_uint32), Uint32, ())
 
-randn() = ccall(dlsym(libmt, :dsfmt_randn), Float64, ())
+randn() = ccall(dlsym(librandom, :dsfmt_randn), Float64, ())
 
 randbit() = randui32()&1
 
@@ -44,10 +46,10 @@ function dsfmt_fill_array_open_open(A::Array{Float64})
         end
     else
         if isodd(n)
-            ccall(dlsym(libmt, :dsfmt_gv_fill_array_open_open), Void, (Ptr{Void}, Int32), A, int32(n-1))
+            ccall(dlsym(librandom, :dsfmt_gv_fill_array_open_open), Void, (Ptr{Void}, Int32), A, int32(n-1))
             A[n] = rand()
         else
-            ccall(dlsym(libmt, :dsfmt_gv_fill_array_open_open), Void, (Ptr{Void}, Int32), A, int32(n))
+            ccall(dlsym(librandom, :dsfmt_gv_fill_array_open_open), Void, (Ptr{Void}, Int32), A, int32(n))
         end
     end
     return A
@@ -97,6 +99,35 @@ end
 # random integer from 1 to n
 randint(n::Int) = randint(one(n), n)
 
+## Normally distributed random numbers using Ziggurat algorithm
+
+ZT_SIZE = 256
+ZT_STATE = Array(Uint32, 628)
+ki = Array(Uint64,  ZT_SIZE)
+ke = Array(Uint64,  ZT_SIZE)
+wi = Array(Float64, ZT_SIZE)
+fi = Array(Float64, ZT_SIZE)
+we = Array(Float64, ZT_SIZE)
+fe = Array(Float64, ZT_SIZE)
+
+zig_randn_init_by_int(x::Uint32) = ccall(dlsym(librandom, :randmtzig_init_by_int), Void, (Uint32, Ptr{Uint32},), x, ZT_STATE)
+
+#zig_randn_init_by_entropy() = ccall(dlsym(librandom, :randmtzig_init_by_entropy), Void, (Ptr{Uint32},), ZT_STATE)
+
+zig_randn() = ccall(dlsym(librandom, :randmtzig_randn), Float64, 
+                    (Ptr{Uint32}, Ptr{Uint64}, Ptr{Uint64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
+                    ZT_STATE, ki, ke, wi, fi, we, fe)
+
+function zig_randn(dims::Dims)
+    A = Array(Float64, dims)
+    ccall(dlsym(librandom, :randmtzig_fill_drandn), Void,
+          (Int32, Ptr{Float64}, Ptr{Uint32}, Ptr{Uint64}, Ptr{Uint64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
+          int32(numel(A)), A, ZT_STATE, ki, ke, wi, fi, we, fe)
+    return A
+end
+
+zig_randn(dims::Size...) = zig_randn(dims)
+
 ## Arrays of random numbers
 
 function rand(dims::Dims)

j/start_image.j

@@ -11,7 +11,7 @@ stderr_stream = fdio(ccall(:jl_stderr, Int32, ()))
 libc = ccall(:jl_load_dynamic_library, Ptr{Void}, (Ptr{Uint8},), C_NULL)
 libm = dlopen("libm")
 libfdm = dlopen("libfdm")
-libmt = dlopen("libMT"); mt_init();
+librandom = dlopen("librandom"); librandom_init();
 libpcre = dlopen("libpcre")
 libBLAS = dlopen("libLAPACK")
 libLAPACK = libBLAS

j/sysimg.j

@@ -161,7 +161,7 @@ compile_hint(cmp, (Int32, Int32))
 compile_hint(min, (Int32, Int32))
 compile_hint(==, (ASCIIString, ASCIIString))
 compile_hint(arg_gen, (ASCIIString,))
-compile_hint(mt_init, ())
+compile_hint(librandom_init, ())
 compile_hint(srand, (ASCIIString, Long))
 compile_hint(open, (ASCIIString, Bool, Bool, Bool, Bool))
 compile_hint(srand, (Uint64,))