1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 package org.apache.commons.rng.examples.jmh.sampling.distribution;
19
20 import org.apache.commons.rng.UniformRandomProvider;
21 import org.apache.commons.rng.sampling.distribution.ContinuousSampler;
22 import org.apache.commons.rng.sampling.distribution.StableSampler;
23 import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
24 import org.apache.commons.rng.simple.RandomSource;
25 import org.openjdk.jmh.annotations.Benchmark;
26 import org.openjdk.jmh.annotations.BenchmarkMode;
27 import org.openjdk.jmh.annotations.Fork;
28 import org.openjdk.jmh.annotations.Measurement;
29 import org.openjdk.jmh.annotations.Mode;
30 import org.openjdk.jmh.annotations.OutputTimeUnit;
31 import org.openjdk.jmh.annotations.Param;
32 import org.openjdk.jmh.annotations.Scope;
33 import org.openjdk.jmh.annotations.Setup;
34 import org.openjdk.jmh.annotations.State;
35 import org.openjdk.jmh.annotations.Warmup;
36
37 import java.util.concurrent.TimeUnit;
38
39 /**
40 * Executes a benchmark to compare the speed of generation of stable random numbers
41 * using different methods.
42 *
43 * <p>The test specifically includes a test of different versions of the sampler
44 * when {@code beta=0}; {@code alpha=1}; or the generic case {@code alpha != 1} and
45 * {@code beta != 0}.
46 */
47 @BenchmarkMode(Mode.AverageTime)
48 @OutputTimeUnit(TimeUnit.NANOSECONDS)
49 @Warmup(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
50 @Measurement(iterations = 5, time = 1, timeUnit = TimeUnit.SECONDS)
51 @State(Scope.Benchmark)
52 @Fork(value = 1, jvmArgs = {"-server", "-Xms128M", "-Xmx128M"})
53 public class StableSamplerPerformance {
54 /** The name of the baseline method. */
55 static final String BASELINE = "baseline";
56
57 /**
58 * Samples from a stable distribution with zero location and unit scale.
59 * Two implementations are available:
60 *
61 * <ol>
62 * <li>Implements the Chambers-Mallows-Stuck (CMS) method from Chambers, et al
63 * (1976) A Method for Simulating Stable Random Variables. Journal of the
64 * American Statistical Association Vol. 71, No. 354, pp. 340-344.
65 * This is a translation of the FORTRAN routine RSTAB; it computes a rearrangement
66 * of the original formula to use double-angle identities and output a 0-parameterization
67 * stable deviate.
68 *
69 * <li>The implementation uses the Chambers-Mallows-Stuck (CMS) method using
70 * the formula proven in Weron (1996) "On the Chambers-Mallows-Stuck method for
71 * simulating skewed stable random variables". Statistics & Probability
72 * Letters. 28 (2): 165–171. This outputs a 1-parameterization stable deviate and
73 * has been modified to the 0-parameterization using a translation.
74 * </ol>
75 *
76 * <p>The code has been copied from o.a.c.rng.sampling.distributions.StableSampler.
77 * The classes have been renamed to StableRandomGenerator to distinguish them.
78 * The copy implementation of the CMS algorithm RSTAB uses Math.tan to compute
79 * tan(x) / x. The main StableSampler implements the CMS RSTAB algorithm using
80 * fast approximation to the Math.tan function.
81 * The implementation of the Weron formula has been copied as it is not possible
82 * to instantiate these samplers directly. They are used for edge case computations.
83 *
84 * @see <a href="https://en.wikipedia.org/wiki/Stable_distribution">Stable
85 * distribution (Wikipedia)</a>
86 * @see <a href="https://doi.org/10.1080%2F01621459.1976.10480344">Chambers et
87 * al (1976) JOASA 71: 340-344</a>
88 * @see <a
89 * href="https://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.46.3280">Weron
90 * (1996) Statistics & Probability Letters. 28: 165–171</a>
91 */
92 abstract static class StableRandomGenerator implements ContinuousSampler {
93 /** The lower support for the distribution. This is the lower bound of {@code (-inf, +inf)}. */
94 static final double LOWER = Double.NEGATIVE_INFINITY;
95 /** The upper support for the distribution. This is the upper bound of {@code (-inf, +inf)}. */
96 static final double UPPER = Double.POSITIVE_INFINITY;
97 /** pi / 2. */
98 static final double PI_2 = Math.PI / 2;
99 /** pi / 4. */
100 static final double PI_4 = Math.PI / 4;
101 /** pi/2 scaled by 2^-53. */
102 static final double PI_2_SCALED = 0x1.0p-54 * Math.PI;
103 /** pi/4 scaled by 2^-53. */
104 static final double PI_4_SCALED = 0x1.0p-55 * Math.PI;
105 /** Underlying source of randomness. */
106 private final UniformRandomProvider rng;
107 /** The exponential sampler. */
108 private final ContinuousSampler expSampler;
109
110
111 /**
112 * @param rng Underlying source of randomness
113 */
114 StableRandomGenerator(final UniformRandomProvider rng) {
115 this.rng = rng;
116 expSampler = ZigguratSampler.Exponential.of(rng);
117 }
118
119 /**
120 * Gets a random value for the omega parameter.
121 * This is an exponential random variable with mean 1.
122 *
123 * @return omega
124 */
125 double getOmega() {
126 return expSampler.sample();
127 }
128
129 /**
130 * Gets a random value for the phi parameter.
131 * This is a uniform random variable in {@code (-pi/2, pi/2)}.
132 *
133 * @return phi
134 */
135 double getPhi() {
136 final double x = (rng.nextLong() >> 10) * PI_2_SCALED;
137 if (x == -PI_2) {
138 return getPhi();
139 }
140 return x;
141 }
142
143 /**
144 * Gets a random value for the phi/2 parameter.
145 * This is a uniform random variable in {@code (-pi/4, pi/4)}.
146 *
147 * @return phi/2
148 */
149 double getPhiBy2() {
150 final double x = (rng.nextLong() >> 10) * PI_4_SCALED;
151 if (x == -PI_4) {
152 return getPhiBy2();
153 }
154 return x;
155 }
156
157 /**
158 * Creates a sampler of a stable distribution with zero location and unit scale.
159 *
160 * <p>WARNING: Parameters are not validated.
161 *
162 * @param rng Generator of uniformly distributed random numbers.
163 * @param alpha Stability parameter. Must be in range {@code (0, 2]}.
164 * @param beta Skewness parameter. Must be in range {@code [-1, 1]}.
165 * @param weron Set to true to use the Weron formulas. Default is CMS.
166 * @return the sampler
167 */
168 static ContinuousSampler of(UniformRandomProvider rng,
169 double alpha,
170 double beta,
171 boolean weron) {
172 // WARNING:
173 // No parameter validation
174 // No Alpha=2 case
175 // No Alpha=1, beta=0 case
176 if (weron) {
177 if (beta == 0.0) {
178 return new Beta0WeronStableRandomGenerator(rng, alpha);
179 }
180 if (alpha == 1.0) {
181 return new Alpha1WeronStableRandomGenerator(rng, alpha);
182 }
183 return new WeronStableRandomGenerator(rng, alpha, beta);
184 }
185
186 if (beta == 0.0) {
187 return new Beta0CMSStableRandomGenerator(rng, alpha);
188 }
189 if (alpha == 1.0) {
190 return new Alpha1CMSStableRandomGenerator(rng, alpha);
191 }
192 return new CMSStableRandomGenerator(rng, alpha, beta);
193 }
194
195 /**
196 * Clip the sample {@code x} to the inclusive limits of the support.
197 *
198 * @param x Sample x.
199 * @param lower Lower bound.
200 * @param upper Upper bound.
201 * @return x in [lower, upper]
202 */
203 static double clipToSupport(double x, double lower, double upper) {
204 if (x < lower) {
205 return lower;
206 }
207 return x < upper ? x : upper;
208 }
209 }
210
211 /**
212 * A baseline for a generator that uses an exponential sample and a uniform deviate in the
213 * range {@code [-pi/2, pi/2]}.
214 */
215 private static class BaselineStableRandomGenerator extends StableRandomGenerator {
216 /**
217 * @param rng Underlying source of randomness
218 */
219 BaselineStableRandomGenerator(UniformRandomProvider rng) {
220 super(rng);
221 }
222
223 @Override
224 public double sample() {
225 return getOmega() * getPhi();
226 }
227 }
228
229 /**
230 * Implement the generic stable distribution case: {@code alpha < 2} and {@code beta != 0},
231 * and {@code alpha != 0}.
232 *
233 * <p>Implements the Weron formula.
234 */
235 private static class WeronStableRandomGenerator extends StableRandomGenerator {
236 /** The lower support for the distribution. */
237 protected final double lower;
238 /** The upper support for the distribution. */
239 protected final double upper;
240 /** Epsilon (1 - alpha). */
241 protected final double eps;
242 /** Epsilon (1 - alpha). */
243 protected final double meps1;
244 /** Cache of expression value used in generation. */
245 protected final double zeta;
246 /** Cache of expression value used in generation. */
247 protected final double atanZeta;
248 /** Cache of expression value used in generation. */
249 protected final double scale;
250 /** 1 / alpha = 1 / (1 - eps). */
251 protected final double inv1mEps;
252 /** (1 / alpha) - 1 = eps / (1 - eps). */
253 protected final double epsDiv1mEps;
254
255 /**
256 * @param rng Underlying source of randomness
257 * @param alpha Stability parameter. Must be in range {@code (0, 2]}.
258 * @param beta Skewness parameter. Must be in range {@code [-1, 1]}.
259 */
260 WeronStableRandomGenerator(UniformRandomProvider rng, double alpha, double beta) {
261 super(rng);
262
263 eps = 1 - alpha;
264 // When alpha < 0.5, 1 - eps == alpha is not always true as the reverse is not exact.
265 // Here we store 1 - eps in place of alpha. Thus eps + (1 - eps) = 1.
266 meps1 = 1 - eps;
267
268 // Compute pre-factors for the Weron formula used during error correction.
269 if (meps1 > 1) {
270 // Avoid calling tan outside the domain limit [-pi/2, pi/2].
271 zeta = beta * Math.tan((2 - meps1) * PI_2);
272 } else {
273 zeta = -beta * Math.tan(meps1 * PI_2);
274 }
275
276 // Do not store xi = Math.atan(-zeta) / meps1 due to floating-point division errors.
277 // Directly store Math.atan(-zeta).
278 atanZeta = Math.atan(-zeta);
279 scale = Math.pow(1 + zeta * zeta, 0.5 / meps1);
280 inv1mEps = 1.0 / meps1;
281 epsDiv1mEps = inv1mEps - 1;
282
283 // Compute the support. This applies when alpha < 1 and beta = +/-1
284 if (alpha < 1 && Math.abs(beta) == 1) {
285 if (beta == 1) {
286 // alpha < 0, beta = 1
287 lower = zeta;
288 upper = UPPER;
289 } else {
290 // alpha < 0, beta = -1
291 lower = LOWER;
292 upper = zeta;
293 }
294 } else {
295 lower = LOWER;
296 upper = UPPER;
297 }
298 }
299
300 @Override
301 public double sample() {
302 final double w = getOmega();
303 final double phi = getPhi();
304
305 // Add back zeta
306 // This creates the parameterization defined in Nolan (1998):
307 // X ~ S0_alpha(s,beta,u0) with s=1, u0=0 for a standard random variable.
308
309 double t1 = Math.sin(meps1 * phi + atanZeta) / Math.pow(Math.cos(phi), inv1mEps);
310 double t2 = Math.pow(Math.cos(eps * phi - atanZeta) / w, epsDiv1mEps);
311
312 // Term t1 and t2 can be zero or infinite.
313 // Used a boxed infinity to avoid inf * 0.
314 double unbox1 = 1;
315 double unbox2 = 1;
316 if (Double.isInfinite(t1)) {
317 t1 = Math.copySign(Double.MAX_VALUE, t1);
318 unbox1 = unbox2 = Double.MAX_VALUE;
319 }
320 if (Double.isInfinite(t2)) {
321 t2 = Math.copySign(Double.MAX_VALUE, t2);
322 unbox1 = unbox2 = Double.MAX_VALUE;
323 }
324 // Note: The order of the product must be maintained to unbox the infinity
325 return clipToSupport(t1 * t2 * unbox1 * unbox2 * scale + zeta, lower, upper);
326 }
327 }
328
329 /**
330 * Implement the {@code alpha == 1} and {@code beta != 0} stable distribution
331 * case.
332 *
333 * <p>Implements the Weron formula.
334 */
335 private static class Alpha1WeronStableRandomGenerator extends StableRandomGenerator {
336 /** Skewness parameter. */
337 private final double beta;
338
339 /**
340 * @param rng Underlying source of randomness
341 * @param beta Skewness parameter. Must be in range {@code [-1, 1]}.
342 */
343 Alpha1WeronStableRandomGenerator(UniformRandomProvider rng, double beta) {
344 super(rng);
345 this.beta = beta;
346 }
347
348 @Override
349 public double sample() {
350 final double w = getOmega();
351 final double phi = getPhi();
352 // Generic stable distribution with alpha = 1.
353 // Note: betaPhi cannot be zero when phi is limited to < pi/2.
354 // This eliminates divide by 0 errors.
355 final double betaPhi = PI_2 + beta * phi;
356 final double x = (betaPhi * Math.tan(phi) -
357 beta * Math.log(PI_2 * w * Math.cos(phi) / betaPhi)) / PI_2;
358 // When w -> 0 this computes +/- infinity.
359 return x;
360 }
361 }
362
363 /**
364 * Implement the {@code alpha < 2} and {@code beta = 0} stable distribution case.
365 *
366 * <p>Implements the Weron formula.
367 */
368 private static class Beta0WeronStableRandomGenerator extends StableRandomGenerator {
369 /** Epsilon (1 - alpha). */
370 protected final double eps;
371 /** Epsilon (1 - alpha). */
372 protected final double meps1;
373 /** 1 / alpha = 1 / (1 - eps). */
374 protected final double inv1mEps;
375 /** (1 / alpha) - 1 = eps / (1 - eps). */
376 protected final double epsDiv1mEps;
377
378 /**
379 * @param rng Underlying source of randomness
380 * @param alpha Stability parameter. Must be in range {@code (0, 2]}.
381 */
382 Beta0WeronStableRandomGenerator(UniformRandomProvider rng, double alpha) {
383 super(rng);
384 eps = 1 - alpha;
385 meps1 = 1 - eps;
386 inv1mEps = 1.0 / meps1;
387 epsDiv1mEps = inv1mEps - 1;
388 }
389
390 @Override
391 public double sample() {
392 final double w = getOmega();
393 final double phi = getPhi();
394 // Compute terms
395 double t1 = Math.sin(meps1 * phi) / Math.pow(Math.cos(phi), inv1mEps);
396 double t2 = Math.pow(Math.cos(eps * phi) / w, epsDiv1mEps);
397
398 double unbox1 = 1;
399 double unbox2 = 1;
400 if (Double.isInfinite(t1)) {
401 t1 = Math.copySign(Double.MAX_VALUE, t1);
402 unbox1 = unbox2 = Double.MAX_VALUE;
403 }
404 if (Double.isInfinite(t2)) {
405 t2 = Math.copySign(Double.MAX_VALUE, t2);
406 unbox1 = unbox2 = Double.MAX_VALUE;
407 }
408 // Note: The order of the product must be maintained to unbox the infinity
409 return t1 * t2 * unbox1 * unbox2;
410 }
411 }
412
413 /**
414 * Implement the generic stable distribution case: {@code alpha < 2} and
415 * {@code beta != 0}. This routine assumes {@code alpha != 1}.
416 *
417 * <p>Implements the CMS formula using the RSTAB algorithm.
418 */
419 static class CMSStableRandomGenerator extends WeronStableRandomGenerator {
420 /** 1/2. */
421 private static final double HALF = 0.5;
422 /** Cache of expression value used in generation. */
423 private final double tau;
424
425 /**
426 * @param rng Underlying source of randomness
427 * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
428 * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
429 */
430 CMSStableRandomGenerator(UniformRandomProvider rng, double alpha, double beta) {
431 super(rng, alpha, beta);
432
433 // Compute the RSTAB pre-factors.
434 // tau = -eps * tan(alpha * Phi0)
435 // with Phi0 = alpha^-1 * arctan(beta * tan(pi alpha / 2)).
436
437 // tau is symmetric around alpha=1
438 // tau -> beta / pi/2 as alpha -> 1
439 // tau -> 0 as alpha -> 2 or 0
440 // Avoid calling tan as the value approaches the domain limit [-pi/2, pi/2].
441 if (eps == 0) {
442 // alpha == 1
443 tau = beta / PI_2;
444 } else if (Math.abs(eps) < HALF) {
445 // 0.5 < alpha < 1.5
446 tau = eps * beta / Math.tan(eps * PI_2);
447 } else {
448 // alpha >= 1.5 or alpha <= 0.5.
449 // Do not call tan with alpha > 1 as it wraps in the domain [-pi/2, pi/2].
450 // Since pi is approximate the symmetry is lost by wrapping.
451 // Keep within the domain using (2-alpha).
452 if (meps1 > 1) {
453 tau = eps * beta * -Math.tan((2 - meps1) * PI_2);
454 } else {
455 tau = eps * beta * Math.tan(meps1 * PI_2);
456 }
457 }
458
459 }
460
461 @Override
462 public double sample() {
463 final double phiby2 = getPhiBy2();
464 final double w = getOmega();
465
466 // Compute as per the RSTAB routine.
467
468 // Generic stable distribution that is continuous as alpha -> 1.
469 // This is a trigonomic rearrangement of equation 4.1 from Chambers et al (1976)
470 // as implemented in the Fortran program RSTAB.
471 // Uses the special functions:
472 // tan2 = tan(x) / x
473 // d2 = (exp(x) - 1) / x
474 // The method is implemented as per the RSTAB routine but using Math.tan
475 // for the tan2 function
476
477 // Compute some tangents
478 final double a = Math.tan(phiby2);
479 final double b = Math.tan(eps * phiby2);
480 final double bb = b == 0 ? 1 : b / (eps * phiby2);
481
482 // Compute some necessary subexpressions
483 final double da = a * a;
484 final double db = b * b;
485 final double a2 = 1 - da;
486 final double a2p = 1 + da;
487 final double b2 = 1 - db;
488 final double b2p = 1 + db;
489
490 // Compute coefficient.
491 final double numerator = b2 + 2 * phiby2 * bb * tau;
492 final double z = a2p * numerator / (w * a2 * b2p);
493
494 // Compute the exponential-type expression
495 final double alogz = Math.log(z);
496 final double d = D2Source.d2(epsDiv1mEps * alogz) * (alogz * inv1mEps);
497
498 // Pre-compute the multiplication factor.
499 final double f = (2 * ((a - b) * (1 + a * b) - phiby2 * tau * bb * (b * a2 - 2 * a))) /
500 (a2 * b2p);
501
502 // Compute the stable deviate:
503 final double x = (1 + eps * d) * f + tau * d;
504
505 // Test the support
506 if (lower < x && x < upper) {
507 return x;
508 }
509
510 // No error correction path here!
511 // Return something so that the test for the support is not optimised away.
512 return 0;
513 }
514 }
515
516 /**
517 * Implement the stable distribution case: {@code alpha == 1} and {@code beta != 0}.
518 *
519 * <p>Implements the same algorithm as the {@link CMSStableRandomGenerator} with
520 * the {@code alpha} assumed to be 1.
521 *
522 * <p>This routine assumes {@code beta != 0}; {@code alpha=1, beta=0} is the Cauchy
523 * distribution case.
524 */
525 static class Alpha1CMSStableRandomGenerator extends StableRandomGenerator {
526 /** Cache of expression value used in generation. */
527 private final double tau;
528
529 /**
530 * @param rng Underlying source of randomness
531 * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
532 */
533 Alpha1CMSStableRandomGenerator(UniformRandomProvider rng, double beta) {
534 super(rng);
535 tau = beta / PI_2;
536 }
537
538 @Override
539 public double sample() {
540 final double phiby2 = getPhiBy2();
541 final double w = getOmega();
542
543 // Compute some tangents
544 final double a = Math.tan(phiby2);
545
546 // Compute some necessary subexpressions
547 final double da = a * a;
548 final double a2 = 1 - da;
549 final double a2p = 1 + da;
550
551 // Compute coefficient.
552 final double z = a2p * (1 + 2 * phiby2 * tau) / (w * a2);
553
554 // Compute the exponential-type expression
555 final double d = Math.log(z);
556
557 // Pre-compute the multiplication factor.
558 final double f = (2 * (a - phiby2 * tau * (-2 * a))) / a2;
559
560 // Compute the stable deviate:
561 return f + tau * d;
562 }
563 }
564
565 /**
566 * Implement the generic stable distribution case: {@code alpha < 2} and {@code beta == 0}.
567 *
568 * <p>Implements the same algorithm as the {@link CMSStableRandomGenerator} with
569 * the {@code beta} assumed to be 0.
570 *
571 * <p>This routine assumes {@code alpha != 1}; {@code alpha=1, beta=0} is the Cauchy
572 * distribution case.
573 */
574 static class Beta0CMSStableRandomGenerator extends Beta0WeronStableRandomGenerator {
575 /**
576 * @param rng Underlying source of randomness
577 * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
578 */
579 Beta0CMSStableRandomGenerator(UniformRandomProvider rng, double alpha) {
580 super(rng, alpha);
581 }
582
583 @Override
584 public double sample() {
585 final double phiby2 = getPhiBy2();
586 final double w = getOmega();
587
588 // Compute some tangents
589 final double a = Math.tan(phiby2);
590 final double b = Math.tan(eps * phiby2);
591 // Compute some necessary subexpressions
592 final double da = a * a;
593 final double db = b * b;
594 final double a2 = 1 - da;
595 final double a2p = 1 + da;
596 final double b2 = 1 - db;
597 final double b2p = 1 + db;
598 // Compute coefficient.
599 final double z = a2p * b2 / (w * a2 * b2p);
600
601 // Compute the exponential-type expression
602 final double alogz = Math.log(z);
603 final double d = D2Source.d2(epsDiv1mEps * alogz) * (alogz * inv1mEps);
604
605 // Pre-compute the multiplication factor.
606 final double f = (2 * ((a - b) * (1 + a * b))) / (a2 * b2p);
607
608 // Compute the stable deviate:
609 final double x = (1 + eps * d) * f;
610
611 // Test the support
612 if (LOWER < x && x < UPPER) {
613 return x;
614 }
615
616 // No error correction path here!
617 // Return something so that the test for the support is not optimised away.
618 return 0;
619 }
620 }
621
622 /**
623 * Defines the {@link RandomSource} for testing a {@link ContinuousSampler}.
624 */
625 public abstract static class SamplerSource {
626 /**
627 * RNG providers.
628 *
629 * <p>Use different speeds.</p>
630 *
631 * @see <a href="https://commons.apache.org/proper/commons-rng/userguide/rng.html">
632 * Commons RNG user guide</a>
633 */
634 @Param({"XO_RO_SHI_RO_128_PP",
635 "MWC_256",
636 "JDK"})
637 private String randomSourceName;
638
639 /**
640 * Gets the source of randomness.
641 *
642 * @return RNG
643 */
644 public UniformRandomProvider getRNG() {
645 return RandomSource.valueOf(randomSourceName).create();
646 }
647
648 /**
649 * @return the sampler.
650 */
651 public abstract ContinuousSampler getSampler();
652 }
653
654 /**
655 * Source for a uniform random deviate in an open interval, e.g. {@code (-0.5, 0.5)}.
656 */
657 @State(Scope.Benchmark)
658 public static class UniformRandomSource extends SamplerSource {
659 /** The lower limit of (-pi/4, pi/4). */
660 private static final double NEG_PI_4 = -Math.PI / 4;
661 /** pi/4 scaled by 2^-53. */
662 private static final double PI_4_SCALED = 0x1.0p-55 * Math.PI;
663
664 /** Method to generate the uniform deviate. */
665 @Param({BASELINE, "nextDoubleNot0", "nextDoubleNot0Recurse",
666 "nextLongNot0", "nextDoubleShifted", "nextLongShifted",
667 "signedShift", "signedShiftPi4", "signedShiftPi4b"})
668 private String method;
669
670 /** The sampler. */
671 private ContinuousSampler sampler;
672
673 /** {@inheritDoc} */
674 @Override
675 public ContinuousSampler getSampler() {
676 return sampler;
677 }
678
679 /** Instantiates sampler. */
680 @Setup
681 public void setup() {
682 final UniformRandomProvider rng = getRNG();
683 if (BASELINE.equals(method)) {
684 sampler = rng::nextDouble;
685 } else if ("nextDoubleNot0".equals(method)) {
686 sampler = () -> {
687 // Sample the 2^53 dyadic rationals in [0, 1) with zero excluded
688 double x;
689 do {
690 x = rng.nextDouble();
691 } while (x == 0);
692 return x - 0.5;
693 };
694 } else if ("nextDoubleNot0Recurse".equals(method)) {
695 sampler = new ContinuousSampler() {
696 @Override
697 public double sample() {
698 // Sample the 2^53 dyadic rationals in [0, 1) with zero excluded.
699 // Use recursion to generate a stack overflow with a bad provider.
700 // This is better than an infinite loop.
701 final double x = rng.nextDouble();
702 if (x == 0) {
703 return sample();
704 }
705 return x - 0.5;
706 }
707 };
708 } else if ("nextLongNot0".equals(method)) {
709 sampler = () -> {
710 // Sample the 2^53 dyadic rationals in [0, 1) with zero excluded
711 long x;
712 do {
713 x = rng.nextLong() >>> 11;
714 } while (x == 0);
715 return x * 0x1.0p-53 - 0.5;
716 };
717 } else if ("nextDoubleShifted".equals(method)) {
718 sampler = () -> {
719 // Sample the 2^52 dyadic rationals in [0, 1) and shift by 2^-53.
720 // No infinite loop but the deviate loses 1 bit of randomness.
721 return 0x1.0p-53 + (rng.nextLong() >>> 12) * 0x1.0p-52 - 0.5;
722 };
723 } else if ("nextLongShifted".equals(method)) {
724 sampler = () -> {
725 // Sample the 2^53 dyadic rationals in [0, 1) but set the lowest
726 // bit. This result in 2^52 dyadic rationals in (0, 1) to avoid 0.
727 return ((rng.nextLong() >>> 11) | 0x1L) * 0x1.0p-53 - 0.5;
728 };
729 } else if ("signedShift".equals(method)) {
730 sampler = new ContinuousSampler() {
731 @Override
732 public double sample() {
733 // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but using a
734 // signed shift of 10 in place of an unsigned shift of 11, and updating the
735 // multiplication factor to 2^-54.
736 final double x = (rng.nextLong() >> 10) * 0x1.0p-54;
737 if (x == -0.5) {
738 // Avoid the extreme of the bounds
739 return sample();
740 }
741 return x;
742 }
743 };
744 } else if ("signedShiftPi4".equals(method)) {
745 // Note: This does generate u in (-0.5, 0.5) which must be scaled by pi
746 // or pi/2 for the CMS algorithm but directly generates in (-pi/4, pi/4).
747 sampler = new ContinuousSampler() {
748 @Override
749 public double sample() {
750 // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but using a
751 // signed shift of 10 in place of an unsigned shift of 11, and updating the
752 // multiplication factor to 2^-54 * pi/2
753 final double x = (rng.nextLong() >> 10) * PI_4_SCALED;
754 if (x == NEG_PI_4) {
755 // Avoid the extreme of the bounds
756 return sample();
757 }
758 return x;
759 }
760 };
761 } else if ("signedShiftPi4b".equals(method)) {
762 // Note: This does generate u in (-0.5, 0.5) which must be scaled by pi
763 // or pi/2 for the CMS algorithm but directly generates in (-pi/4, pi/4).
764 sampler = new ContinuousSampler() {
765 @Override
766 public double sample() {
767 final long x = rng.nextLong();
768 if (x == Long.MIN_VALUE) {
769 // Avoid the extreme of the bounds
770 return sample();
771 }
772 // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but using a
773 // signed shift of 10 in place of an unsigned shift of 11, and updating the
774 // multiplication factor to 2^-54 * pi/2
775 return (x >> 10) * PI_4_SCALED;
776 }
777 };
778 } else {
779 throw new IllegalStateException("Unknown method: " + method);
780 }
781 }
782 }
783
784 /**
785 * Source for testing implementations of tan(x) / x. The function must work on a value
786 * in the range {@code [0, pi/4]}.
787 *
788 * <p>The tan(x) / x function is required for the trigonomic rearrangement of the CMS
789 * formula to create a stable random variate.
790 */
791 @State(Scope.Benchmark)
792 public static class TanSource extends SamplerSource {
793 /** pi / 2. */
794 private static final double PI_2 = Math.PI / 2;
795 /** pi / 4. */
796 private static final double PI_4 = Math.PI / 4;
797 /** 4 / pi. */
798 private static final double FOUR_PI = 4 / Math.PI;
799
800 /** Method to generate tan(x) / x. */
801 @Param({BASELINE, "tan", "tan4283", "tan4288", "tan4288b", "tan4288c"})
802 private String method;
803
804 /** The sampler. */
805 private ContinuousSampler sampler;
806
807 /** {@inheritDoc} */
808 @Override
809 public ContinuousSampler getSampler() {
810 return sampler;
811 }
812
813 /** Instantiates sampler. */
814 @Setup
815 public void setup() {
816 final UniformRandomProvider rng = getRNG();
817 if (BASELINE.equals(method)) {
818 sampler = () -> {
819 // A value in [-pi/4, pi/4]
820 return PI_2 * (rng.nextDouble() - 0.5);
821 };
822 } else if ("tan".equals(method)) {
823 sampler = () -> {
824 final double x = PI_2 * (rng.nextDouble() - 0.5);
825 // Require tan(0) / 0 = 1 and not NaN
826 return x == 0 ? 1.0 : Math.tan(x) / x;
827 };
828 } else if ("tan4283".equals(method)) {
829 sampler = () -> {
830 final double x = PI_2 * (rng.nextDouble() - 0.5);
831 return x * tan4283(x);
832 };
833 } else if ("tan4288".equals(method)) {
834 sampler = () -> {
835 final double x = PI_2 * (rng.nextDouble() - 0.5);
836 return x * tan4288(x);
837 };
838 } else if ("tan4288b".equals(method)) {
839 sampler = () -> {
840 final double x = PI_2 * (rng.nextDouble() - 0.5);
841 return x * tan4288b(x);
842 };
843 } else if ("tan4288c".equals(method)) {
844 sampler = () -> {
845 final double x = PI_2 * (rng.nextDouble() - 0.5);
846 return x * tan4288c(x);
847 };
848 } else {
849 throw new IllegalStateException("Unknown tan method: " + method);
850 }
851 }
852
853 // Mean (max) ULP is very low
854 // 2^30 random samples in [0, pi / 4)
855 // tan(x) tan(x) / x
856 // tan4283 93436.25534446817 201185 : 68313.16171793547 128079
857 // tan4288 0.5898815048858523 4 : 0.40416045393794775 3
858 // tan4288b 0.8608690425753593 8 : 0.6117749745026231 5
859 // tan4288c 0.5905972588807344 4 : 0.4047176940366626 3
860
861 /**
862 * Evaluate {@code tan(x) / x}.
863 *
864 * <p>For {@code x} in the range {@code [0, pi/4]} this ideally should return
865 * a value in the range {@code [1, 4 / pi]}. However due to lack of precision
866 * this is not the case for the extremes of the interval. It is not recommended
867 * to use this function in the CMS algorithm when computing the routine with
868 * double precision. The original RSTAB routine is for single precision floats.
869 *
870 * <p>Uses tan 4283 from Hart, JF et al (1968) Computer Approximations,
871 * New York: John Wiley & Sons, Inc.
872 *
873 * <p>This is the original method used in the RSTAB function by Chambers et al (1976).
874 * It has been updated to use the full precision constants provided in Hart.
875 *
876 * @param x the x
877 * @return {@code tan(x) / x}.
878 */
879 static double tan4283(double x) {
880 double xa = Math.abs(x);
881 if (xa > PI_4) {
882 return Math.tan(x) / x;
883 }
884 // Approximation 4283 from Hart et al (1968, P. 251).
885 // Max relative error = 1e-10.66.
886 // When testing verses Math.tan(x) the mean ULP difference was 93436.3
887 final double p0 = 0.129221035031569917e3;
888 final double p1 = -0.8876623770211723e1;
889 final double p2 = 0.52864445522248e-1;
890 final double q0 = 0.164529331810168605e3;
891 final double q1 = -0.45132056100598961e2;
892 // q2 = 1.0
893
894 xa = xa / PI_4;
895 final double xx = xa * xa;
896
897 // Polynomial implemented as per Chambers (1976) using the nested form.
898 return (p0 + xx * (p1 + xx * p2)) / (PI_4 * (q0 + xx * (q1 + xx /* * q2 */)));
899 }
900
901 /**
902 * Evaluate {@code tan(x) / x}.
903 *
904 * <p>For {@code x} in the range {@code [0, pi/4]} this returns
905 * a value in the range {@code [1, 4 / pi]}.
906 *
907 * <p>Uses tan 4288 from Hart, JF et al (1968) Computer Approximations,
908 * New York: John Wiley & Sons, Inc.
909 *
910 * <p>This is a higher precision method provided in Hart.
911 * It has the properties that {@code tan(x) / x = 1, x = 0} and
912 * {@code tan(x) = 1, x = pi / 4}.
913 *
914 * @param x the x
915 * @return {@code tan(x) / x}.
916 */
917 static double tan4288(double x) {
918 double xa = Math.abs(x);
919 if (xa > PI_4) {
920 return Math.tan(x) / x;
921 }
922
923 // Approximation 4288 from Hart et al (1968, P. 252).
924 // Max relative error = 1e-26.68 (for tan(x))
925 // When testing verses Math.tan(x) the mean ULP difference was 0.589882
926 final double p0 = -0.5712939549476836914932149599e+10;
927 final double p1 = +0.4946855977542506692946040594e+9;
928 final double p2 = -0.9429037070546336747758930844e+7;
929 final double p3 = +0.5282725819868891894772108334e+5;
930 final double p4 = -0.6983913274721550913090621370e+2;
931
932 final double q0 = -0.7273940551075393257142652672e+10;
933 final double q1 = +0.2125497341858248436051062591e+10;
934 final double q2 = -0.8000791217568674135274814656e+8;
935 final double q3 = +0.8232855955751828560307269007e+6;
936 final double q4 = -0.2396576810261093558391373322e+4;
937 // q5 = 1.0
938
939 xa = xa * FOUR_PI;
940 final double xx = xa * xa;
941
942 // Polynomial implemented as per Chambers (1976) using the nested form.
943 return (p0 + xx * (p1 + xx * (p2 + xx * (p3 + xx * p4)))) /
944 (PI_4 * (q0 + xx * (q1 + xx * (q2 + xx * (q3 + xx * (q4 + xx /* * q5 */))))));
945 }
946
947 /**
948 * Evaluate {@code tan(x) / x}.
949 *
950 * <p>For {@code x} in the range {@code [0, pi/4]} this returns
951 * a value in the range {@code [1, 4 / pi]}.
952 *
953 * <p>Uses tan 4288 from Hart, JF et al (1968) Computer Approximations,
954 * New York: John Wiley & Sons, Inc.
955 *
956 * @param x the x
957 * @return {@code tan(x) / x}.
958 */
959 static double tan4288b(double x) {
960 double xa = Math.abs(x);
961 if (xa > PI_4) {
962 return Math.tan(x) / x;
963 }
964
965 final double p0 = -0.5712939549476836914932149599e+10;
966 final double p1 = +0.4946855977542506692946040594e+9;
967 final double p2 = -0.9429037070546336747758930844e+7;
968 final double p3 = +0.5282725819868891894772108334e+5;
969 final double p4 = -0.6983913274721550913090621370e+2;
970
971 final double q0 = -0.7273940551075393257142652672e+10;
972 final double q1 = +0.2125497341858248436051062591e+10;
973 final double q2 = -0.8000791217568674135274814656e+8;
974 final double q3 = +0.8232855955751828560307269007e+6;
975 final double q4 = -0.2396576810261093558391373322e+4;
976 // q5 = 1.0
977
978 xa = xa * FOUR_PI;
979 // Rearrange the polynomial to the power form.
980 // Allows parallel computation of terms.
981 // This is faster but has been tested to have lower accuracy verses Math.tan.
982 // When testing verses Math.tan(x) the mean ULP difference was 0.860869
983 final double x2 = xa * xa;
984 final double x4 = x2 * x2;
985 final double x6 = x4 * x2;
986 final double x8 = x4 * x4;
987
988 return (p0 + x2 * p1 + x4 * p2 + x6 * p3 + x8 * p4) /
989 (PI_4 * (q0 + x2 * q1 + x4 * q2 + x6 * q3 + x8 * q4 + x8 * x2));
990 }
991
992 /**
993 * Evaluate {@code tan(x) / x}.
994 *
995 * <p>For {@code x} in the range {@code [0, pi/4]} this returns
996 * a value in the range {@code [1, 4 / pi]}.
997 *
998 * <p>Uses tan 4288 from Hart, JF et al (1968) Computer Approximations,
999 * New York: John Wiley & Sons, Inc.
1000 *
1001 * @param x the x
1002 * @return {@code tan(x) / x}.
1003 */
1004 static double tan4288c(double x) {
1005 if (Math.abs(x) > PI_4) {
1006 return Math.tan(x) / x;
1007 }
1008
1009 final double p0 = -0.5712939549476836914932149599e+10;
1010 final double p1 = +0.4946855977542506692946040594e+9;
1011 final double p2 = -0.9429037070546336747758930844e+7;
1012 final double p3 = +0.5282725819868891894772108334e+5;
1013 final double p4 = -0.6983913274721550913090621370e+2;
1014
1015 final double q0 = -0.7273940551075393257142652672e+10;
1016 final double q1 = +0.2125497341858248436051062591e+10;
1017 final double q2 = -0.8000791217568674135274814656e+8;
1018 final double q3 = +0.8232855955751828560307269007e+6;
1019 final double q4 = -0.2396576810261093558391373322e+4;
1020 // q5 = 1.0
1021
1022 final double xi = x * FOUR_PI;
1023 // Rearrange the polynomial to the power form.
1024 // Allows parallel computation of terms.
1025 // This is faster and has been tested to have accuracy similar to the nested form.
1026 // When testing verses Math.tan(x) the mean ULP difference was 0.590597
1027 final double x2 = xi * xi;
1028 final double x4 = x2 * x2;
1029 final double x6 = x4 * x2;
1030 final double x8 = x4 * x4;
1031
1032 // Reverse summation order to have least significant terms first.
1033 return (x8 * p4 + x6 * p3 + x4 * p2 + x2 * p1 + p0) /
1034 (PI_4 * (x8 * x2 + x8 * q4 + x6 * q3 + x4 * q2 + x2 * q1 + q0));
1035 }
1036 }
1037
1038 /**
1039 * Source for testing implementations of {@code (exp(x) - 1) / x}. The function must
1040 * work on a value in the range {@code [-inf, +inf]}.
1041 *
1042 * <p>The d2 function is required for the trigonomic rearrangement of the CMS
1043 * formula to create a stable random variate.
1044 *
1045 * <p>For testing the value x is generated as a uniform deviate with a range around
1046 * zero ({@code [-scale/2, +scale/2]}) using a configurable scale parameter.
1047 */
1048 @State(Scope.Benchmark)
1049 public static class D2Source extends SamplerSource {
1050 /** Method to generate (exp(x) - 1) / x. */
1051 @Param({BASELINE, "expm1", "expm1b", "exp", "hybrid"})
1052 private String method;
1053 /** Scale for the random value x. */
1054 @Param({"0.12345", "1.2345", "12.345", "123.45"})
1055 private double scale;
1056
1057 /** The sampler. */
1058 private ContinuousSampler sampler;
1059
1060 /** {@inheritDoc} */
1061 @Override
1062 public ContinuousSampler getSampler() {
1063 return sampler;
1064 }
1065
1066 /** Instantiates sampler. */
1067 @Setup
1068 public void setup() {
1069 final double s = scale;
1070 final UniformRandomProvider rng = getRNG();
1071 if (BASELINE.equals(method)) {
1072 sampler = () -> s * (rng.nextDouble() - 0.5);
1073 } else if ("expm1".equals(method)) {
1074 sampler = () -> expm1(s * (rng.nextDouble() - 0.5));
1075 } else if ("expm1b".equals(method)) {
1076 sampler = () -> expm1b(s * (rng.nextDouble() - 0.5));
1077 } else if ("exp".equals(method)) {
1078 sampler = () -> exp(s * (rng.nextDouble() - 0.5));
1079 } else if ("hybrid".equals(method)) {
1080 sampler = () -> hybrid(s * (rng.nextDouble() - 0.5));
1081 } else {
1082 throw new IllegalStateException("Unknown d2 method: " + method);
1083 }
1084 }
1085
1086 /**
1087 * Evaluate {@code (exp(x) - 1) / x}. For {@code x} in the range {@code [-inf, inf]} returns
1088 * a result in {@code [0, inf]}.
1089 *
1090 * <ul>
1091 * <li>For {@code x=-inf} this returns {@code 0}.
1092 * <li>For {@code x=0} this returns {@code 1}.
1093 * <li>For {@code x=inf} this returns {@code inf}.
1094 * </ul>
1095 *
1096 * <p> This corrects {@code 0 / 0} and {@code inf / inf} division from
1097 * {@code NaN} to either {@code 1} or the upper bound respectively.
1098 *
1099 * @param x value to evaluate
1100 * @return {@code (exp(x) - 1) / x}.
1101 */
1102 static double d2(double x) {
1103 return hybrid(x);
1104 }
1105
1106 /**
1107 * Evaluate {@code (exp(x) - 1) / x}.
1108 *
1109 * @param x value to evaluate
1110 * @return {@code (exp(x) - 1) / x}.
1111 */
1112 static double expm1(double x) {
1113 // Here we use a conditional to detect both edge cases, which are then corrected.
1114 final double d2 = Math.expm1(x) / x;
1115 if (Double.isNaN(d2)) {
1116 // Correct edge cases.
1117 if (x == 0) {
1118 return 1.0;
1119 }
1120 // x must have been +infinite or NaN
1121 return x;
1122 }
1123 return d2;
1124 }
1125
1126 /**
1127 * Evaluate {@code (exp(x) - 1) / x}.
1128 *
1129 * @param x value to evaluate
1130 * @return {@code (exp(x) - 1) / x}.
1131 */
1132 static double expm1b(double x) {
1133 // Edge cases
1134 if (x == 0) {
1135 return 1;
1136 }
1137 if (x == Double.POSITIVE_INFINITY) {
1138 return x;
1139 }
1140 return Math.expm1(x) / x;
1141 }
1142
1143 /**
1144 * Evaluate {@code (exp(x) - 1) / x}.
1145 *
1146 * @param x value to evaluate
1147 * @return {@code (exp(x) - 1) / x}.
1148 */
1149 static double exp(double x) {
1150 // No use of expm1.
1151 // Here we use a conditional to detect both edge cases, which are then corrected.
1152 final double d2 = (Math.exp(x) - 1) / x;
1153 if (Double.isNaN(d2)) {
1154 // Correct edge cases.
1155 if (x == 0) {
1156 return 1.0;
1157 }
1158 // x must have been +infinite or NaN
1159 return x;
1160 }
1161 return d2;
1162 }
1163
1164 /**
1165 * Evaluate {@code (exp(x) - 1) / x}.
1166 *
1167 * @param x value to evaluate
1168 * @return {@code (exp(x) - 1) / x}.
1169 */
1170 static double hybrid(double x) {
1171 // This is the threshold where use of expm1 and exp consistently
1172 // compute different results by more than 0.5 ULP.
1173 if (Math.abs(x) < 0.5) {
1174 if (x == 0) {
1175 return 1;
1176 }
1177 return Math.expm1(x) / x;
1178 }
1179 // No use of expm1. Accuracy as x moves away from 0 is not required as the result
1180 // is divided by x and the accuracy of the final result is the same.
1181 if (x == Double.POSITIVE_INFINITY) {
1182 return x;
1183 }
1184 return (Math.exp(x) - 1) / x;
1185 }
1186 }
1187
1188 /**
1189 * Baseline with an exponential deviate and a uniform deviate.
1190 */
1191 @State(Scope.Benchmark)
1192 public static class BaselineSamplerSource extends SamplerSource {
1193 /** The sampler. */
1194 private ContinuousSampler sampler;
1195
1196 /** {@inheritDoc} */
1197 @Override
1198 public ContinuousSampler getSampler() {
1199 return sampler;
1200 }
1201
1202 /** Instantiates sampler. */
1203 @Setup
1204 public void setup() {
1205 sampler = new BaselineStableRandomGenerator(getRNG());
1206 }
1207 }
1208
1209 /**
1210 * The sampler to use for testing. Defines the RandomSource and the type of
1211 * stable distribution sampler.
1212 */
1213 public abstract static class StableSamplerSource extends SamplerSource {
1214 /** The sampler type. */
1215 @Param({"CMS", "CMS_weron", "CMS_tan"})
1216 private String samplerType;
1217
1218 /** The sampler. */
1219 private ContinuousSampler sampler;
1220
1221 /**
1222 * The alpha value.
1223 *
1224 * @return alpha
1225 */
1226 abstract double getAlpha();
1227
1228 /**
1229 * The beta value.
1230 *
1231 * @return beta
1232 */
1233 abstract double getBeta();
1234
1235 /** {@inheritDoc} */
1236 @Override
1237 public ContinuousSampler getSampler() {
1238 return sampler;
1239 }
1240
1241 /** Instantiates sampler. */
1242 @Setup
1243 public void setup() {
1244 final UniformRandomProvider rng = getRNG();
1245 if ("CMS".equals(samplerType)) {
1246 sampler = StableSampler.of(rng, getAlpha(), getBeta());
1247 } else if ("CMS_weron".equals(samplerType)) {
1248 sampler = StableRandomGenerator.of(rng, getAlpha(), getBeta(), true);
1249 } else if ("CMS_tan".equals(samplerType)) {
1250 sampler = StableRandomGenerator.of(rng, getAlpha(), getBeta(), false);
1251 } else {
1252 throw new IllegalStateException("Unknown sampler: " + samplerType);
1253 }
1254 }
1255 }
1256
1257 /**
1258 * Sampling with {@code alpha = 1} and {@code beta != 0}.
1259 */
1260 @State(Scope.Benchmark)
1261 public static class Alpha1StableSamplerSource extends StableSamplerSource {
1262 /** The beta value. */
1263 @Param({"0.1", "0.4", "0.7"})
1264 private double beta;
1265
1266 /** {@inheritDoc} */
1267 @Override
1268 double getAlpha() {
1269 return 1;
1270 }
1271
1272 /** {@inheritDoc} */
1273 @Override
1274 double getBeta() {
1275 return beta;
1276 }
1277 }
1278
1279 /**
1280 * Sampling with {@code alpha != 1} and {@code beta = 0}.
1281 */
1282 @State(Scope.Benchmark)
1283 public static class Beta0StableSamplerSource extends StableSamplerSource {
1284 /** The alpha value. Use a range around 1. */
1285 @Param({
1286 //"0.2", "0.3", "0.4", "0.5", "0.6", "0.7", "0.8", "0.9", "0.99", "1.1", "1.2", "1.3", "1.4", "1.5"
1287 // The Weron formula is fast at alpha=0.5; possibly due to the two power
1288 // functions using a power of 1 and 2 respectively.
1289 //"0.45", "0.48", "0.5", "0.52", "0.55"
1290 "0.3", "0.5", "0.7", "0.9", "1.1", "1.3", "1.5", "1.7"
1291 })
1292 private double alpha;
1293
1294 /** {@inheritDoc} */
1295 @Override
1296 double getAlpha() {
1297 return alpha;
1298 }
1299
1300 /** {@inheritDoc} */
1301 @Override
1302 double getBeta() {
1303 return 0;
1304 }
1305 }
1306
1307 /**
1308 * Sampling with {@code alpha != 1} and {@code beta != 0}.
1309 */
1310 @State(Scope.Benchmark)
1311 public static class GeneralStableSamplerSource extends Beta0StableSamplerSource {
1312 /** The beta value. */
1313 @Param({"0.1", "0.4", "0.7"})
1314 private double beta;
1315
1316 /** {@inheritDoc} */
1317 @Override
1318 double getBeta() {
1319 return beta;
1320 }
1321 }
1322
1323 /**
1324 * Test methods for producing a uniform deviate in an open interval, e.g. {@code (-0.5, 0.5)}.
1325 * This is a requirement of the stable distribution CMS algorithm.
1326 *
1327 * @param source Source of randomness.
1328 * @return the {@code double} value
1329 */
1330 @Benchmark
1331 public double nextUniformDeviate(UniformRandomSource source) {
1332 return source.getSampler().sample();
1333 }
1334
1335 /**
1336 * Test methods for producing {@code tan(x) / x} in the range {@code [0, pi/4]}.
1337 * This is a requirement of the stable distribution CMS algorithm.
1338 *
1339 * @param source Source of randomness.
1340 * @return the {@code tan(x)} value
1341 */
1342 @Benchmark
1343 public double nextTan(TanSource source) {
1344 return source.getSampler().sample();
1345 }
1346
1347 /**
1348 * Test methods for producing {@code (exp(x) - 1) / x}.
1349 * This is a requirement of the stable distribution CMS algorithm.
1350 *
1351 * @param source Source of randomness.
1352 * @return the {@code (exp(x) - 1) / x} value
1353 */
1354 @Benchmark
1355 public double nextD2(D2Source source) {
1356 return source.getSampler().sample();
1357 }
1358
1359 /**
1360 * Baseline for any sampler that uses an exponential deviate and a uniform deviate.
1361 *
1362 * @param source Source of randomness.
1363 * @return the {@code double} value
1364 */
1365 @Benchmark
1366 public double sampleBaseline(BaselineSamplerSource source) {
1367 return source.getSampler().sample();
1368 }
1369
1370 /**
1371 * Run the stable sampler with {@code alpha = 1}.
1372 *
1373 * @param source Source of randomness.
1374 * @return the sample value
1375 */
1376 @Benchmark
1377 public double sampleAlpha1(Alpha1StableSamplerSource source) {
1378 return source.getSampler().sample();
1379 }
1380
1381 /**
1382 * Run the stable sampler with {@code beta = 0}.
1383 *
1384 * @param source Source of randomness.
1385 * @return the sample value
1386 */
1387 @Benchmark
1388 public double sampleBeta0(Beta0StableSamplerSource source) {
1389 return source.getSampler().sample();
1390 }
1391
1392 /**
1393 * Run the stable sampler.
1394 *
1395 * @param source Source of randomness.
1396 * @return the sample value
1397 */
1398 @Benchmark
1399 public double sample(GeneralStableSamplerSource source) {
1400 return source.getSampler().sample();
1401 }
1402 }