gaussr   xtlang

(bind-func gaussr
  (let ((phase 0))
    (lambda ()
      (let ((u1:double (random))
            (u2:double (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))
                (gaussr)
                (* v1 (sqrt (/ (* -2.0 (log s)) s))))
            (begin (set! phase (- 1 phase))
                   (* v2 (sqrt (/ (* -2.0 (log s)) s)))))))))