From commits-return-60584-archive-asf-public=cust-asf.ponee.io@commons.apache.org Thu Jan 4 15:43:16 2018 Return-Path: X-Original-To: archive-asf-public@eu.ponee.io Delivered-To: archive-asf-public@eu.ponee.io Received: from cust-asf.ponee.io (cust-asf.ponee.io [163.172.22.183]) by mx-eu-01.ponee.io (Postfix) with ESMTP id ED164180657 for ; Thu, 4 Jan 2018 15:43:16 +0100 (CET) Received: by cust-asf.ponee.io (Postfix) id DD18B160C2B; Thu, 4 Jan 2018 14:43:16 +0000 (UTC) Delivered-To: archive-asf-public@cust-asf.ponee.io Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by cust-asf.ponee.io (Postfix) with SMTP id 07C57160C28 for ; Thu, 4 Jan 2018 15:43:15 +0100 (CET) Received: (qmail 21636 invoked by uid 500); 4 Jan 2018 14:43:15 -0000 Mailing-List: contact commits-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list commits@commons.apache.org Received: (qmail 21621 invoked by uid 99); 4 Jan 2018 14:43:15 -0000 Received: from git1-us-west.apache.org (HELO git1-us-west.apache.org) (140.211.11.23) by apache.org (qpsmtpd/0.29) with ESMTP; Thu, 04 Jan 2018 14:43:15 +0000 Received: by git1-us-west.apache.org (ASF Mail Server at git1-us-west.apache.org, from userid 33) id F3664DFF75; Thu, 4 Jan 2018 14:43:14 +0000 (UTC) Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit From: erans@apache.org To: commits@commons.apache.org Date: Thu, 04 Jan 2018 14:43:14 -0000 Message-Id: <439294d4b7ae40b6aa7d93902a34bfdf@git.apache.org> X-Mailer: ASF-Git Admin Mailer Subject: [1/3] commons-rng git commit: RNG-37: Review of "Ziggurat" implementation. Repository: commons-rng Updated Branches: refs/heads/master 3b3f19beb -> 8375c9ada RNG-37: Review of "Ziggurat" implementation. Avoid "magic" numbers. Adjust ranges: 64 bits ("long" and "double" here) vs 32 bits ("int" and "float" in original C code). Add missing calls to the RNG (macro in the original C code). Make branching explicit. Removed "n" suffix or prefix in variable and method names. Project: http://git-wip-us.apache.org/repos/asf/commons-rng/repo Commit: http://git-wip-us.apache.org/repos/asf/commons-rng/commit/75c7e43c Tree: http://git-wip-us.apache.org/repos/asf/commons-rng/tree/75c7e43c Diff: http://git-wip-us.apache.org/repos/asf/commons-rng/diff/75c7e43c Branch: refs/heads/master Commit: 75c7e43c7ea98ca710713c53d0506aad0da9cc42 Parents: 3b3f19b Author: Gilles Authored: Thu Jan 4 12:29:11 2018 +0100 Committer: Gilles Committed: Thu Jan 4 12:29:11 2018 +0100 ---------------------------------------------------------------------- .../rng/sampling/distribution/SamplerBase.java | 7 + .../ZigguratNormalizedGaussianSampler.java | 166 ++++++++++--------- 2 files changed, 94 insertions(+), 79 deletions(-) ---------------------------------------------------------------------- http://git-wip-us.apache.org/repos/asf/commons-rng/blob/75c7e43c/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/SamplerBase.java ---------------------------------------------------------------------- diff --git a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/SamplerBase.java b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/SamplerBase.java index 5972e71..eb74f75 100644 --- a/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/SamplerBase.java +++ b/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/SamplerBase.java @@ -55,6 +55,13 @@ public class SamplerBase { return rng.nextInt(max); } + /** + * @return a random {@code long} value. + */ + protected long nextLong() { + return rng.nextLong(); + } + /** {@inheritDoc} */ @Override public String toString() { http://git-wip-us.apache.org/repos/asf/commons-rng/blob/75c7e43c/commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/ZigguratNormalizedGaussianSampler.java ---------------------------------------------------------------------- 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 5a47b2a..7836bb3 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 @@ -24,64 +24,63 @@ import org.apache.commons.rng.UniformRandomProvider; * Marsaglia and Tsang "Ziggurat" method for sampling from a Gaussian * distribution with mean 0 and standard deviation 1. * - *

- * Implementation is based on this paper - *

