commons-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From er...@apache.org
Subject [1/3] commons-rng git commit: RNG-37: Review of "Ziggurat" implementation.
Date Thu, 04 Jan 2018 14:43:14 GMT
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 <erans@apache.org>
Authored: Thu Jan 4 12:29:11 2018 +0100
Committer: Gilles <erans@apache.org>
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</a> for sampling from a Gaussian
  * distribution with mean 0 and standard deviation 1.
  *
- * <p>
- * Implementation is based on this <a href="http://www.jstatsoft.org/article/view/v005i08/ziggurat.pdf">paper</a>
- * </p>
+ * The algorithm is explained in this
+ * <a href="http://www.jstatsoft.org/article/view/v005i08/ziggurat.pdf">paper</a>
+ * 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);
     }
 }


Mime
View raw message