From 0b5c8fc0dd21318d0254b1356bc0abad088e74dd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?H=C3=B6rmet=20Yiltiz?= <yiltiz@apple.com>
Date: Thu, 11 Jan 2024 01:19:11 -0500
Subject: [PATCH 1/3] Feature: add Gaussian sampler

---
 spork/randgen.janet | 38 ++++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)

diff --git a/spork/randgen.janet b/spork/randgen.janet
index 0c02130..eeee468 100644
--- a/spork/randgen.janet
+++ b/spork/randgen.janet
@@ -28,6 +28,44 @@
   (def diff (- end start))
   (+ start (math/floor (* diff (rand-uniform)))))
 
+(defn rand-gaussian
+  "Get a random sample from the standard Gaussian distribution. 
+  Optionall specify the mean m and the standard deviation sd. E.g.:
+  (rand-gaussian) # => 0.1324...
+  (rand-gaussian 5 0.1) # => 5.3397...
+  "
+  [&opt m sd]
+  (default [m sd] [0 1])
+  (defn scale [x] (+ m (* sd x)))
+
+  (def p (rand-uniform))
+
+  (def q (rand-uniform))
+
+  (math/rng-uniform ))
+
+  # We use the Box-Muller transform
+  (let [rho (math/sqrt (* -2 (math/log q)))
+        theta (* 2 math/pi p)
+        _box    (* rho (math/cos theta))
+        _muller (* rho (math/sin theta))
+        box    (scale _box)
+        muller (scale _muller)]
+
+    (yield box)
+    muller
+    ))
+
+(defn sample-n 
+  "Generate n samples based on the random sampler f. E.g.
+  (sample-n |(rand-int 0 3) 4) # => @[0 1 2 0]
+  (sample-n |(rand-uniform) 4)
+  (sample-n |(rand-gaussian 5 0.1) 4)
+  "
+  [f n]
+  (take n (generate [_ :iterate true] 
+            (f))))
+
 (defn rand-index
   "Get a random numeric index of an indexed data structure"
   [xs]

From da6d029145261e1548d59e4a5b8e7e711c3abf7d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?H=C3=B6rmet=20Yiltiz?= <yiltiz@apple.com>
Date: Thu, 11 Jan 2024 12:03:07 -0500
Subject: [PATCH 2/3] add tests

---
 spork/randgen.janet      | 14 +++++---------
 test/suite-randgen.janet | 20 ++++++++++++++++++++
 2 files changed, 25 insertions(+), 9 deletions(-)
 create mode 100644 test/suite-randgen.janet

diff --git a/spork/randgen.janet b/spork/randgen.janet
index eeee468..eb73294 100644
--- a/spork/randgen.janet
+++ b/spork/randgen.janet
@@ -35,14 +35,12 @@
   (rand-gaussian 5 0.1) # => 5.3397...
   "
   [&opt m sd]
-  (default [m sd] [0 1])
+  (default m 0)
+  (default sd 1)
   (defn scale [x] (+ m (* sd x)))
 
-  (def p (rand-uniform))
-
-  (def q (rand-uniform))
-
-  (math/rng-uniform ))
+  (def p (math/rng-uniform (get-rng)))
+  (def q (math/rng-uniform (get-rng)))
 
   # We use the Box-Muller transform
   (let [rho (math/sqrt (* -2 (math/log q)))
@@ -53,8 +51,7 @@
         muller (scale _muller)]
 
     (yield box)
-    muller
-    ))
+    muller))
 
 (defn sample-n 
   "Generate n samples based on the random sampler f. E.g.
@@ -123,4 +120,3 @@
   [weights & paths]
   ~(case (,rand-weights ,weights)
      ,;(array/concat @[] ;(map tuple (range (length paths)) paths))))
-
diff --git a/test/suite-randgen.janet b/test/suite-randgen.janet
new file mode 100644
index 0000000..b11e7f1
--- /dev/null
+++ b/test/suite-randgen.janet
@@ -0,0 +1,20 @@
+(use ../spork/test)
+(use ../spork/randgen)
+
+(start-suite)
+
+(assert-docs "/spork/randgen")
+
+(def delta 0.000000000000009)
+
+(assert (do
+          (set-seed 1)
+          (def expected @[5.00917440069385 
+                          5.0318756480189 
+                          4.97222830484943
+                          5.10069692026214])
+          (def actual
+            (sample-n |(rand-gaussian 5 0.1) 4))
+          (and (all |(> delta (math/abs (- $0 $1)))
+                    expected actual)))
+        "sample-rand-gaussian")

From 018bab4e09be19888f3922eba355b220cd6b6b1b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?H=C3=B6rmet=20Yiltiz?= <yiltiz@apple.com>
Date: Thu, 11 Jan 2024 12:06:52 -0500
Subject: [PATCH 3/3] reformat with spork/fmt/format-file

---
 spork/randgen.janet | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/spork/randgen.janet b/spork/randgen.janet
index eb73294..6d0fc9a 100644
--- a/spork/randgen.janet
+++ b/spork/randgen.janet
@@ -45,22 +45,22 @@
   # We use the Box-Muller transform
   (let [rho (math/sqrt (* -2 (math/log q)))
         theta (* 2 math/pi p)
-        _box    (* rho (math/cos theta))
+        _box (* rho (math/cos theta))
         _muller (* rho (math/sin theta))
-        box    (scale _box)
+        box (scale _box)
         muller (scale _muller)]
 
     (yield box)
     muller))
 
-(defn sample-n 
+(defn sample-n
   "Generate n samples based on the random sampler f. E.g.
   (sample-n |(rand-int 0 3) 4) # => @[0 1 2 0]
   (sample-n |(rand-uniform) 4)
   (sample-n |(rand-gaussian 5 0.1) 4)
   "
   [f n]
-  (take n (generate [_ :iterate true] 
+  (take n (generate [_ :iterate true]
             (f))))
 
 (defn rand-index