;; gaussian random number X ~ N(0,1)
;; algorithm by Marsaglia http://c-faq.com/lib/gaussian.html
(bind-func dsp_randn
(let ((phase:i64 0))
(lambda ()
(let ((u1:SAMPLE (random))
(u2:SAMPLE (random))
(v1 (- (* 2.0 u1) 1.0))
(v2 (- (* 2.0 u2) 1.0))
(s (+ (* v1 v1) (* v2 v2))))
(if (= phase 0)
(if (or (> s 1.0) (= s 0.0) (= s 1.0))
(dsp_randn)
(* v1 (sqrt (/ (* -2.0 (log s)) s))))
(begin (set! phase (- 1 phase))
(* v2 (sqrt (/ (* -2.0 (log s)) s)))))))))