Translate ziggurat algorithms for normal and exponential variates to Julia.

[andreasnoack]

This pull request removes the C file that computes normal variates in favor of a translation to Julia. On julia.mit.edu the implementation in Julia was 25 pct. faster which I can't really explain why. On my MacBook, they are approximately equally fast.

The tables necessary for calculating the normal variates are made in BigFloat and finally translated to Uint64 and Float64 thereby making them identical on 32 and 64 bit systems without storing the hardcoded tables as we done for about a week.

The server is busy with other computations, and therefore I haven't run BigCrush on the normal and exponential variates yet, but I have tried the test that the rand fails, see #6464. Surprisingly, the normal variates don't fail the test.

[JeffBezanson]

I still think I'd rather hard-code the tables; if nothing else this will make it easier to build more of julia without bigfloat support.

[andreasnoack]

Okay. Is there a point in letting the generating function stay such that we can see how they are created?

I can see that Travis doesn't like this. I have removed the random directory because it only included dSMFT. Can someone with Travis knowledge help me here?

[StefanKarpinski]

How about checking the tables against the generating code as a BigFloat test? That way we keep the code around, make sure it keeps working, and don't cause the build process to depend in BigFloats.

On Apr 11, 2014, at 1:04 PM, Andreas Noack Jensen [email protected] wrote:

Okay. Is there a point in letting the generating function stay such that we can see how they are created?

I can see that Travis don't like this. I have removed the random directory because it only included dSMFT. Can someone with Travis knowledge help me here?

—
Reply to this email directly or view it on GitHub.

[staticfloat]

Travis is sad because make is trying to resolve the dsfmt-$(DSFMT_VER)/dSFMT.c dependency before it resolves the dsfmt-$(DSFMT_VER)/config.status dependency. This diff should fix your problems:

diff --git a/deps/Makefile b/deps/Makefile
index f9b5278..6f17fd5 100644
--- a/deps/Makefile
+++ b/deps/Makefile
@@ -645,7 +645,7 @@ dsfmt-$(DSFMT_VER)/config.status: dsfmt-$(DSFMT_VER).tar.gz
        $(TAR) -C dsfmt-$(DSFMT_VER) --strip-components 1 -xf dsfmt-$(DSFMT_VER).tar.gz && \
        cd dsfmt-$(DSFMT_VER) && patch < ../dSFMT.h.patch && patch < ../dSFMT.c.patch
        echo 1 > $@
-$(LIBRANDOM_OBJ_SOURCE): dsfmt-$(DSFMT_VER)/dSFMT.c dsfmt-$(DSFMT_VER)/config.status
+$(LIBRANDOM_OBJ_SOURCE): dsfmt-$(DSFMT_VER)/config.status dsfmt-$(DSFMT_VER)/dSFMT.c
        $(CC) $(CPPFLAGS) $(LIBRANDOM_CFLAGS) $(LDFLAGS) dsfmt-$(DSFMT_VER)/dSFMT.c -o librandom.$(SHLIB_EXT) && \
        $(INSTALL_NAME_CMD)librandom.$(SHLIB_EXT) librandom.$(SHLIB_EXT)
 $(LIBRANDOM_OBJ_TARGET): $(LIBRANDOM_OBJ_SOURCE)

[ViralBShah]

Certainly good to have the generation code in there as well.

[andreasnoack]

I have a problem with defining pointers to the arrays storing the ziggurat tables. First of all, I am in doubt if this ought to be faster that indexing the array with @inbounds, but in practice it made the code exponential code significantly faster and the normal code slightly/insignificantly faster. The problem is that after moving the code from an external file to Base, something goes wrong with the pointers. When printed, they don't show a hexadecimal value and Julia crashes when I try to load a value from them.

When using normal array indexing, lights are green(thank you @staticfloat) and both normal and exponential pass the BigCrush tests. I have included the table creation as a test for the hardcoded tables.

[StefanKarpinski]

@inbounds doesn't avoid the null pointer checks when array objects could be undefined. Not sure if that's the issue, but that's a problem that I've encountered with this sort of thing.

[andreasnoack]

Correction: The pointers do show a hexadecimal number, i.e. zero, so for some reason they become null pointers when the code is loaded from Base, but not when loaded from an external file.

[JeffBezanson]

A Ptr is serialized as NULL, since saving the address would not do any good. On each startup, the exact addresses of objects can be different.

[JeffBezanson]

When the code is compiled and run in the same session, the JIT will be able to inline the pointer constants into the generated code, which should explain the speedup. Unfortunately there is no way to obtain this optimization in pre-compiled code. We could add the equivalent of @ccallable for data, allowing the pointers to be read via cglobal.

[andreasnoack]

@JeffBezanson, Thank you for the explanation. It makes sense.

I have removed the pointers and the normals are still much faster in the Julia implementation on julia.mit.edu, but the exponentials are slightly slower without the pointers.

The definitions with and without an RNG are almost identical. Is it possible to write a loop with @eval where I define functions with and without an argument? I have a loop tuple like

for t in (:nothing, :x)
    @eval begin
    f($t)=2
    end
end

but it doesn't define a version without an argument. Is it possible to make a definition like that?

[toivoh]

I believe that something like

for t in ([], [:x])
    @eval begin
    f($(t...))=2
    end
end

should do the trick.

[andreasnoack]

It works. Thanks @toivoh.

[ViralBShah]

Is this good to merge now?

[andreasnoack]

Yes, I think so.