You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by ah...@apache.org on 2021/08/24 17:13:37 UTC
[commons-rng] 04/05: Add ziggurat exponential sampler with no
recursion to the benchmark
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 a39ce5d8317fdd87bca9d80710b2774bd636f1b8
Author: aherbert <ah...@apache.org>
AuthorDate: Tue Aug 24 17:42:55 2021 +0100
Add ziggurat exponential sampler with no recursion to the benchmark
---
.../distribution/ZigguratSamplerPerformance.java | 143 ++++++++++++++++++++-
1 file changed, 140 insertions(+), 3 deletions(-)
diff --git a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ZigguratSamplerPerformance.java b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ZigguratSamplerPerformance.java
index 36ba206..ed1531c 100644
--- a/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ZigguratSamplerPerformance.java
+++ b/commons-rng-examples/examples-jmh/src/main/java/org/apache/commons/rng/examples/jmh/sampling/distribution/ZigguratSamplerPerformance.java
@@ -95,6 +95,8 @@ public class ZigguratSamplerPerformance {
private static final String MOD_EXPONENTIAL_INLINING = "ModExponentialInlining";
/** The name for the {@link ModifiedZigguratExponentialSamplerLoop}. */
private static final String MOD_EXPONENTIAL_LOOP = "ModExponentialLoop";
+ /** The name for the {@link ModifiedZigguratExponentialSamplerLoop2}. */
+ private static final String MOD_EXPONENTIAL_LOOP2 = "ModExponentialLoop2";
/** The name for the {@link ModifiedZigguratExponentialSamplerRecursion}. */
private static final String MOD_EXPONENTIAL_RECURSION = "ModExponentialRecursion";
/** The name for the {@link ModifiedZigguratExponentialSamplerIntMap}. */
@@ -214,7 +216,8 @@ public class ZigguratSamplerPerformance {
MOD_GAUSSIAN_INLINING_SIMPLE_OVERHANGS, MOD_GAUSSIAN_INT_MAP, MOD_GAUSSIAN_512,
// Experimental McFarland Gaussian ziggurat samplers
MOD_EXPONENTIAL2, MOD_EXPONENTIAL_SIMPLE_OVERHANGS, MOD_EXPONENTIAL_INLINING,
- MOD_EXPONENTIAL_LOOP, MOD_EXPONENTIAL_RECURSION, MOD_EXPONENTIAL_INT_MAP, MOD_EXPONENTIAL_512})
+ MOD_EXPONENTIAL_LOOP, MOD_EXPONENTIAL_LOOP2,
+ MOD_EXPONENTIAL_RECURSION, MOD_EXPONENTIAL_INT_MAP, MOD_EXPONENTIAL_512})
private String type;
/** The sampler. */
@@ -273,6 +276,8 @@ public class ZigguratSamplerPerformance {
return new ModifiedZigguratExponentialSamplerInlining(rng);
} else if (MOD_EXPONENTIAL_LOOP.equals(type)) {
return new ModifiedZigguratExponentialSamplerLoop(rng);
+ } else if (MOD_EXPONENTIAL_LOOP2.equals(type)) {
+ return new ModifiedZigguratExponentialSamplerLoop2(rng);
} else if (MOD_EXPONENTIAL_RECURSION.equals(type)) {
return new ModifiedZigguratExponentialSamplerRecursion(rng);
} else if (MOD_EXPONENTIAL_INT_MAP.equals(type)) {
@@ -324,8 +329,8 @@ public class ZigguratSamplerPerformance {
MOD_GAUSSIAN2, MOD_GAUSSIAN_SIMPLE_OVERHANGS, MOD_GAUSSIAN_INLINING,
MOD_GAUSSIAN_INLINING_SIMPLE_OVERHANGS, MOD_GAUSSIAN_INT_MAP, MOD_GAUSSIAN_512,
// Experimental McFarland Gaussian ziggurat samplers
- MOD_EXPONENTIAL2, MOD_EXPONENTIAL_SIMPLE_OVERHANGS, MOD_EXPONENTIAL_INLINING,
- MOD_EXPONENTIAL_LOOP, MOD_EXPONENTIAL_RECURSION, MOD_EXPONENTIAL_INT_MAP, MOD_EXPONENTIAL_512})
+ MOD_EXPONENTIAL_LOOP, MOD_EXPONENTIAL_LOOP2,
+ MOD_EXPONENTIAL_RECURSION, MOD_EXPONENTIAL_INT_MAP, MOD_EXPONENTIAL_512})
private String type;
/** The size. */
@@ -2747,6 +2752,138 @@ public class ZigguratSamplerPerformance {
*
* <p>Uses the algorithm from McFarland, C.D. (2016).
*
+ * <p>This implementation separates sampling of the main ziggurat and sampling from the edge
+ * into different methods. This allows inlining of the main sample method.
+ *
+ * <p>The sampler will output different values due to the use of a bit shift to generate
+ * unsigned integers. This removes the requirement to load the mask MAX_INT64
+ * and ensures the method is under 35 bytes.
+ *
+ * <p>Tail sampling outside of the main sample method is performed in a loop. No recursion
+ * is used. The first random deviate is recycled if the sample if from the edge.
+ */
+ static class ModifiedZigguratExponentialSamplerLoop2
+ extends ModifiedZigguratExponentialSampler {
+
+ /**
+ * @param rng Generator of uniformly distributed random numbers.
+ */
+ ModifiedZigguratExponentialSamplerLoop2(UniformRandomProvider rng) {
+ super(rng);
+ }
+
+ /** {@inheritDoc} */
+ @Override
+ public double sample() {
+ // Ideally this method byte code size should be below -XX:MaxInlineSize
+ // (which defaults to 35 bytes). This compiles to 35 bytes.
+
+ final long x = nextLong();
+ // Float multiplication squashes these last 8 bits, so they can be used to sample i
+ final int i = ((int) x) & 0xff;
+
+ if (i < I_MAX) {
+ // Early exit.
+ // This branch is called about 0.984379 times per call into createSample.
+ // Note: Frequencies have been empirically measured for the first call to
+ // createSample; recursion due to retries have been ignored. Frequencies sum to 1.
+ return X[i] * (x >>> 1);
+ }
+
+ // Recycle x as the upper 56 bits have not been used.
+
+ return edgeSample(x);
+ }
+
+ /**
+ * Create the sample from the edge of the ziggurat.
+ *
+ * <p>This method has been extracted to fit the main sample method within 35 bytes (the
+ * default size for a JVM to inline a method).
+ *
+ * @param xx Initial random deviate
+ * @return a sample
+ */
+ private double edgeSample(long xx) {
+ int j = expSampleA();
+ if (j != 0) {
+ // Overhang frequency = 0.0151056
+ return expOverhang(j, xx);
+ }
+ // Tail frequency = 0.000515503
+
+ // Perform a new sample and add it to the start of the tail.
+ double x0 = X_0;
+ for (;;) {
+ final long x = nextLong();
+ final int i = ((int) x) & 0xff;
+
+ if (i < I_MAX) {
+ // Early exit.
+ return x0 + X[i] * (x >>> 1);
+ }
+ // Edge of the ziggurat
+ j = expSampleA();
+ if (j != 0) {
+ return x0 + expOverhang(j, x);
+ }
+ // Another tail sample
+ x0 += X_0;
+ }
+ }
+
+ /**
+ * Draws a PRN from overhang.
+ *
+ * <p>This does not use recursion.
+ *
+ * @param j Index j (must be {@code > 0})
+ * @param xx Initial random deviate
+ * @return the sample
+ */
+ protected double expOverhang(int j, long xx) {
+ // Recycle the initial random deviate.
+ // Shift right to make an unsigned long.
+ for (long ux = xx >>> 1;; ux = randomInt63()) {
+ // To sample a unit right-triangle:
+ // U_x <- min(U_1, U_2)
+ // distance <- | U_1 - U_2 |
+ // U_y <- 1 - (U_x + distance)
+ long uDistance = randomInt63() - ux;
+ if (uDistance < 0) {
+ uDistance = -uDistance;
+ ux -= uDistance;
+ }
+ // _FAST_PRNG_SAMPLE_X(xj, ux)
+ final double x = fastPrngSampleX(j, ux);
+ if (uDistance >= IE_MAX) {
+ // Frequency (per call into createSample): 0.0126230
+ // Frequency (per call into expOverhang): 0.823328
+ // Early Exit: x < y - epsilon
+ return x;
+ }
+ // Frequency per call into createSample:
+ // Return = 0.00248262
+ // Recursion = 0.000226013
+ // Frequency per call into expOverhang:
+ // Return = 0.161930
+ // Recursion = 0.0147417
+
+ // _FAST_PRNG_SAMPLE_Y(j, pow(2, 63) - (ux + uDistance))
+ // Long.MIN_VALUE is used as an unsigned int with value 2^63:
+ // uy = Long.MIN_VALUE - (ux + uDistance)
+ if (fastPrngSampleY(j, Long.MIN_VALUE - (ux + uDistance)) <= Math.exp(-x)) {
+ return x;
+ }
+ }
+ }
+ }
+
+ /**
+ * Modified Ziggurat method for sampling from an exponential distribution.
+ *
+ * <p>Uses the algorithm from McFarland, C.D. (2016).
+ *
* <p>This implementation separates sampling of the main ziggurat and the recursive
* sampling from the edge into different methods.
*/