seems to work even if pretty damn slow

search.pl

@@ -458,10 +458,6 @@ normal_sample([P | L], Rho) :-
   normal_sample(L, Rho).
 
 
-normal_random(X) :-
-  random(X).
-
-
 generate_random_seed :-
   (
     toplevel:random_seed(Seed)

util.pl

@@ -53,7 +53,8 @@
     write_expression_to_c/2,
     write_expression_to_c_with_parentheses/2,
     write_condition_to_c/2,
-    write_condition_to_c_with_parentheses/2
+    write_condition_to_c_with_parentheses/2,
+    normal_random/1
   ]).
 
 :- use_module(library(process)).
@@ -885,3 +886,20 @@ write_condition_to_c_with_parentheses(Condition, CallBack) :-
   write('('),
   write_condition_to_c(Condition, CallBack),
   write(')').
+
+
+% Box-Muller method for random sampling a normal distribution
+:- dynamic(done_normal/1).
+
+
+normal_random(X) :-
+  retract(done_normal(X)),
+  !.
+
+normal_random(X) :-
+  random(U),
+  random(V),
+  LN is sqrt(- 2 * log(U)),
+  X is LN * cos(2*pi*V),
+  Y is LN * sin(2*pi*V),
+  assertz(done_normal(Y)).