Improved Gaussian random number generation from ThreadLocalRandom in methods of RandomVariates class

CHANGELOG.md

@@ -4,11 +4,12 @@ All notable changes to this project will be documented in this file.
 The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
-## [Unreleased] - 2024-05-27
+## [Unreleased] - 2024-05-30
 
 ### Added
 
 ### Changed
+* Improved Gaussian random number generation from ThreadLocalRandom in methods of RandomVariates class (internal change).
 
 ### Deprecated
 

src/main/java/org/cicirello/math/rand/RandomVariates.java

@@ -1,6 +1,6 @@
 /*
  * rho mu - A Java library of randomization enhancements and other math utilities.
- * Copyright (C) 2017-2023 Vincent A. Cicirello, <https://www.cicirello.org/>.
+ * Copyright (C) 2017-2024 Vincent A. Cicirello, <https://www.cicirello.org/>.
  *
  * This file is part of the rho mu library.
  *
@@ -120,17 +120,24 @@ public static double nextCauchy(double scale, RandomGenerator r) {
 
   /**
    * Generates a random number from a Gaussian distribution with mean mu and standard deviation
-   * sigma. This method uses the library's current most-efficient algorithm for Gaussian random
-   * number generation, which may change in future releases. If you require a guarantee of the
-   * algorithm used, then see the API for the classes that implement specific Gaussian algorithms.
-   * {@link ThreadLocalRandom} is used as the source of randomness.
+   * sigma.
+   *
+   * <p>{@link ThreadLocalRandom} is used as the source of randomness. However, it does not directly
+   * use ThreadLocalRandom's nextGaussian, which is the slow polar method. Instead, it uses the
+   * approach described in the following paper to use the RandomGenerator interface's default
+   * implementation of McFarland's modified ziggurat algorithm, which is much faster.
+   *
+   * <p>Vincent A. Cicirello. 2024. <a href="https://reports.cicirello.org/24/009/">Fast Gaussian
+   * Distributed Pseudorandom Number Generation in Java via the Ziggurat Algorithm</a>. Technical
+   * Report ALG-24-009, Cicirello.org, May 2024. <a
+   * href="https://reports.cicirello.org/24/009/ALG-24-009.pdf">[PDF]</a>
    *
    * @param mu The mean of the Gaussian.
    * @param sigma The standard deviation of the Gaussian.
    * @return A random number from a Gaussian distribution with mean mu and standard deviation sigma.
    */
   public static double nextGaussian(double mu, double sigma) {
-    return mu + ZigguratGaussian.nextGaussian(sigma);
+    return gaussianThreadLocalRandom.nextGaussian(mu, sigma);
   }
 
   /**
@@ -145,21 +152,28 @@ public static double nextGaussian(double mu, double sigma) {
    * @return A random number from a Gaussian distribution with mean mu and standard deviation sigma.
    */
   public static double nextGaussian(double mu, double sigma, RandomGenerator r) {
-    return mu + ZigguratGaussian.nextGaussian(sigma, r);
+    return mu + sigma * ZigguratGaussian.nextGaussian(r);
   }
 
   /**
    * Generates a random number from a Gaussian distribution with mean 0 and standard deviation
-   * sigma. This method uses the library's current most-efficient algorithm for Gaussian random
-   * number generation, which may change in future releases. If you require a guarantee of the
-   * algorithm used, then see the API for the classes that implement specific Gaussian algorithms.
-   * {@link ThreadLocalRandom} is used as the source of randomness.
+   * sigma.
+   *
+   * <p>{@link ThreadLocalRandom} is used as the source of randomness. However, it does not directly
+   * use ThreadLocalRandom's nextGaussian, which is the slow polar method. Instead, it uses the
+   * approach described in the following paper to use the RandomGenerator interface's default
+   * implementation of McFarland's modified ziggurat algorithm, which is much faster.
+   *
+   * <p>Vincent A. Cicirello. 2024. <a href="https://reports.cicirello.org/24/009/">Fast Gaussian
+   * Distributed Pseudorandom Number Generation in Java via the Ziggurat Algorithm</a>. Technical
+   * Report ALG-24-009, Cicirello.org, May 2024. <a
+   * href="https://reports.cicirello.org/24/009/ALG-24-009.pdf">[PDF]</a>
    *
    * @param sigma The standard deviation of the Gaussian.
    * @return A random number from a Gaussian distribution with mean 0 and standard deviation sigma.
    */
   public static double nextGaussian(double sigma) {
-    return ZigguratGaussian.nextGaussian(sigma);
+    return gaussianThreadLocalRandom.nextGaussian(0, sigma);
   }
 
   /**
@@ -173,20 +187,26 @@ public static double nextGaussian(double sigma) {
    * @return A random number from a Gaussian distribution with mean 0 and standard deviation sigma.
    */
   public static double nextGaussian(double sigma, RandomGenerator r) {
-    return ZigguratGaussian.nextGaussian(sigma, r);
+    return sigma * ZigguratGaussian.nextGaussian(r);
   }
 
   /**
    * Generates a random number from a Gaussian distribution with mean 0 and standard deviation 1.
-   * This method uses the library's current most-efficient algorithm for Gaussian random number
-   * generation, which may change in future releases. If you require a guarantee of the algorithm
-   * used, then see the API for the classes that implement specific Gaussian algorithms. {@link
-   * ThreadLocalRandom} is used as the pseudorandom number generator for the source of randomness.
+   *
+   * <p>{@link ThreadLocalRandom} is used as the source of randomness. However, it does not directly
+   * use ThreadLocalRandom's nextGaussian, which is the slow polar method. Instead, it uses the
+   * approach described in the following paper to use the RandomGenerator interface's default
+   * implementation of McFarland's modified ziggurat algorithm, which is much faster.
+   *
+   * <p>Vincent A. Cicirello. 2024. <a href="https://reports.cicirello.org/24/009/">Fast Gaussian
+   * Distributed Pseudorandom Number Generation in Java via the Ziggurat Algorithm</a>. Technical
+   * Report ALG-24-009, Cicirello.org, May 2024. <a
+   * href="https://reports.cicirello.org/24/009/ALG-24-009.pdf">[PDF]</a>
    *
    * @return A random number from a Gaussian distribution with mean 0 and standard deviation 1.
    */
   public static double nextGaussian() {
-    return ZigguratGaussian.nextGaussian();
+    return gaussianThreadLocalRandom.nextGaussian();
   }
 
   /**
@@ -202,6 +222,31 @@ public static double nextGaussian(RandomGenerator r) {
     return ZigguratGaussian.nextGaussian(r);
   }
 
+  /*
+   * For implementing nextGaussian with ThreadLocalRandom:
+   *
+   * Beginning in Java 17, the RandomGenerator interface provides a default
+   * implementation of nextGaussian that uses the very fast modified ziggurat
+   * algorithm. However, ThreadLocalRandom extends the legacy Random class, which
+   * instead uses the very slow polar method because the Random class's linear
+   * congruential method doesn't meet the Java 17's modified ziggurat's requirements
+   * on the quality of the randomness of low-order bits of random longs. The
+   * ThreadLocalRandom actually does meet those requirements because it uses a
+   * stronger pseudorandom number generator. This static field provides a minimal
+   * implementation of the RandomGenerator interface, using ThreadLocalRandom to
+   * generate random longs, and otherwise inheriting the default implementations
+   * of all other interface methods, including nextGaussian. Use the nextGaussian
+   * method on this instance in the methods that generate Gaussians with
+   * ThreadLocalRandom.
+   */
+  private static final RandomGenerator gaussianThreadLocalRandom =
+      new RandomGenerator() {
+        @Override
+        public long nextLong() {
+          return ThreadLocalRandom.current().nextLong();
+        }
+      };
+
   /*
    * INTERNAL METHODS FOR nextCauchy start here.
    *

src/main/java/org/cicirello/math/rand/ZigguratGaussian.java

@@ -49,7 +49,6 @@
 
 package org.cicirello.math.rand;
 
-import java.util.concurrent.ThreadLocalRandom;
 import java.util.random.RandomGenerator;
 
 /**
@@ -62,12 +61,22 @@
  * of the class. If interested, see the source code comments, which highlights any differences
  * between this Java implementation and the C implementation on which it is based.
  *
+ * <p>See the following report for more information about this implementation along with
+ * experimental results:
+ *
+ * <ul>
+ *   <li>Vincent A. Cicirello. 2024. <a href="https://reports.cicirello.org/24/009/">Fast Gaussian
+ *       Distributed Pseudorandom Number Generation in Java via the Ziggurat Algorithm</a>.
+ *       Technical Report ALG-24-009, Cicirello.org, May 2024. <a
+ *       href="https://reports.cicirello.org/24/009/ALG-24-009.pdf">[PDF]</a>
+ * </ul>
+ *
  * <p>This Java implementation originated as part of an effort to speed up the runtime of a parallel
  * genetic algorithm (PGA). The PGA in question evolved its control parameters (i.e., crossover and
  * mutation rates, etc) using Gaussian mutation. The only Gaussian implementation within the Java
  * API, at that time, was the polar method (nextGaussian method of the {@link java.util.Random
- * Random} and {@link ThreadLocalRandom} classes, however the polar method is quite slow relative to
- * other newer available alternatives, such as the Ziggurat method.
+ * Random} class, however the polar method is quite slow relative to other newer available
+ * alternatives, such as the Ziggurat method.
  *
  * <p>You can find some experimental data comparing the performance of a sequential genetic
  * algorithm (GA) using this implementation of the Ziggurat method for Gaussian mutation vs using
@@ -104,25 +113,18 @@ final class ZigguratGaussian {
    * A Few Java Implementation Notes:
    *
    * - I wanted to support generating Gaussian distributed random
-   *   numbers via the Ziggurat method using any of the following
-   *   Java classes as the underlying pseudorandom number generator:
-   *   Random, ThreadLocalRandom, and SplittableRandom. Aside from SecureRandom,
-   *   these were all of the random number generator classes available
-   *   in the Java API at the time this was originally implemented for Java 8.
+   *   numbers via the Ziggurat method using any of the Java 8
+   *   pseudorandom number generator classes.
    *
-   * - One option to accomplish that would be to have three classes,
-   *   each extending one of those three and either overriding
-   *   nextGaussian (for the two classes that have such a method),
-   *   or in the case of SplittableRandom adding such a method.
-   *   There are problems with this.  First, SplittableRandom is a
-   *   final class so extending is not an option.  Second, although
-   *   ThreadLocalRandom is not final, it is impractical to extend
-   *   none-the-less due to how it manages each local thread's
+   * - One option considered was extending each of those classes and
+   *   overriding nextGaussian. However, SplittableRandom is a final
+   *   class, and although ThreadLocalRandom is not final, it is
+   *   impractical to extend due to how it manages each local thread's
    *   ThreadLocalRandom instance, relying on the static method current().
    *   Overriding ThreadLocalRandom would likely involve recreating most
    *   of the class's functionality due to this.  So although this option
-   *   would be fine for the Random class, I desired a common approach
-   *   to support all three.
+   *   would be fine for the Random class and SecureRandom class, I desired
+   *   a common approach to support all three.
    *
    * - A second option would be to have three classes, one corresponding to
    *   each of the three Java library classes above.  In the cases of
@@ -162,6 +164,30 @@ final class ZigguratGaussian {
    *   by using this interface. That is, rather than the previous one method that takes
    *   a Random and another method that takes a SplittableRandom, we now have a single
    *   method that takes a RandomGenerator.
+   *
+   * COMMENT UPDATE (5/30/2024):
+   * - Also as of Java 17, McFarland's modified ziggurat algorithm, which is faster than
+   *   my port of the original ziggurat from C, is available as the default implementation
+   *   in the RandomGenerator interface. However, it has requirements that the
+   *   Random class does not meet, so the Random class in Java 17 continues to use the
+   *   slow polar method. SecureRandom and ThreadLocalRandom both extend Random and as a
+   *   consequence inherit the slow polar method, despite both meeting requirements of
+   *   the modified ziggurat. So, my implementation of the original ziggurat still has
+   *   a meaningful purpose for SecureRandom, as well as for Random (the original ziggurat)
+   *   does not have the same statistical requirements of the underlying pseudorandom
+   *   number generator. However, there is a way of gaining access to the Java 17
+   *   modified ziggurat for the ThreadLocalRandom class, which is now done in the
+   *   RandomVariates class. See the following report for more detailed explanations of
+   *   the relevant issues:
+   *
+   *   Vincent A. Cicirello. 2024. <a href="https://reports.cicirello.org/24/009/">Fast Gaussian
+   *   Distributed Pseudorandom Number Generation in Java via the Ziggurat Algorithm</a>.
+   *   Technical Report ALG-24-009, Cicirello.org, May 2024. <a
+   *   href="https://reports.cicirello.org/24/009/ALG-24-009.pdf">[PDF]</a>
+   *
+   * - Thus, as of 5/30/2024, we are stripping out all of the ThreadLocalRandom related
+   *   parts of this class as they are no longer needed. Since this class is an internal
+   *   class, this is not a breaking change.
    */
 
   /*
@@ -391,48 +417,13 @@ private ZigguratGaussian() {}
     1.83813550477e-07, 1.92166040885e-07, 2.05295471952e-07, 2.22600839893e-07
   };
 
-  /**
-   * Generates a random number from a Gaussian distribution with mean 0 and standard deviation,
-   * sigma, of your choosing. {@link ThreadLocalRandom} is used as the pseudorandom number generator
-   * for the source of randomness.
-   *
-   * @param sigma The standard deviation of the Gaussian.
-   * @return A random number from a Gaussian distribution with mean 0 and standard deviation sigma.
-   */
-  public static double nextGaussian(double sigma) {
-    return sigma * nextGaussian(ThreadLocalRandom.current());
-  }
-
-  /**
-   * Generates a random number from a Gaussian distribution with mean 0 and standard deviation,
-   * sigma, of your choosing.
-   *
-   * @param sigma The standard deviation of the Gaussian.
-   * @param r The pseudorandom number generator to use for the source of randomness.
-   * @return A random number from a Gaussian distribution with mean 0 and standard deviation sigma.
-   */
-  public static double nextGaussian(double sigma, RandomGenerator r) {
-    return sigma * nextGaussian(r);
-  }
-
-  /**
-   * Generates a random number from a Gaussian distribution with mean 0 and standard deviation 1.
-   * {@link ThreadLocalRandom} is used as the pseudorandom number generator for the source of
-   * randomness.
-   *
-   * @return A random number from a Gaussian distribution with mean 0 and standard deviation 1.
-   */
-  public static double nextGaussian() {
-    return nextGaussian(ThreadLocalRandom.current());
-  }
-
   /**
    * Generates a random number from a Gaussian distribution with mean 0 and standard deviation 1.
    *
    * @param r The pseudorandom number generator to use for the source of randomness.
    * @return A random number from a Gaussian distribution with mean 0 and standard deviation 1.
    */
-  public static double nextGaussian(RandomGenerator r) {
+  static double nextGaussian(RandomGenerator r) {
     double x, y;
     int sign;
 

src/test/java/org/cicirello/math/rand/DefaultGaussianTests.java → src/test/java/org/cicirello/math/rand/GaussianTests.java

@@ -1,6 +1,6 @@
 /*
  * rho mu - A Java library of randomization enhancements and other math utilities.
- * Copyright 2017-2022 Vincent A. Cicirello, <https://www.cicirello.org/>.
+ * Copyright 2017-2024 Vincent A. Cicirello, <https://www.cicirello.org/>.
  *
  * This file is part of the rho mu library.
  *
@@ -29,8 +29,8 @@
 import java.util.SplittableRandom;
 import org.junit.jupiter.api.*;
 
-/** JUnit test cases for the methods of the library's default Gaussian class. */
-public class DefaultGaussianTests {
+/** JUnit test cases for the nextGaussian methods of RandomVariates. */
+public class GaussianTests {
 
   // Test cases use chi square goodness of fit.  This constant
   // can be used to adjust the number of samples used for this test.