[commons-rng] 03/03: Update to match ZigguratExponentialSampler

[ah...@apache.org]

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 <ah...@apache.org>
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);
     }