This is an automated email from the ASF dual-hosted git repository. aherbert pushed a commit to branch master in repository https://gitbox.apache.org/repos/asf/commons-rng.git
commit a1a95bf66335df6679c92d4753eefd928afd809b Author: aherbert <[email protected]> AuthorDate: Sun Jul 4 08:54:04 2021 +0100 Update to match ZigguratExponentialSampler - Remove unnecessary nested statements - Move fill table code inside the static block - Update the branch frequencies - Rename gauss to pdf - Primitive cast before the mask is applied to obtain the index --- .../ZigguratNormalizedGaussianSampler.java | 73 +++++++++++----------- 1 file changed, 37 insertions(+), 36 deletions(-) diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.java index a8b3e5d..7b62be8 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.java @@ -43,52 +43,54 @@ public class ZigguratNormalizedGaussianSampler private static final double R = 3.442619855899; /** Inverse of R. */ private static final double ONE_OVER_R = 1 / R; - /** Rectangle area. */ - private static final double V = 9.91256303526217e-3; - /** 2^63. */ - private static final double MAX = Math.pow(2, 63); - /** 2^-63. */ - private static final double ONE_OVER_MAX = 1d / MAX; - /** Number of entries. */ - private static final int LEN = 128; - /** Index of last entry. */ - private static final int LAST = LEN - 1; + /** Index of last entry in the tables (which have a size that is a power of 2). */ + private static final int LAST = 127; /** Auxiliary table. */ - private static final long[] K = new long[LEN]; + private static final long[] K; /** Auxiliary table. */ - private static final double[] W = new double[LEN]; + private static final double[] W; /** Auxiliary table. */ - private static final double[] F = new double[LEN]; + private static final double[] F; + /** Underlying source of randomness. */ private final UniformRandomProvider rng; static { // Filling the tables. + // Rectangle area. + final double v = 9.91256303526217e-3; + // Direction support uses the sign bit so the maximum magnitude from the long is 2^63 + final double max = Math.pow(2, 63); + final double oneOverMax = 1d / max; + + K = new long[LAST + 1]; + W = new double[LAST + 1]; + F = new double[LAST + 1]; double d = R; double t = d; - double fd = gauss(d); - final double q = V / fd; + double fd = pdf(d); + final double q = v / fd; - K[0] = (long) ((d / q) * MAX); + K[0] = (long) ((d / q) * max); K[1] = 0; - W[0] = q * ONE_OVER_MAX; - W[LAST] = d * ONE_OVER_MAX; + W[0] = q * oneOverMax; + W[LAST] = d * oneOverMax; F[0] = 1; F[LAST] = fd; for (int i = LAST - 1; i >= 1; i--) { - d = Math.sqrt(-2 * Math.log(V / d + fd)); - fd = gauss(d); + d = Math.sqrt(-2 * Math.log(v / d + fd)); + fd = pdf(d); - K[i + 1] = (long) ((d / t) * MAX); + K[i + 1] = (long) ((d / t) * max); t = d; F[i] = fd; - W[i] = d * ONE_OVER_MAX; + W[i] = d * oneOverMax; } } @@ -103,12 +105,11 @@ public class ZigguratNormalizedGaussianSampler @Override public double sample() { final long j = rng.nextLong(); - final int i = (int) (j & LAST); + final int i = ((int) j) & LAST; if (Math.abs(j) < K[i]) { return j * W[i]; - } else { - return fix(j, i); } + return fix(j, i); } /** {@inheritDoc} */ @@ -140,24 +141,24 @@ public class ZigguratNormalizedGaussianSampler final double out = R + x; return hz > 0 ? out : -out; - } else { - // Wedge of other strips. - // This branch is called about 0.027323 times per sample. - if (F[iz] + rng.nextDouble() * (F[iz - 1] - F[iz]) < gauss(x)) { - return x; - } else { - // Try again. - // This branch is called about 0.012362 times per sample. - return sample(); - } } + // Wedge of other strips. + // This branch is called about 0.027323 times per sample. + if (F[iz] + rng.nextDouble() * (F[iz - 1] - F[iz]) < pdf(x)) { + return x; + } + // Try again. + // This branch is called about 0.012362 times per sample. + return sample(); } /** + * Compute the Gaussian probability density function {@code f(x) = e^-0.5x^2}. + * * @param x Argument. * @return \( e^{-\frac{x^2}{2}} \) */ - private static double gauss(double x) { + private static double pdf(double x) { return Math.exp(-0.5 * x * x); }
