RNG-36: Ensure precondition. Not both elements of the generated pair can be zero.
Project: http://git-wip-us.apache.org/repos/asf/commons-rng/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-rng/commit/b1b44a5c Tree: http://git-wip-us.apache.org/repos/asf/commons-rng/tree/b1b44a5c Diff: http://git-wip-us.apache.org/repos/asf/commons-rng/diff/b1b44a5c Branch: refs/heads/master Commit: b1b44a5ce686c4b19b6c642fe38c24b9f33d5805 Parents: e591a65 Author: Gilles <er...@apache.org> Authored: Fri Apr 7 16:21:13 2017 +0200 Committer: Gilles <er...@apache.org> Committed: Fri Apr 7 16:21:13 2017 +0200 ---------------------------------------------------------------------- .../MarsagliaNormalizedGaussianSampler.java | 27 ++++++++++---------- 1 file changed, 13 insertions(+), 14 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-rng/blob/b1b44a5c/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaNormalizedGaussianSampler.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaNormalizedGaussianSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaNormalizedGaussianSampler.java index 6c86452..0ef6880 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaNormalizedGaussianSampler.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/MarsagliaNormalizedGaussianSampler.java @@ -43,7 +43,6 @@ public class MarsagliaNormalizedGaussianSampler /** {@inheritDoc} */ @Override public double sample() { - final double random; if (Double.isNaN(nextGaussian)) { // Rejection scheme for selecting a pair that lies within the unit circle. SAMPLE: while (true) { @@ -52,30 +51,30 @@ public class MarsagliaNormalizedGaussianSampler final double y = 2 * nextDouble() - 1; final double r2 = x * x + y * y; - if (r2 < 1) { - // Pair (x, y) is within unit circle. + if (r2 > 1 || r2 == 0) { + // Pair is not within the unit circle: Generate another one. + continue SAMPLE; + } - final double alpha = Math.sqrt(-2 * Math.log(r2) / r2); + // Pair (x, y) is within unit circle. + final double alpha = Math.sqrt(-2 * Math.log(r2) / r2); - // Return the first element of the generated pair. - random = alpha * x; + // Keep second element of the pair for next invocation. + nextGaussian = alpha * y; - // Keep second element of the pair for next invocation. - nextGaussian = alpha * y; - break SAMPLE; - } - // Pair is not within the unit circle: Generate another one. + // Return the first element of the generated pair. + return alpha * x; } } else { // Use the second element of the pair (generated at the // previous invocation). - random = nextGaussian; + final double r = nextGaussian; // Both elements of the pair have been used. nextGaussian = Double.NaN; - } - return random; + return r; + } } /** {@inheritDoc} */