+ * The algorithm is explained in this + * paper + * and this implementation has been adapted from the C code provided therein. * * @since 1.1 */ public class ZigguratNormalizedGaussianSampler extends SamplerBase implements NormalizedGaussianSampler { + /** Start of tail. */ + 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; + /** Auxiliary table. */ + private static final long[] K = new long[LEN]; + /** Auxiliary table. */ + private static final double[] W = new double[LEN]; + /** Auxiliary table. */ + private static final double[] F = new double[LEN]; - // Generates values from Gaussian (normal) probability distribution - // It uses two tables, integers KN and reals WN. Some 99% of the time, - // the required x is produced by: - // generate a random 32-bit integer j and let i be the index formed from - // the rightmost 8 bits of j. If j < k_i return x = j * w_i. + static { + // Filling the tables. - /** Auxiliary table (see code) */ - private static final int[] KN = new int[128]; - /** Auxiliary table (see code) */ - private static final double[] WN = new double[128]; - /** Auxiliary table (see code) */ - private static final double[] FN = new double[128]; + double d = R; + double t = d; + double fd = gauss(d); + final double q = V / fd; - static { - // Filling the tables. - // k_0 = 2^32 * r * f(dn) / vn - // k_i = 2^32 * ( x_{i-1} / x_i ) - // w_0 = 0.5^32 * vn / f(dn) - // w_i = 0.5^32 * x_i - // f(dn) = exp(-0.5 * dn * dn) - // where "dn" is the rightmost x_i - // "vn" is the area of the rectangle - final double m = 2147483648.0; // 2^31. - - // Provides z(r) = 0, where z(r): x_255 = r, vn = r * f(r) + integral_r^inf f(x) dx - final double vn = 9.91256303526217e-3; - - double dn = 3.442619855899; - double tn = dn; - double e = Math.exp(-0.5 * dn * dn); - final double q = vn / e; - - KN[0] = (int) ((dn / q) * m); - KN[1] = 0; - - WN[0] = q / m; - WN[127] = dn / m; - - FN[0] = 1; - FN[127] = e; - - for (int i = 126; i >= 1; i--){ - e = Math.exp(-0.5 * dn * dn); - dn = Math.sqrt(-2 * Math.log(vn / dn + e)); - KN[i+1] = (int) ((dn / tn) * m); - tn = dn; - FN[i] = e; - WN[i] = dn / m; + K[0] = (long) ((d / q) * MAX); + K[1] = 0; + + W[0] = q * ONE_OVER_MAX; + W[LAST] = d * ONE_OVER_MAX; + + 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); + + K[i + 1] = (long) ((d / t) * MAX); + t = d; + + F[i] = fd; + + W[i] = d * ONE_OVER_MAX; } } @@ -95,9 +94,19 @@ public class ZigguratNormalizedGaussianSampler /** {@inheritDoc} */ @Override public double sample() { - final int j = nextInt(); - final int i = j & 127; - return (j < KN[i]) ? j * WN[i] : nfix(j, i); + final long j = nextLong(); + final int i = (int) (j & LAST); + if (Math.abs(j) < K[i]) { + return j * W[i]; + } else { + return fix(j, i); + } + } + + /** {@inheritDoc} */ + @Override + public String toString() { + return "Ziggurat normalized Gaussian deviate [" + super.toString() + "]"; } /** @@ -107,43 +116,42 @@ public class ZigguratNormalizedGaussianSampler * @param iz Corresponding to hz cell's number. * @return the requested random value. */ - private double nfix(int hz, - int iz) { - // The start of the right tail. - final double r = 3.442619855899; - - double uni; + private double fix(long hz, + int iz) { double x; double y; while (true) { - uni = 0.5 + hz * 0.2328306e-9; - x = hz * WN[iz]; - // iz == 0 handles the base strip. + x = hz * W[iz]; if (iz == 0) { - // Return x from the tail. + // Base strip. do { - y = -Math.log(uni); - x = y * 0.2904764; - uni = 0.5 + nextInt() * 0.2328306e-9; + y = -Math.log(nextDouble()); + x = -Math.log(nextDouble()) * ONE_OVER_R; } while (y + y < x * x); - return (hz > 0) ? r + x : -r - x; - } - // iz > 0 handles the wedges of other strips. - if (FN[iz] + uni * (FN[iz - 1] - FN[iz]) < Math.exp(-0.5 * x * x)) { - return x; - } - hz = nextInt(); - iz = hz & 127; - if (Math.abs(hz) < KN[iz]) { - return hz * WN[iz]; + + final double out = R + x; + return hz > 0 ? out : -out; + } else { + // Wedge of other strips. + if (F[iz] + nextDouble() * (F[iz - 1] - F[iz]) < gauss(x)) { + return x; + } else { + hz = nextLong(); + iz = (int) (hz & LAST); + if (Math.abs(hz) < K[iz]) { + return hz * W[iz]; + } + } } } } - /** {@inheritDoc} */ - @Override - public String toString() { - return "Ziggurat normalized Gaussian deviate [" + super.toString() + "]"; + /** + * @param x Argument. + * @return \( e^{-\frac{x^2}{2}} \) + */ + private static double gauss(double x) { + return Math.exp(-0.5 * x * x); } }