(defparameter *last-one* 0d0)
(declaim (type double-float *last-one*))
(defparameter *has-last-one* nil)
(declaim (type boolean *has-last-one*))

(declaim (ftype (function () double-float) normal-distribution-1))
(defun normal-distribution-1 ()
  "Slow old fashioned normal distribution."
  (declare (optimize (speed 3) (safety 0)))
  (if *has-last-one*
      (progn
	(setf *has-last-one* nil)
	*last-one*)
      (let ((x1 0d0) (x2 0d0) (w 0d0) (y1 0d0))
	(declare (type double-float x1 x2 w y1))
	(loop :do
	   (setf x1 (* 2.0 (- (random 1d0) 1d0))
		 x2 (* 2.0 (- (random 1d0) 1d0))
		 w (+ (* x1 x1) (* x2 x2)))
	   :while (or (> w 1) (= w 0)))
	(setf w (sqrt (/ (* -2d0 (log w)) w))
	      y1 (* x1 w)
	      *last-one* (* x2 w)
	      *has-last-one* t)
	(the double-float y1))))

(defmacro with-slow-normal-dist (() &body body)
  "Wrap around calls of normal-distribution-1 to be thread “safe”."
  `(let ((*last-one* 0d0)
	 (*has-last-one* nil))
     ,@body))

;; Stolen from:
;; https://github.com/tpapp/cl-random/blob/master/src/univariate.lisp
;; Which doesn't seem to have a license.

(defun normal-distribution-2 ()
  "Draw a random number from N(0,1)."
  ;; Method from Leva (1992). This is considered much better/faster than the
  ;; Box-Muller method.
  (declare (optimize (speed 3) (safety 0))
	   #+sbcl (sb-ext:muffle-conditions sb-ext:compiler-note))
  (tagbody
   top
     (let* ((u (random 1d0))
	    (v (* 1.7156d0 (- (random 1d0) 0.5d0)))
	    (x (- u 0.449871d0))
	    (y (+ (abs v) 0.386595d0))
	    (q (+ (expt x 2) (* y (- (* 0.19600d0 y) (* 0.25472d0 x))))))
       (if (and (> q 0.27597d0)
		(or (> q 0.27846d0)
		    (plusp (+ (expt v 2) (* 4 (expt u 2) (log u))))))
	   (go top)
	   (return-from normal-distribution-2 (/ v u))))))