Harvey> (defun one-normal-rand ()
Harvey> (let* ((u (random 1d0))
Harvey> (u1 (random 1d0))
Harvey> (v (* SQRT2OVERE (- (* 2d0 u1) 1d0)))
Harvey> (x (/ v u))
Harvey> (y (/ (* x x) 4d0)))
Harvey> (if (and (> y (- 1d0 u))
Harvey> (> y (- (log u))))
Harvey> (one-normal-rand)
Harvey> x)))
I think random numbers is quite expensive; generating the gaussians
pairwise gives an efficiency of pi/4:
(defun gausspair (&optional (state *random-state*))
(declare (optimize (speed 3) (safety 0)))
(loop
for u of-type double-float = (- (random 2.d0 state) 1.d0)
for v of-type double-float = (- (random 2.d0 state) 1.d0)
for r of-type double-float = (+ (expt u 2) (expt v 2))
until (< r 1.d0)
finally (let* ((l (log r))
(d (sqrt (/ (* -2.d0 l) r))))
(declare (type double-float l d))
(return (values (* u d) (* v d))))))
Ole
