Unverified Commit b8bd1033 authored by aarvid's avatar aarvid Committed by Attila Lendvai

Fixed bug in gaussian-random when min or max is nil but not both.

parent 1f3ae3aa
...@@ -18,31 +18,35 @@ normal distribution around 0.0d0. ...@@ -18,31 +18,35 @@ normal distribution around 0.0d0.
Sufficiently positive MIN or negative MAX will cause the algorithm used to Sufficiently positive MIN or negative MAX will cause the algorithm used to
take a very long time. If MIN is positive it should be close to zero, and take a very long time. If MIN is positive it should be close to zero, and
similarly if MAX is negative it should be close to zero." similarly if MAX is negative it should be close to zero."
(labels ((gauss () (macrolet
(loop ((valid (x)
for x1 = (- (random 2.0d0) 1.0d0) `(<= (or min ,x) ,x (or max ,x)) ))
for x2 = (- (random 2.0d0) 1.0d0) (labels
for w = (+ (expt x1 2) (expt x2 2)) ((gauss ()
when (< w 1.0d0) (loop
do (let ((v (sqrt (/ (* -2.0d0 (log w)) w)))) for x1 = (- (random 2.0d0) 1.0d0)
(return (values (* x1 v) (* x2 v)))))) for x2 = (- (random 2.0d0) 1.0d0)
(guard (x min max) for w = (+ (expt x1 2) (expt x2 2))
(unless (<= min x max) when (< w 1.0d0)
(tagbody do (let ((v (sqrt (/ (* -2.0d0 (log w)) w))))
:retry (return (values (* x1 v) (* x2 v))))))
(multiple-value-bind (x1 x2) (gauss) (guard (x)
(when (<= min x1 max) (unless (valid x)
(setf x x1) (tagbody
(go :done)) :retry
(when (<= min x2 max) (multiple-value-bind (x1 x2) (gauss)
(setf x x2) (when (valid x1)
(go :done)) (setf x x1)
(go :retry)) (go :done))
:done)) (when (valid x2)
x)) (setf x x2)
(multiple-value-bind (g1 g2) (gauss) (go :done))
(values (guard g1 (or min g1) (or max g1)) (go :retry))
(guard g2 (or min g2) (or max g2)))))) :done))
x))
(multiple-value-bind
(g1 g2) (gauss)
(values (guard g1) (guard g2))))))
(declaim (inline iota)) (declaim (inline iota))
(defun iota (n &key (start 0) (step 1)) (defun iota (n &key (start 0) (step 1))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment