/*
    * Licensed to the Apache Software Foundation (ASF) under one or more
    * contributor license agreements.  See the NOTICE file distributed with
    * this work for additional information regarding copyright ownership.
    * The ASF licenses this file to You under the Apache License, Version 2.0
    * (the "License"); you may not use this file except in compliance with
    * the License.  You may obtain a copy of the License at
    *
    *      http://www.apache.org/licenses/LICENSE-2.0
   *
   * Unless required by applicable law or agreed to in writing, software
   * distributed under the License is distributed on an "AS IS" BASIS,
   * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   * See the License for the specific language governing permissions and
   * limitations under the License.
   */
  package org.apache.commons.rng.sampling.distribution;
  
  import java.util.function.Supplier;
  import org.apache.commons.rng.UniformRandomProvider;
  import org.apache.commons.rng.core.source64.SplitMix64;
  import org.apache.commons.rng.sampling.RandomAssert;
  import org.apache.commons.rng.sampling.distribution.StableSampler.Beta0CMSStableSampler;
  import org.apache.commons.rng.sampling.distribution.StableSampler.Beta0WeronStableSampler;
  import org.apache.commons.rng.sampling.distribution.StableSampler.CMSStableSampler;
  import org.apache.commons.rng.sampling.distribution.StableSampler.SpecialMath;
  import org.apache.commons.rng.sampling.distribution.StableSampler.WeronStableSampler;
  import org.apache.commons.rng.simple.RandomSource;
  import org.junit.jupiter.api.Assertions;
  import org.junit.jupiter.api.Test;
  
  /**
   * Tests for the class {@link StableSampler}.
   *
   * <p>Note: Samples from the stable distribution are tested in
   * {@link ContinuousSamplerParametricTest}.
   *
   * <p>This contains tests for the assumptions made by the {@link StableSampler} implementation
   * of the Chambers-Mallows-Stuck (CMS) method as described in
   * Chambers, Mallows &amp; Stuck (1976) "A Method for Simulating Stable Random Variables".
   * Journal of the American Statistical Association. 71 (354): 340–344.
   *
   * <p>The test class contains copy implementations of the routines in the {@link StableSampler}
   * to test the algorithms with various parameters. This avoids excess manipulation
   * of the RNG provided to the stable sampler to test edge cases and also allows
   * calling the algorithm with values that are eliminated by the sampler (e.g. u=0).
   *
   * <p>Some tests of the sampler are performed that manipulate the underlying RNG to create
   * extreme values for the random deviates. This hits edges cases where the computation has
   * to be corrected.
   */
  class StableSamplerTest {
      /** pi / 2. */
      private static final double PI_2 = Math.PI / 2;
      /** pi / 4. */
      private static final double PI_4 = Math.PI / 4;
      /** pi/4 scaled by 2^-53. */
      private static final double PI_4_SCALED = 0x1.0p-55 * Math.PI;
      /** The interval between successive values of a uniform variate u.
       * This is the gap between the 2^53 dyadic rationals in [0, 1). */
      private static final double DU = 0x1.0p-53;
      /** The smallest non-zero sample from the ZigguratSampler.Exponential sampler. */
      private static final double SMALL_W = 6.564735882096453E-19;
      /** A tail sample from the ZigguratSampler.Exponential after 1 recursions of the sample method.  */
      private static final double TAIL_W = 7.569274694148063;
      /** A largest sample from the ZigguratSampler.Exponential after 4 recursions of the sample method.  */
      private static final double LARGE_W = 4 * TAIL_W;
      /** The smallest value for alpha where 1 - (1-alpha) = alpha. */
      private static final double SMALLEST_ALPHA = 1.0 - Math.nextDown(1.0);
  
      private static final double VALID_ALPHA = 1.23;
      private static final double VALID_BETA = 0.23;
      private static final double VALID_GAMMA = 2.34;
      private static final double VALID_DELTA = 3.45;
  
      @Test
      void testAlphaZeroThrows() {
          assertConstructorThrows(0.0, VALID_BETA, VALID_GAMMA, VALID_DELTA);
      }
  
      @Test
      void testAlphaBelowZeroThrows() {
          assertConstructorThrows(Math.nextDown(0.0), VALID_BETA, VALID_GAMMA, VALID_DELTA);
      }
  
      @Test
      void testAlphaTooCloseToZeroThrows() {
          // The realistic range for alpha is not Double.MIN_VALUE.
          // The number 1 - alpha must not be 1.
          // This is valid
          final UniformRandomProvider rng = new SplitMix64(0L);
          StableSampler s = StableSampler.of(rng, SMALLEST_ALPHA, VALID_BETA, VALID_GAMMA, VALID_DELTA);
          Assertions.assertNotNull(s);
  
          // Smaller than this is still above zero but 1 - alpha == 1
          final double alphaTooSmall = SMALLEST_ALPHA / 2;
          Assertions.assertNotEquals(0.0, alphaTooSmall, "Expected alpha to be positive");
          Assertions.assertEquals(1.0, 1 - alphaTooSmall, "Expected rounding to 1");
  
         // Because alpha is effectively zero this will throw
         assertConstructorThrows(alphaTooSmall, VALID_BETA, VALID_GAMMA, VALID_DELTA);
     }
 
     @Test
     void testAlphaAboveTwoThrows() {
         assertConstructorThrows(Math.nextUp(2.0), VALID_BETA, VALID_GAMMA, VALID_DELTA);
     }
 
     @Test
     void testAlphaNaNThrows() {
         assertConstructorThrows(Double.NaN, VALID_BETA, VALID_GAMMA, VALID_DELTA);
     }
 
     @Test
     void testBetaBelowMinusOneThrows() {
         assertConstructorThrows(VALID_ALPHA, Math.nextDown(-1.0), VALID_GAMMA, VALID_DELTA);
     }
 
     @Test
     void testBetaAboveOneThrows() {
         assertConstructorThrows(VALID_ALPHA, Math.nextUp(1.0), VALID_GAMMA, VALID_DELTA);
     }
 
     @Test
     void testBetaNaNThrows() {
         assertConstructorThrows(VALID_ALPHA, Double.NaN, VALID_GAMMA, VALID_DELTA);
     }
 
     @Test
     void testGammaNotStrictlyPositiveThrows() {
         assertConstructorThrows(VALID_ALPHA, VALID_BETA, 0.0, VALID_DELTA);
     }
 
     @Test
     void testGammaInfThrows() {
         assertConstructorThrows(VALID_ALPHA, VALID_BETA, Double.POSITIVE_INFINITY, VALID_DELTA);
     }
 
     @Test
     void testGammaNaNThrows() {
         assertConstructorThrows(VALID_ALPHA, VALID_BETA, Double.NaN, VALID_DELTA);
     }
 
     @Test
     void testDeltaInfThrows() {
         assertConstructorThrows(VALID_ALPHA, VALID_BETA, VALID_GAMMA, Double.POSITIVE_INFINITY);
     }
 
     @Test
     void testDeltaNegInfThrows() {
         assertConstructorThrows(VALID_ALPHA, VALID_BETA, VALID_GAMMA, Double.NEGATIVE_INFINITY);
     }
 
     @Test
     void testDeltaNaNThrows() {
         assertConstructorThrows(VALID_ALPHA, VALID_BETA, VALID_GAMMA, Double.NaN);
     }
 
     /**
      * Asserts the stable sampler factory constructor throws an {@link IllegalArgumentException}.
      *
      * @param alpha Stability parameter. Must be in range {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in range {@code [-1, 1]}.
      * @param gamma Scale parameter. Must be strictly positive and finite.
      * @param delta Location parameter. Must be finite.
      */
     private static void assertConstructorThrows(double alpha, double beta, double gamma, double delta) {
         final UniformRandomProvider rng = new SplitMix64(0L);
         Assertions.assertThrows(IllegalArgumentException.class,
             () -> StableSampler.of(rng, alpha, beta, gamma, delta));
     }
 
     /**
      * Assumption test:
      * Test the limits of the value {@code tau} at the extreme limits of {@code alpha}.
      * The expression is evaluated against the original CMS algorithm. The method
      * has been updated to ensure symmetry around zero.
      *
      * <p>The test demonstrates that tau can be zero even when beta is not zero. Thus
      * the choice of a beta=0 sampler must check tau and not beta.
      */
     @Test
     void testTauLimits() {
         // At the limit of beta, tau ranges from 2/pi to 0 as alpha moves away from 1.
         final double beta = 1;
 
         // alpha -> 2: tau -> 0
         // alpha -> 0: tau -> 0
         Assertions.assertEquals(0.0, CMSStableSampler.getTau(2, beta));
         Assertions.assertEquals(0.0, CMSStableSampler.getTau(0, beta));
 
         // Full range over 0 to 2.
         for (int i = 0; i <= 512; i++) {
             // This is a power of 2 so the symmetric test uses an exact mirror
             final double alpha = (double) i / 256;
             final double tau = CMSStableSampler.getTau(alpha, beta);
             final double expected = getTauOriginal(alpha, beta);
             Assertions.assertEquals(expected, tau, 1e-15);
 
             // Symmetric
             Assertions.assertEquals(tau, CMSStableSampler.getTau(2 - alpha, beta));
         }
 
         // alpha -> 1: tau -> beta / (pi / 2) = 0.6366
         final double limit = beta / PI_2;
         Assertions.assertEquals(limit, CMSStableSampler.getTau(1, beta));
         for (double alpha : new double[] {1.01, 1 + 1e-6, 1, 1 - 1e-6, 0.99}) {
             final double tau = CMSStableSampler.getTau(alpha, beta);
             final double expected = getTauOriginal(alpha, beta);
             Assertions.assertEquals(expected, tau, 1e-15);
             // Approach the limit
             Assertions.assertEquals(limit, tau, Math.abs(1 - alpha) + 1e-15);
         }
 
         // It can be zero if beta is zero or close to zero when alpha != 1.
         // This requires we check tau==0 instead of beta==0 to switch to
         // a beta = 0 sampler.
         Assertions.assertEquals(0.0, CMSStableSampler.getTau(1.3, 0.0));
         Assertions.assertEquals(0.0, CMSStableSampler.getTau(1.5, Double.MIN_VALUE));
         Assertions.assertNotEquals(0.0, CMSStableSampler.getTau(1.0, Double.MIN_VALUE));
 
         // The sign of beta determines the sign of tau.
         Assertions.assertEquals(0.5, CMSStableSampler.getTau(1.5, beta));
         Assertions.assertEquals(0.5, CMSStableSampler.getTau(0.5, beta));
         Assertions.assertEquals(-0.5, CMSStableSampler.getTau(1.5, -beta));
         Assertions.assertEquals(-0.5, CMSStableSampler.getTau(0.5, -beta));
 
         // Check monototic at the transition point to switch to a different computation.
         final double tau1 = CMSStableSampler.getTau(Math.nextDown(1.5), 1);
         final double tau2 = CMSStableSampler.getTau(1.5, 1);
         final double tau3 = CMSStableSampler.getTau(Math.nextUp(1.5), 1);
         Assertions.assertTrue(tau1 > tau2);
         Assertions.assertTrue(tau2 > tau3);
         // Test symmetry at the transition
         Assertions.assertEquals(tau1, CMSStableSampler.getTau(2 - Math.nextDown(1.5), 1));
         Assertions.assertEquals(tau2, CMSStableSampler.getTau(0.5, 1));
         Assertions.assertEquals(tau3, CMSStableSampler.getTau(2 - Math.nextUp(1.5), 1));
     }
 
     /**
      * Gets tau using the original method from the CMS algorithm implemented in the
      * program RSTAB. This does not use {@link SpecialMath#tan2(double)} but uses
      * {@link Math#tan(double)} to implement {@code tan(x) / x}.
      *
      * @param alpha alpha
      * @param beta the beta
      * @return tau
      */
     private static double getTauOriginal(double alpha, double beta) {
         final double eps = 1 - alpha;
         // Compute RSTAB prefactor
         double tau;
 
         // Use the method from Chambers et al (1976).
         // TAN2(x) = tan(x) / x
         // PIBY2 = pi / 2
         // Comments are the FORTRAN code from the RSTAB routine.
 
         if (eps > -0.99) {
             // TAU = BPRIME / (TAN2(EPS * PIBY2) * PIBY2)
             final double tan2 = eps == 0 ? 1 : Math.tan(eps * PI_2) / (eps * PI_2);
             tau = beta / (tan2 * PI_2);
         } else {
             // TAU = BPRIME * PIBY2 * EPS * (1.-EPS) * TAN2 ((1. -EPS) * PIBY2)
             final double meps1 = 1 - eps;
             final double tan2 = Math.tan(meps1 * PI_2) / (meps1 * PI_2);
             tau = beta * PI_2 * eps * meps1 * tan2;
         }
 
         return tau;
     }
 
     /**
      * Assumption test:
      * Test the value {@code a2} is not zero. Knowing {@code a2} is not zero simplifies
      * correction of non-finite results from the CMS algorithm.
      */
     @Test
     void testA2IsNotZero() {
         // The extreme limit of the angle phiby2. This is ignored by the sampler
         // as it can result in cancellation of terms and invalid results.
         final double p0 = getU(Long.MIN_VALUE);
         Assertions.assertEquals(-PI_4, p0);
 
         // These are the limits to generate (-pi/4, pi/4)
         final double p1 = getU(Long.MIN_VALUE + (1 << 10));
         final double p2 = getU(Long.MAX_VALUE);
         Assertions.assertNotEquals(-PI_4, p1);
         Assertions.assertNotEquals(PI_4, p2);
         Assertions.assertEquals(-PI_4 + PI_4 * DU, p1);
         Assertions.assertEquals(PI_4 - PI_4 * DU, p2);
 
         for (double phiby2 : new double[] {p1, p2}) {
             // phiby2 in (-pi/4, pi/4)
             // a in (-1, 1)
             final double a = phiby2 * SpecialMath.tan2(phiby2);
             Assertions.assertEquals(Math.copySign(Math.nextDown(1.0), phiby2), a);
             final double da = a * a;
             final double a2 = 1 - da;
             // The number is close to but not equal to zero
             Assertions.assertNotEquals(0.0, a2);
             // The minimum value of a2 is 2.220E-16 = 2^-52
             Assertions.assertEquals(0x1.0p-52, a2);
         }
     }
 
     /**
      * Assumption test:
      * Test the value of the numerator used to compute z. If this is negative then
      * computation of log(z) creates a NaN. This effect occurs when the uniform
      * random deviate u is either 0 or 1 and beta is -1 or 1. The effect is reduced
      * when u is in the range {@code (0, 1)} but not eliminated. The test
      * demonstrates: (a) the requirement to check z during the sample method when
      * {@code alpha!=1}; and (b) when {@code alpha=1} then z cannot be zero when u
      * is in the open interval {@code (0, 1)}.
      */
     @Test
     void testZIsNotAlwaysAboveZero() {
         // A long is used to create phi/2:
         // The next to limit values for the phi/2
         final long x00 = Long.MIN_VALUE;
         final long x0 = Long.MIN_VALUE + (1 << 10);
         final long x1 = Long.MAX_VALUE;
         Assertions.assertEquals(-PI_4, getU(x00));
         Assertions.assertEquals(-PI_4 + DU * PI_4, getU(x0));
         Assertions.assertEquals(PI_4 - DU * PI_4, getU(x1));
         // General case numerator:
         // b2 + 2 * phiby2 * bb * tau
         // To generate 0 numerator requires:
         // b2 == -2 * phiby2 * bb * tau
         //
         // The expansion of the terms is:
         // 1 - tan^2(eps * phi/2)
         // == -phi * tan2(eps * phi/2) *        beta
         //                               -----------------------
         //                               tan2(eps * pi/2) * pi/2
         //
         // == -2 * phi * tan2(eps * phi/2) * beta
         //    -----------------------------------
         //          pi * tan2(eps * pi/2)
         //
         // == -2 * phi * tan(eps * phi/2) / (eps * phi/2) * beta
         //    --------------------------------------------------
         //          pi * tan(eps * pi/2) / (eps * pi/2)
         //
         // if phi/2 = pi/4, x = eps * phi/2:
         // == -2 * pi/4 * tan(x) / (x) * beta
         //    ---------------------------------
         //          pi * tan(2x) / (2x)
         //
         // This is a known double-angle identity for tan:
         // 1 - tan^2(x) == -2 tan(x) * beta
         //                 ---------
         //                  tan(2x)
         // Thus if |beta|=1 with opposite sign to phi, and |phi|=pi/2
         // the numerator is zero for all alpha.
         // This is not true due to floating-point
         // error but the following cases are known to exhibit the result.
 
         Assertions.assertEquals(0.0, computeNumerator(0.859375, 1, x00));
         // Even worse to have negative as the log(-ve) = nan
         Assertions.assertTrue(0.0 > computeNumerator(0.9375, 1, x00));
         Assertions.assertTrue(0.0 > computeNumerator(1.90625, 1, x00));
 
         // As phi reduces in magnitude the equality fails.
         // The numerator=0 can often be corrected
         // with the next random variate from the range limit.
         Assertions.assertTrue(0.0 < computeNumerator(0.859375, 1, x0));
         Assertions.assertTrue(0.0 < computeNumerator(0.9375, 1, x0));
         Assertions.assertTrue(0.0 < computeNumerator(1.90625, 1, x0));
 
         // WARNING:
         // Even when u is not at the limit floating point error can still create
         // a bad numerator. This is rare but shows we must still detect this edge
         // case.
         Assertions.assertTrue(0.0 > computeNumerator(0.828125, 1, x0));
         Assertions.assertTrue(0.0 > computeNumerator(1.291015625, -1, x1));
 
         // beta=0 case the numerator reduces to b2:
         // b2 = 1 - tan^2((1-alpha) * phi/2)
         // requires tan(x)=1; x=pi/4.
         // Note: tan(x) = x * SpecialMath.tan2(x) returns +/-1 for u = +/-pi/4.
         Assertions.assertEquals(-1, SpecialMath.tan2(getU(x00)) * getU(x00));
         // Using the next value in the range this is not an issue.
         // The beta=0 sampler does not have to check for z=0.
         Assertions.assertTrue(-1 < SpecialMath.tan2(getU(x0)) * getU(x0));
         Assertions.assertTrue(1 > SpecialMath.tan2(getU(x1)) * getU(x1));
         // Use alpha=2 so 1-alpha (eps) is at the limit
         final double beta = 0;
         Assertions.assertEquals(0.0, computeNumerator(2, beta, x00));
         Assertions.assertTrue(0.0 < computeNumerator(2, beta, x0));
         Assertions.assertTrue(0.0 < computeNumerator(2, beta, x1));
         Assertions.assertTrue(0.0 < computeNumerator(Math.nextDown(2), beta, x0));
         Assertions.assertTrue(0.0 < computeNumerator(Math.nextDown(2), beta, x1));
 
         // alpha=1 case the numerator reduces to:
         // 1 + 2 * phi/2 * tau
         // A zero numerator requires:
         // 2 * phiby2 * tau = -1
         //
         // tau = 2 * beta / pi
         // phiby2 = -pi / (2 * 2 * beta)
         // beta = 1 => phiby2 = -pi/4
         // beta = -1 => phiby2 = pi/4
         // The alpha=1 sampler does not have to check for z=0 if phiby2 excludes -pi/4.
         final double alpha = 1;
         Assertions.assertEquals(0.0, computeNumerator(alpha, 1, x00));
         // Next value of u computes above zero
         Assertions.assertTrue(0.0 < computeNumerator(alpha, 1, x0));
         Assertions.assertTrue(0.0 < computeNumerator(alpha, -1, x1));
         // beta < 1 => u < 0
         // beta > -1 => u > 1
         // z=0 not possible with any other beta
         Assertions.assertTrue(0.0 < computeNumerator(alpha, Math.nextUp(-1), x00));
     }
 
     /**
      * Compute the numerator value for the z coefficient in the CMS algorithm.
      *
      * @param alpha Stability parameter. Must be in range {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in range {@code [-1, 1]}.
      * @param x The random long used to generate the uniform deviate in the range {@code [-pi/4, pi/4)}.
      * @return numerator
      */
     private static double computeNumerator(double alpha, double beta, long x) {
         final double phiby2 = getU(x);
         final double eps = 1 - alpha;
         final double tau = CMSStableSampler.getTau(alpha, beta);
 
         final double bb = SpecialMath.tan2(eps * phiby2);
         final double b = eps * phiby2 * bb;
         // Compute some necessary subexpressions
         final double db = b * b;
         final double b2 = 1 - db;
         // Compute z coefficient numerator.
         return b2 + 2 * phiby2 * bb * tau;
     }
 
     /**
      * Assumption test:
      * Test the CMS algorithm can compute the value {@code d} without creating a NaN
      * when {@code z} is any non-zero finite value. When the value {@code z} is zero or infinite
      * the computation may multiply infinity by zero and create NaN.
      */
     @Test
     void testComputeDWhenZIsFiniteNonZero() {
         final double[] zs = {Double.MIN_VALUE, Double.MAX_VALUE};
 
         final double[] alphas = {2, 1.5, 1 + 1e-6, 1, 1 - 1e-6, 0.5, 0.01, 1e-10, SMALLEST_ALPHA};
         for (final double alpha : alphas) {
             // Finite z
             for (final double z : zs) {
                 // The result may be infinite, but not NaN
                 Assertions.assertNotEquals(Double.NaN, computeD(alpha, z));
             }
 
             // May be invalid with z=0 or z=inf as some combinations multiply the
             // infinity by 0 to create NaN.
 
             // When z=0, log(z) = -inf, d = d2(sign(1-alpha) * -inf) * -inf
             final double d0 = computeD(alpha, 0);
             if (alpha < 1) {
                 // d2(-inf) * -inf = 0 * -inf = NaN
                 Assertions.assertEquals(Double.NaN, d0);
             } else if (alpha == 1) {
                 // d2(0 * -inf) -> NaN
                 Assertions.assertEquals(Double.NaN, d0);
             } else {
                 // alpha > 1
                 // d2(inf) * -inf = -inf
                 Assertions.assertEquals(Double.NEGATIVE_INFINITY, d0);
             }
 
             // When z=inf, log(z) = inf, d = d2(sign(1-alpha) * inf) * inf
             final double di = computeD(alpha, Double.POSITIVE_INFINITY);
             if (alpha < 1) {
                 // d2(inf) * inf = inf
                 Assertions.assertEquals(Double.POSITIVE_INFINITY, di);
             } else if (alpha == 1) {
                 // d2(0 * inf) -> NaN
                 Assertions.assertEquals(Double.NaN, di);
             } else {
                 // alpha > 1
                 // d2(-inf) * inf = 0 * inf = NaN
                 Assertions.assertEquals(Double.NaN, di);
             }
         }
     }
 
     /**
      * Compute the {@code d} value in the CMS algorithm.
      *
      * @param alpha alpha
      * @param z z
      * @return d
      */
     private static double computeD(double alpha, double z) {
         final double alogz = Math.log(z);
         final double eps = 1 - alpha;
         final double meps1 = 1 - eps;
         return SpecialMath.d2(eps * alogz / meps1) * (alogz / meps1);
     }
 
     /**
      * Assumption test:
      * Test the sin(alpha * phi + atan(-zeta)) term can be zero.
      * This applies to the Weron formula.
      */
     @Test
     void testSinAlphaPhiMinusAtanZeta() {
         // Note sin(alpha * phi + atan(-zeta)) is zero when:
         // alpha * phi = -atan(-zeta)
         // tan(-alpha * phi) = -zeta
         //                   = beta * tan(alpha * pi / 2)
         // beta = tan(-alpha * phi) / tan(alpha * pi / 2)
         // Find a case where the result is zero...
         for (double alpha : new double[] {0.25, 0.125}) {
             for (double phi : new double[] {PI_4, PI_4 / 2}) {
                 double beta = Math.tan(-alpha * phi) / Math.tan(alpha * PI_2);
                 double zeta = -beta * Math.tan(alpha * PI_2);
                 double atanZeta = Math.atan(-zeta);
                 Assertions.assertEquals(0.0, alpha * phi + atanZeta);
             }
         }
     }
 
     /**
      * Assumption test:
      * Test the cos(phi - alpha * (phi + xi)) term is positive.
      * This applies to the Weron formula.
      */
     @Test
     void testCosPhiMinusAlphaPhiXi() {
         // This is the extreme of cos(x) that should be used
         final double cosPi2 = Math.cos(PI_2);
         // The function is symmetric
         Assertions.assertEquals(cosPi2, Math.cos(-PI_2));
         // As pi is an approximation then the cos value is not exactly 0
         Assertions.assertTrue(cosPi2 > 0);
 
         final UniformRandomProvider rng = RandomAssert.createRNG();
 
         // The term is mirrored around 1 so use extremes between 1 and 0
         final double[] alphas = {1, Math.nextDown(1), 0.99, 0.5, 0.1, 0.05, 0.01, DU};
         // Longs to generate extremes for the angle phi. This is mirrored
         // by negation is the assert method so use values to create phi in [0, pi/2).
         final long[] xs = {0, 1 << 10, Long.MIN_VALUE >>> 1, Long.MAX_VALUE};
         for (final double alpha : alphas) {
             for (final long x : xs) {
                 assertCosPhiMinusAlphaPhiXi(alpha, x);
                 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
             }
             for (int j = 0; j < 1000; j++) {
                 final long x = rng.nextLong();
                 assertCosPhiMinusAlphaPhiXi(alpha, x);
                 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
             }
         }
         // Random alpha
         for (int i = 0; i < 1000; i++) {
             final double alpha = rng.nextDouble();
             for (final long x : xs) {
                 assertCosPhiMinusAlphaPhiXi(alpha, x);
                 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
             }
             for (int j = 0; j < 1000; j++) {
                 final long x = rng.nextLong();
                 assertCosPhiMinusAlphaPhiXi(alpha, x);
                 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
             }
         }
 
         // Enumerate alpha
         for (int i = 0; i <= 1023; i++)  {
             final double alpha = (double) i / 1023;
             for (final long x : xs) {
                 assertCosPhiMinusAlphaPhiXi(alpha, x);
                 assertCosPhiMinusAlphaPhiXi(2 - alpha, x);
             }
         }
     }
 
     /**
      * Assert the cos(phi - alpha * (phi + xi)) term is positive.
      * This asserts the term (phi - alpha * (phi + xi)) is in the interval (-pi/2, pi/2) when
      * beta is the extreme of +/-1.
      *
      * @param alpha alpha
      * @param x the long used to create the uniform deviate
      */
     private static void assertCosPhiMinusAlphaPhiXi(double alpha, long x) {
         // Update for symmetry around alpha = 1
         final double eps = 1 - alpha;
         final double meps1 = 1 - eps;
 
         // zeta = -beta * tan(alpha * pi / 2)
         // xi = atan(-zeta) / alpha
         // Compute phi - alpha * (phi + xi).
         // This value must be in (-pi/2, pi/2).
         // The term expands to:
         // phi - alpha * (phi + xi)
         // = phi - alpha * phi - atan(-zeta)
         // = (1-alpha) * phi - atan(-zeta)
         // When beta = +/-1,
         // atanZeta = +/-alpha * pi/2  if alpha < 1
         // atanZeta = +/-(2-alpha) * pi/2  if alpha > 1
         // alpha=1 => always +/-pi/2
         // alpha=0,2 => always +/-phi
         // Values in between use the addition:
         // (1-alpha) * phi +/- alpha * pi/2
         // Since (1-alpha) is exact and alpha = 1 - (1-alpha) the addition
         // cannot exceed pi/2.
 
         // Avoid the round trip using tan and arctan when beta is +/- 1
         // zeta = -beta * Math.tan(alpha * pi / 2);
         // atan(-zeta) = alpha * pi / 2
 
         final double alphaPi2;
         if (meps1 > 1) {
             // Avoid calling tan outside the domain limit [-pi/2, pi/2].
             alphaPi2 = -(2 - meps1) * PI_2;
         } else {
             alphaPi2 = meps1 * PI_2;
         }
 
         // Compute eps * phi +/- alpha * pi / 2
         // Test it is in the interval (-pi/2, pi/2)
         double phi = getU(x) * 2;
         double value = eps * phi + alphaPi2;
         Assertions.assertTrue(value <= PI_2);
         Assertions.assertTrue(value >= -PI_2);
         value = eps * phi - alphaPi2;
         Assertions.assertTrue(value <= PI_2);
         Assertions.assertTrue(value >= -PI_2);
 
         // Mirror the deviate
         phi = -phi;
         value = eps * phi + alphaPi2;
         Assertions.assertTrue(value <= PI_2);
         Assertions.assertTrue(value >= -PI_2);
         value = eps * phi - alphaPi2;
         Assertions.assertTrue(value <= PI_2);
         Assertions.assertTrue(value >= -PI_2);
     }
     /**
      * Assumption test:
      * Test the sin(alpha * phi) term is only zero when phi is zero.
      * This applies to the Weron formula when {@code beta = 0}.
      */
     @Test
     void testSinAlphaPhi() {
         // Smallest non-zero phi.
         // getU creates in the domain (-pi/4, pi/4) so double the angle.
         for (final double phi : new double[] {getU(-1) * 2, getU(1 << 10) * 2}) {
             final double x = Math.sin(SMALLEST_ALPHA * phi);
             Assertions.assertNotEquals(0.0, x);
             // Value is actually:
             Assertions.assertEquals(1.9361559566769725E-32, Math.abs(x));
         }
     }
 
     /**
      * Assumption test:
      * Test functions to compute {@code (exp(x) - 1) / x}. This tests the use of
      * {@link Math#expm1(double)} and {@link Math#exp(double)} to determine if the switch
      * point to the high precision version is monotonic.
      */
     @Test
     void testExpM1() {
         // Test monotonic at the switch point
         Assertions.assertEquals(d2(0.5), d2b(0.5));
         // When positive x -> 0 the value smaller bigger.
         Assertions.assertTrue(d2(Math.nextDown(0.5)) <= d2b(0.5));
         Assertions.assertEquals(d2(-0.5), d2b(-0.5));
         // When negative x -> 0 the value gets bigger.
         Assertions.assertTrue(d2(-Math.nextDown(0.5)) >= d2b(-0.5));
         // Potentially the next power of 2 could be used based on ULP errors but
         // the switch is not monotonic.
         Assertions.assertFalse(d2(Math.nextDown(0.25)) <= d2b(0.25));
     }
 
     /**
      * This is not a test.
      *
      * <p>This outputs a report of the mean ULP difference between
      * using {@link Math#expm1(double)} and {@link Math#exp(double)} to evaluate
      * {@code (exp(x) - 1) / x}. This helps choose the switch point to avoid the computationally
      * expensive expm1 function.
      */
     //@Test
     void expm1ULPReport() {
         // Create random doubles with a given exponent. Compute the mean and max ULP difference.
         final UniformRandomProvider rng = RandomAssert.createRNG();
         // For a quicker report set to <= 2^20.
         final int size = 1 << 30;
         // Create random doubles using random bits in the 52-bit mantissa.
         final long mask = (1L << 52) - 1;
         // Note:
         // The point at which there *should* be no difference between the two is when
         // exp(x) - 1 == exp(x). This will occur at exp(x)=2^54, x = ln(2^54) = 37.43.
         Assertions.assertEquals(((double) (1L << 54)) - 1, (double) (1L << 54));
         // However since expm1 and exp are only within 1 ULP of the exact result differences
         // still occur above this threshold.
         // 2^6 = 64; 2^-4 = 0.0625
         for (int signedExp = 6; signedExp >= -4; signedExp--) {
             // The exponent must be unsigned so + 1023 to the signed exponent
             final long exp = (signedExp + 1023L) << 52;
             // Test we are creating the correct numbers
             Assertions.assertEquals(signedExp, Math.getExponent(Double.longBitsToDouble(exp)));
             Assertions.assertEquals(signedExp, Math.getExponent(Double.longBitsToDouble((-1 & mask) | exp)));
             // Get the average and max ulp
             long sum1 = 0;
             long sum2 = 0;
             long max1 = 0;
             long max2 = 0;
             for (int i = size; i-- > 0;) {
                 final long bits = rng.nextLong() & mask;
                 final double x = Double.longBitsToDouble(bits | exp);
                 final double x1 = d2(x);
                 final double x2 = d2b(x);
                 final double x1b = d2(-x);
                 final double x2b = d2b(-x);
                 final long ulp1 = Math.abs(Double.doubleToRawLongBits(x1) - Double.doubleToRawLongBits(x2));
                 final long ulp2 = Math.abs(Double.doubleToRawLongBits(x1b) - Double.doubleToRawLongBits(x2b));
                 sum1 += ulp1;
                 sum2 += ulp2;
                 if (max1 < ulp1) {
                     max1 = ulp1;
                 }
                 if (max2 < ulp2) {
                     max2 = ulp2;
                 }
             }
             // CHECKSTYLE: stop Regexp
             System.out.printf("%-6s   %2d   %-24s (%d)   %-24s (%d)%n",
                 Double.longBitsToDouble(exp), signedExp,
                 (double) sum1 / size, max1, (double) sum2 / size, max2);
             // CHECKSTYLE: resume Regexp
         }
     }
 
     /**
      * Evaluate {@code (exp(x) - 1) / x} using {@link Math#expm1(double)}.
      * For {@code x} in the range {@code [-inf, inf]} returns
      * a result in {@code [0, inf]}.
      *
      * <ul>
      * <li>For {@code x=-inf} this returns {@code 0}.
      * <li>For {@code x=0} this returns {@code 1}.
      * <li>For {@code x=inf} this returns {@code inf}.
      * </ul>
      *
      * <p> This corrects {@code 0 / 0} and {@code inf / inf} division from
      * {@code NaN} to either {@code 1} or the upper bound respectively.
      *
      * @param x value to evaluate
      * @return {@code (exp(x) - 1) / x}.
      */
     private static double d2(double x) {
         // Here we use a conditional to detect both edge cases, which are then corrected.
         final double d2 = Math.expm1(x) / x;
         if (Double.isNaN(d2)) {
             // Correct edge cases.
             if (x == 0) {
                 return 1.0;
             }
             // x must have been +infinite or NaN
             return x;
         }
         return d2;
     }
 
     /**
      * Evaluate {@code (exp(x) - 1) / x} using {@link Math#exp(double)}.
      * For {@code x} in the range {@code [-inf, inf]} returns
      * a result in {@code [0, inf]}.
      *
      * <ul>
      * <li>For {@code x=-inf} this returns {@code 0}.
      * <li>For {@code x=0} this returns {@code 1}.
      * <li>For {@code x=inf} this returns {@code inf}.
      * </ul>
      *
      * <p> This corrects {@code 0 / 0} and {@code inf / inf} division from
      * {@code NaN} to either {@code 1} or the upper bound respectively.
      *
      * @param x value to evaluate
      * @return {@code (exp(x) - 1) / x}.
      */
     private static double d2b(double x) {
         // Here we use a conditional to detect both edge cases, which are then corrected.
         final double d2 = (Math.exp(x) - 1) / x;
         if (Double.isNaN(d2)) {
             // Correct edge cases.
             if (x == 0) {
                 return 1.0;
             }
             // x must have been +infinite or NaN
             return x;
         }
         return d2;
     }
 
     /**
      * Test the special d2 function returns {@code (exp(x) - 1) / x}.
      * The limits of the function are {@code [0, inf]} and it should return 1 when x=0.
      */
     @Test
     void testD2() {
         for (final double x : new double[] {Double.MAX_VALUE, Math.log(Double.MAX_VALUE), 10, 5, 1, 0.5, 0.1, 0.05, 0.01}) {
             Assertions.assertEquals(Math.expm1(x) / x, SpecialMath.d2(x), 1e-15);
             Assertions.assertEquals(Math.expm1(-x) / -x, SpecialMath.d2(-x), 1e-15);
         }
 
         // Negative infinity computes without correction
         Assertions.assertEquals(0.0, Math.expm1(Double.NEGATIVE_INFINITY) / Double.NEGATIVE_INFINITY);
         Assertions.assertEquals(0.0, SpecialMath.d2(Double.NEGATIVE_INFINITY));
 
         // NaN is returned (i.e. no correction)
         Assertions.assertEquals(Double.NaN, SpecialMath.d2(Double.NaN));
 
         // Edge cases for z=0 or z==inf require correction
         Assertions.assertEquals(Double.NaN, Math.expm1(0) / 0.0);
         Assertions.assertEquals(Double.NaN, Math.expm1(Double.POSITIVE_INFINITY) / Double.POSITIVE_INFINITY);
         // Corrected in the special function
         Assertions.assertEquals(1.0, SpecialMath.d2(0.0));
         Assertions.assertEquals(Double.POSITIVE_INFINITY, SpecialMath.d2(Double.POSITIVE_INFINITY));
     }
 
     /**
      * Test the tan2 function returns {@code tan(x) / x}.
      */
     @Test
     void testTan2() {
         // Test the value of tan(x) when the angle is generated in the open interval (-pi/4, pi/4)
         for (final long x : new long[] {Long.MIN_VALUE + (1 << 10), Long.MAX_VALUE}) {
             final double phiby2 = getU(x);
             Assertions.assertEquals(PI_4 - DU * PI_4, Math.abs(phiby2));
             final double a = phiby2 * SpecialMath.tan2(phiby2);
             // Check this is not 1
             Assertions.assertNotEquals(1, Math.abs(a));
             Assertions.assertTrue(Math.abs(a) < 1.0);
         }
 
         // At pi/4 the function reverts to Math.tan(x) / x. Test through the transition.
         final double pi = Math.PI;
         for (final double x : new double[] {pi, pi / 2, pi / 3.99, pi / 4, pi / 4.01, pi / 8, pi / 16}) {
             final double y = Math.tan(x) / x;
             Assertions.assertEquals(y, SpecialMath.tan2(x), Math.ulp(y));
         }
 
         // Test this closely matches the JDK tan function.
         // Test uniformly between 0 and pi / 4.
         // Count the errors with the ULP difference.
         // Get max ULP and mean ULP. Do this for both tan(x) and tan(x)/x functions.
         final UniformRandomProvider rng = RandomSource.XO_RO_SHI_RO_128_PP.create(0x1647816481684L);
         int count = 0;
         long ulp = 0;
         long max = 0;
         long ulp2 = 0;
         long max2 = 0;
         for (int i = 0; i < 1000; i++) {
             final double x = rng.nextDouble() * PI_4;
             count++;
             final double tanx = Math.tan(x);
             final double tan2x = SpecialMath.tan2(x);
             // Test tan(x)
             double y = x * tan2x;
             if (y != tanx) {
                 final long u = Math.abs(Double.doubleToRawLongBits(tanx) - Double.doubleToRawLongBits(y));
                 if (max < u) {
                     max = u;
                 }
                 ulp += u;
                 // Within 4 ulp. Note tan(x) is within 1 ulp of the result. So this
                 // is max 5 ulp from the result.
                 Assertions.assertEquals(tanx, y, 4 * Math.ulp(tanx));
             }
             // Test tan(x) / x
             y = tanx / x;
             if (y != tan2x) {
                 final long u = Math.abs(Double.doubleToRawLongBits(tan2x) - Double.doubleToRawLongBits(y));
                 if (max2 < u) {
                     max2 = u;
                 }
                 ulp2 += u;
                 // Within 3 ulp.
                 Assertions.assertEquals(y, tan2x, 3 * Math.ulp(y));
             }
         }
         // Mean (max) ULP is very low
         // 2^30 random samples in [0, pi / 4)
         //          tan(x)                    tan(x) / x
         // tan4283  93436.25534446817         201185 : 68313.16171793547         128079
         // tan4288c     0.5905972588807344         4 :     0.4047176940366626         3
         Assertions.assertTrue((double) ulp / count < 0.6, "Mean ULP to tan(x) is too high");
         Assertions.assertTrue((double) ulp2 / count < 0.45, "Mean ULP to tan(x) / x is too high");
         // If the value is under 1 then the sampler will break due to cancellation errors.
         Assertions.assertEquals(1.0, SpecialMath.tan2(0.0), "Must be exact tan(x) / x at x=0");
         Assertions.assertEquals(4 / Math.PI, SpecialMath.tan2(PI_4), Math.ulp(4 / Math.PI));
         Assertions.assertEquals(1.0, PI_4 * SpecialMath.tan2(PI_4), Math.ulp(1.0));
         // If this is above 1 then the sampler will break. Test at the switch point pi/4.
         Assertions.assertTrue(1.0 >= PI_4 * SpecialMath.tan2(PI_4));
         Assertions.assertTrue(1.0 >= PI_4 * SpecialMath.tan2(Math.nextDown(PI_4)));
         // Monotonic function at the transition
         Assertions.assertTrue(SpecialMath.tan2(Math.nextUp(PI_4)) >= SpecialMath.tan2(PI_4));
     }
 
     /**
      * Assumption test:
      * Demonstrate the CMS algorithm matches the Weron formula when {@code alpha != 1}.
      * This shows the two are equivalent; they should match as the formulas are rearrangements.
      */
     @Test
     void testSamplesWithAlphaNot1() {
         // Use non-extreme parameters. beta and u are negated so use non-redundant values
         final double[] alphas = {0.3, 0.9, 1.1, 1.5};
         final double[] betas = {-1, -0.5, -0.3, 0};
         final double[] ws = {0.1, 1, 3};
         final double[] us = {0.1, 0.25, 0.5, 0.8};
 
         final double relative = 1e-5;
         final double absolute = 1e-10;
         for (final double alpha : alphas) {
             for (final double beta : betas) {
                 for (final double w : ws) {
                     for (final double u : us) {
                         final double x = sampleCMS(alpha, beta, w, u);
                         final double y = sampleWeronAlphaNot1(alpha, beta, w, u);
                         Assertions.assertEquals(x, y, Math.max(absolute, Math.abs(x) * relative));
                         // Test symmetry
                         final double z = sampleCMS(alpha, -beta, w, 1 - u);
                         Assertions.assertEquals(x, -z, 0.0);
                     }
                 }
             }
         }
     }
 
     /**
      * Assumption test:
      * Demonstrate the CMS algorithm matches the Weron formula when {@code alpha == 1}.
      * This shows the two are equivalent; they should match as the formulas are rearrangements.
      */
     @Test
     void testSamplesWithAlpha1() {
         // Use non-extreme parameters. beta and u are negated so use non-redundant values
         final double[] betas = {-1, -0.5, -0.3, 0};
         final double[] ws = {0.1, 1, 3};
         final double[] us = {0.1, 0.25, 0.5, 0.8};
 
         final double relative = 1e-5;
         final double absolute = 1e-10;
         final double alpha = 1;
         for (final double beta : betas) {
             for (final double w : ws) {
                 for (final double u : us) {
                     final double x = sampleCMS(alpha, beta, w, u);
                     final double y = sampleWeronAlpha1(beta, w, u);
                     Assertions.assertEquals(x, y, Math.max(absolute, Math.abs(x) * relative));
                     // Test symmetry
                     final double z = sampleCMS(alpha, -beta, w, 1 - u);
                     Assertions.assertEquals(x, -z, 0.0);
                 }
             }
         }
     }
 
     /**
      * Assumption test:
      * Demonstrate the CMS formula is continuous as {@code alpha -> 1}.
      * Demonstrate the Weron formula is not continuous as {@code alpha -> 1}.
      */
     @Test
     void testConvergenceWithAlphaCloseTo1() {
         final double[] betas = {-1, -0.5, 0, 0.3, 1};
         final double[] ws = {0.1, 1, 10};
         final double[] us = {0.1, 0.25, 0.5, 0.8};
         final int steps = 30;
 
         // Start with alpha not close to 0. The value 0.0625 is a power of 2 so is scaled
         // exactly by dividing by 2. With 30 steps this ranges from 2^-4 to 2^-34 leaving alpha:
         // 1.0625 -> 1.0000000000582077 or
         // 0.9375 -> 0.9999999999417923.
         for (double deltaStart : new double[] {-0.0625, 0.0625}) {
             // As alpha approaches 1 the value should approach the value when alpha=0.
             // Count the number of times it get further away with a change of alpha.
             int cmsCount = 0;
             int weronCount = 0;
 
             for (final double beta : betas) {
                 for (final double w : ws) {
                     for (final double u : us) {
                         // CMS formulas
                         double x0 = sampleCMS(1, beta, w, u);
                         Assertions.assertTrue(Double.isFinite(x0), "Target must be finite");
 
                         // Sample should approach x0 as alpha approaches 1
                         double delta = deltaStart;
                         double dx = Math.abs(x0 - sampleCMS(1 + delta, beta, w, u));
                         for (int i = 0; i < steps; i++) {
                             delta /= 2;
                             final double dx2 = Math.abs(x0 - sampleCMS(1 + delta, beta, w, u));
                             if (dx2 > dx) {
                                 cmsCount++;
                             }
                             dx = dx2;
                         }
 
                         // Weron formulas
                         x0 = sampleWeronAlpha1(beta, w, u);
                         Assertions.assertTrue(Double.isFinite(x0), "Target must be finite");
 
                         // Sample should approach x0 as alpha approaches 1
                         delta = deltaStart;
                         dx = Math.abs(x0 - sampleWeronAlphaNot1(1 + delta, beta, w, u));
                         for (int i = 0; i < steps; i++) {
                             delta /= 2;
                             final double dx2 = Math.abs(x0 - sampleWeronAlphaNot1(1 + delta, beta, w, u));
                             if (dx2 > dx) {
                                 weronCount++;
                             }
                             dx = dx2;
                         }
                     }
                 }
             }
 
             // The CMS formala monotonically converges
             Assertions.assertEquals(0, cmsCount);
             // The weron formula does not monotonically converge
             // (difference to the target can be bigger when alpha moves closer to 1).
             Assertions.assertTrue(weronCount > 200);
         }
     }
 
     /**
      * Test extreme inputs to the CMS algorithm where {@code alpha != 1} and/or
      * {@code beta != 0}. These demonstrate cases where the parameters and the
      * random variates will create non-finite samples. The test checks that the Weron
      * formula can create an appropriate sample for all cases where the CMS formula fails.
      */
     @Test
     void testExtremeInputsToSample() {
         // Demonstrate instability when w = 0
         Assertions.assertEquals(Double.NaN, sampleCMS(1.3, 0.7, 0, 0.25));
         Assertions.assertTrue(Double.isFinite(sampleCMS(1.3, 0.7, SMALL_W, 0.25)));
 
         // Demonstrate instability when u -> 0 or 1, and |beta| = 1
         Assertions.assertEquals(Double.NaN, sampleCMS(1.1, 1.0, 0.1, 0));
         Assertions.assertTrue(Double.isFinite(sampleCMS(1.1, 1.0, 0.1, DU)));
 
         // Demonstrate instability when alpha -> 0
 
         // Small alpha does not tolerate very small w.
         Assertions.assertEquals(Double.NaN, sampleCMS(0.01, 0.7, SMALL_W, 0.5));
 
         // Very small alpha does not tolerate u approaching 0 or 1 (depending on the
         // skew)
         Assertions.assertEquals(Double.NaN, sampleCMS(1e-5, 0.7, 1.0, 1e-4));
         Assertions.assertEquals(Double.NaN, sampleCMS(1e-5, -0.7, 1.0, 1 - 1e-4));
 
         final double[] alphas = {Math.nextDown(2), 1.3, 1.1, Math.nextUp(1), 1, Math.nextDown(1), 0.7, 0.1, 0.05, 0.01, 0x1.0p-16};
         final double[] betas = {1, 0.9, 0.001, 0};
         // Avoid zero for the exponential sample.
         // Test the smallest non-zero sample from the ArhensDieter exponential sampler,
         // and the largest sample.
         final double[] ws = {0, SMALL_W, 0.001, 1, 10, LARGE_W};
         // The algorithm requires a uniform deviate in (0, 1).
         // Use extremes of the 2^53 dyadic rationals in (0, 1) up to the symmetry limit
         // (i.e. 0.5).
         final double[] us = {DU, 2 * DU, 0.0001, 0.5 - DU, 0.5};
 
         int nan1 = 0;
 
         for (final double alpha : alphas) {
             for (final double beta : betas) {
                 if (alpha == 1 && beta == 0) {
                     // Ignore the Cauchy case
                     continue;
                 }
                 // Get the support of the distribution. This is not -> +/-infinity
                 // when alpha < 1 and beta = +/-1.
                 final double[] support = getSupport(alpha, beta);
                 final double lower = support[0];
                 final double upper = support[1];
                 for (final double w : ws) {
                     for (final double u : us) {
                         final double x1 = sampleCMS(alpha, beta, w, u);
                         final double x2 = sampleWeron(alpha, beta, w, u);
 
                         if (Double.isNaN(x1)) {
                             nan1++;
                         }
                         // The edge-case corrected Weron formula should not fail
                         Assertions.assertNotEquals(Double.NaN, x2);
 
                         // Check symmetry of each formula.
                         // Use a delta of zero to allow equality of 0.0 and -0.0.
                         Assertions.assertEquals(x1, 0.0 - sampleCMS(alpha, -beta, w, 1 - u), 0.0);
                         Assertions.assertEquals(x2, 0.0 - sampleWeron(alpha, -beta, w, 1 - u), 0.0);
 
                         if (Double.isInfinite(x1) && x1 != x2) {
                             // Check the Weron correction for extreme samples.
                             // The result should be at the correct *finite* support bounds.
                             // Note: This applies when alpha < 1 and beta = +/-1.
                             Assertions.assertTrue(lower <= x2 && x2 <= upper);
                         }
                     }
                 }
             }
         }
 
         // The CMS algorithm is expected to fail some cases
         Assertions.assertNotEquals(0, nan1);
     }
 
     /**
      * Create a sample from a stable distribution. This is an implementation of the CMS
      * algorithm to allow exploration of various input values. The algorithm matches that
      * in the {@link CMSStableSampler} with the exception that the uniform variate
      * is provided in {@code (0, 1)}, not{@code (-pi/4, pi/4)}.
      *
      * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
      * @param w Exponential variate
      * @param u Uniform variate
      * @return the sample
      */
     private static double sampleCMS(double alpha, double beta, double w, double u) {
         final double phiby2 = PI_2 * (u - 0.5);
         final double eps = 1 - alpha;
         // Do not use alpha in place of 1 - eps. When alpha < 0.5, 1 - eps == alpha is not
         // always true as the reverse is not exact.
         final double meps1 = 1 - eps;
 
         // Compute RSTAB prefactor
         final double tau = CMSStableSampler.getTau(alpha, beta);
 
         // Generic stable distribution that is continuous as alpha -> 1.
         // This is a trigonomic rearrangement of equation 4.1 from Chambers et al (1976)
         // as implemented in the Fortran program RSTAB.
         // Uses the special functions:
         // tan2 = tan(x) / x
         // d2 = (exp(x) - 1) / x
         // Here tan2 is implemented using an high precision approximation.
 
         // Compute some tangents
         // Limits for |phi/2| < pi/4
         // a in (-1, 1)
         final double a = phiby2 * SpecialMath.tan2(phiby2);
         // bb in [1, 4/pi)
         final double bb = SpecialMath.tan2(eps * phiby2);
         // b in (-1, 1)
         final double b = eps * phiby2 * bb;
         // Compute some necessary subexpressions
         final double da = a * a;
         final double db = b * b;
         // a2 in (0, 1]
         final double a2 = 1 - da;
         // a2p in [1, 2)
         final double a2p = 1 + da;
         // b2 in (0, 1]
         final double b2 = 1 - db;
         // b2p in [1, 2)
         final double b2p = 1 + db;
         // Compute coefficient.
         // Note:
         // Avoid z <= 0 to avoid log(z) as negative infinity or nan.
         // This occurs when |phiby2| -> +/-pi/4 and |beta| -> 1.
         // Problems:
         // numerator=0 => z=0
         // denominator=0 => z=inf
         // numerator=denominator=0 => z=nan
         // 1. w or a2 are zero so the denominator is zero. w can be a rare exponential sample.
         // a2 -> zero if the uniform deviate is 0 or 1 and angle is |pi/4|.
         // If w -> 0 and |u-0.5| -> 0.5 then the product of w * a2 can be zero.
         // 2. |eps|=1, phiby2=|pi/4| => bb=4/pi, b=1, b2=0; if tau=0 then the numerator is zero.
         // This requires beta=0.
 
         final double z = a2p * (b2 + 2 * phiby2 * bb * tau) / (w * a2 * b2p);
         // Compute the exponential-type expression
         final double alogz = Math.log(z);
         final double d = SpecialMath.d2(eps * alogz / meps1) * (alogz / meps1);
 
         // Compute stable
         return (1 + eps * d) *
                 (2 * ((a - b) * (1 + a * b) - phiby2 * tau * bb * (b * a2 - 2 * a))) /
                 (a2 * b2p) + tau * d;
     }
 
     /**
      * Create a sample from a stable distribution. This is an implementation of the Weron formula.
      * The formula has been modified when alpha != 1 to return the 0-parameterization result and
      * correct extreme samples to +/-infinity.
      *
      * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
      * @param w Exponential variate
      * @param u Uniform variate
      * @return the sample
      * @see <a href="https://doi.org/10.1016%2F0167-7152%2895%2900113-1">Weron, R (1996).
      * "On the Chambers-Mallows-Stuck method for simulating skewed stable random variables".
      * Statistics &amp; Probability Letters. 28 (2): 165–171.</a>
      */
     private static double sampleWeron(double alpha, double beta, double w, double u) {
         return alpha == 1 ? sampleWeronAlpha1(beta, w, u) : sampleWeronAlphaNot1(alpha, beta, w, u);
     }
 
     /**
      * Create a sample from a stable distribution. This is an implementation of the
      * Weron {@code alpha != 1} formula. The formula has been modified to return the
      * 0-parameterization result and correct extreme samples to +/-infinity. The
      * algorithm matches that in the {@link WeronStableSampler} with the exception that
      * the uniform variate is provided in {@code (0, 1)}, not{@code (-pi/2, pi/2)}.
      *
      * <p>Due to the increasingly large shift (up to 1e16) as {@code alpha -> 1}
      * that is used to move the result to the 0-parameterization the samples around
      * the mode of the distribution have large cancellation and a reduced number of
      * bits in the sample value.
      *
      * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
      * @param w Exponential variate
      * @param u Uniform variate
      * @return the sample
      * @see <a href="https://doi.org/10.1016%2F0167-7152%2895%2900113-1">Weron, R
      * (1996).
      * "On the Chambers-Mallows-Stuck method for simulating skewed stable random variables".
      * Statistics &amp; Probability Letters. 28 (2): 165–171.</a>
      */
     private static double sampleWeronAlphaNot1(double alpha, double beta, double w, double u) {
         // Update for symmetry around alpha = 1
         final double eps = 1 - alpha;
         final double meps1 = 1 - eps;
 
         double zeta;
         if (meps1 > 1) {
             zeta = beta * Math.tan((2 - meps1) * PI_2);
         } else {
             zeta = -beta * Math.tan(meps1 * PI_2);
         }
 
         final double scale = Math.pow(1 + zeta * zeta, 0.5 / meps1);
         final double invAlpha = 1.0 / meps1;
         final double invAlphaM1 = invAlpha - 1;
 
         final double phi = Math.PI * (u - 0.5);
 
         // Generic stable distribution.
 
         // Note: Adding back zeta creates the 0-parameterization defined in Nolan (1998):
         // X ~ S0_alpha(s,beta,u0) with s=1, u0=0 for a standard random variable.
         // As alpha -> 1 the translation zeta to create the stable deviate
         // in the 0-parameterization is increasingly large as tan(pi/2) -> infinity.
         // The max translation is approximately 1e16.
         // Without this translation the stable deviate is in the 1-parameterization
         // and the function is not continuous with respect to alpha.
         // Due to the large zeta when alpha -> 1 the number of bits of the output variable
         // are very low due to cancellation.
 
         // As alpha -> 0 or 2 then zeta -> 0 and cancellation is not relevant.
         // The formula can be modified for infinite terms to compute a result for extreme
         // deviates u and w when the CMS formula fails.
 
         // Note the following term is subject to floating point error:
         // final double xi = Math.atan(-zeta) / alpha;
         // final double alphaPhiXi = alpha * (phi + xi);
         // This is required: cos(phi - alphaPhiXi) > 0 => phi - alphaPhiXi in (-pi/2, pi/2).
         // Thus we compute atan(-zeta) and use it to compute two terms:
         // [1] alpha * (phi + xi) = alpha * (phi + atan(-zeta) / alpha) = alpha * phi + atan(-zeta)
         // [2] phi - alpha * (phi + xi) = phi - alpha * phi - atan(-zeta) = (1-alpha) * phi - atan(-zeta)
         final double atanZeta = Math.atan(-zeta);
 
         // Compute terms
         // Either term can be infinite or 0. Certain parameters compute 0 * inf.
         // t1=inf occurs alpha -> 0.
         // t1=0 occurs when beta = tan(-alpha * phi) / tan(alpha * pi / 2).
         // t2=inf occurs when w -> 0 and alpha -> 0.
         // t2=0 occurs when alpha -> 0 and phi -> pi/2.
         // Detect zeros and return as zeta.
 
         // Note sin(alpha * phi + atanZeta) is zero when:
         // alpha * phi = -atan(-zeta)
         // tan(-alpha * phi) = -zeta
         //                   = beta * tan(alpha * pi / 2)
         // Since |phi| < pi/2 this requires beta to have an opposite sign to phi
         // and a magnitude < 1. This is possible and in this case avoid a possible
         // 0 / 0 by setting the result as if term t1=0 and the result is zeta.
         double t1 = Math.sin(meps1 * phi + atanZeta);
         if (t1 == 0) {
             return zeta;
         }
         // Since cos(phi) is in (0, 1] this term will not create a
         // large magnitude to create t1 = 0.
         t1 /= Math.pow(Math.cos(phi), invAlpha);
 
         // Iff Math.cos(eps * phi - atanZeta) is zero then 0 / 0 can occur if w=0.
         // Iff Math.cos(eps * phi - atanZeta) is below zero then NaN will occur
         // in the power function. These cases are avoided by u=(0,1) and direct
         // use of arctan(-zeta).
         final double t2 = Math.pow(Math.cos(eps * phi - atanZeta) / w, invAlphaM1);
         if (t2 == 0) {
             return zeta;
         }
 
         return t1 * t2 * scale + zeta;
     }
 
     /**
      * Create a sample from a stable distribution. This is an implementation of the Weron
      * {@code alpha == 1} formula. The algorithm matches that
      * in the {@link Alpha1StableSampler} with the exception that
      * the uniform variate is provided in {@code (0, 1)}, not{@code (-pi/2, pi/2)}.
      *
      * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
      * @param w Exponential variate
      * @param u Uniform variate
      * @return the sample
      * @see <a href="https://doi.org/10.1016%2F0167-7152%2895%2900113-1">Weron, R (1996).
      * "On the Chambers-Mallows-Stuck method for simulating skewed stable random variables".
      * Statistics &amp; Probability Letters. 28 (2): 165–171.</a>
      */
     private static double sampleWeronAlpha1(double beta, double w, double u) {
         // phi in (-pi/2, pi/2)
         final double phi = Math.PI * (u - 0.5);
 
         // Generic stable distribution with alpha = 1
         final double betaPhi = PI_2 + beta * phi;
         return (betaPhi * Math.tan(phi) -
                beta * Math.log(PI_2 * w * Math.cos(phi) / betaPhi)) / PI_2;
     }
 
     /*******************************/
     /* Tests for the StableSampler */
     /*******************************/
 
     /**
      * Test the general CMS sampler when the random generator outputs create
      * deviates that cause the value {@code z} to be negative.
      */
     @Test
     void testSamplesWithZBelow0() {
         // Call the CMS algorithm with u->1; phi/2 -> pi/4.
         // The value with all bits set generates phi/2 -> pi/4.
         // Add a long to create a big value for w of 5.
         // The parameters create cancellation in the numerator of z to create a negative z.
         final long[] longs = {Long.MAX_VALUE, -6261465550279131136L};
 
         final double phiby2 = PI_4 - PI_4 * DU;
         final double w = 5.0;
         assertUWSequence(new double[] {
             phiby2, w,
         }, longs);
 
         // The alpha parameter has been identified via a search with beta=-1.
         // See testZIsNotAlwaysAboveZero()
         final double alpha = 1.291015625;
         final double beta = -1;
         Assertions.assertTrue(0.0 > computeNumerator(alpha, beta, Long.MAX_VALUE));
 
         // z will be negative. Repeat computation assumed to be performed by the sampler.
         // This ensures the test should be updated if the sampler implementation changes.
         final double eps = 1 - alpha;
         final double tau = CMSStableSampler.getTau(alpha, beta);
         final double a = phiby2 * SpecialMath.tan2(phiby2);
         final double bb = SpecialMath.tan2(eps * phiby2);
         final double b = eps * phiby2 * bb;
         final double da = a * a;
         final double db = b * b;
         final double a2 = 1 - da;
         final double a2p = 1 + da;
         final double b2 = 1 - db;
         final double b2p = 1 + db;
         final double z = a2p * (b2 + 2 * phiby2 * bb * tau) / (w * a2 * b2p);
         Assertions.assertTrue(0.0 > z);
 
         final StableSampler sampler = StableSampler.of(createRngWithSequence(longs), alpha, beta);
         // It should not be NaN or infinite
         Assertions.assertTrue(Double.isFinite(sampler.sample()), "Sampler did not recover");
     }
 
     /**
      * Test the general CMS sampler when the random generator outputs create
      * deviates that cause the value {@code z} to be infinite.
      */
     @Test
     void testSamplesWithZInfinite() {
         // Call the CMS algorithm with w=0 (and phi/2 is not extreme).
         final long[] longs = {Long.MIN_VALUE >>> 1, 0};
 
         assertUWSequence(new double[] {
             PI_4 / 2, 0,
         }, longs);
 
         for (final double alpha : new double[] {0.789, 1, 1.23}) {
             // Test all directions
             for (final double beta : new double[] {-0.56, 0, 0.56}) {
                 // Ignore Cauchy case which does not use the exponential deviate
                 if (alpha == 1 && beta == 0) {
                     continue;
                 }
                 final StableSampler sampler = StableSampler.of(createRngWithSequence(longs), alpha, beta);
                 final double x = sampler.sample();
                 // It should not be NaN
                 Assertions.assertFalse(Double.isNaN(x), "Sampler did not recover");
                 if (beta != 0) {
                     // The sample is extreme so should be at a limit of the support
                     if (alpha < 0) {
                         // Effectively +/- infinity
                         Assertions.assertEquals(Math.copySign(Double.POSITIVE_INFINITY, beta), x);
                     } else if (alpha > 1) {
                         // At the distribution mean
                         final double[] support = getSupport(alpha, beta);
                         final double mu = support[2];
                         Assertions.assertEquals(mu, x);
                     }
                 }
             }
         }
     }
 
     /**
      * Test the CMS sampler when the random generator outputs create
      * deviates that cause the value {@code d} to be infinite.
      */
     @Test
     void testSamplesWithDInfinite() {
         // beta != 0 but with low skew to allow the direction switch in
         // phi/2 to create opposite directions.
         testSamplesWithDInfinite(0.01);
         testSamplesWithDInfinite(-0.01);
     }
 
     /**
      * Test the {@code beta=0} CMS sampler when the random generator outputs create
      * deviates that cause the value {@code d} to be infinite.
      */
     @Test
     void testBeta0SamplesWithDInfinite() {
         testSamplesWithDInfinite(0.0);
     }
 
     /**
      * Test the CMS sampler when the random generator outputs create deviates that
      * cause the value {@code d} to be infinite. This applies to the general sampler
      * or the sampler with {@code beta=0}.
      *
      * @param beta beta (should be close or equal to zero to allow direction changes due to the
      * angle phi/2 to be detected)
      */
     private static void testSamplesWithDInfinite(double beta) {
         // Set-up the random deviate u to be close to -pi/4 (low), pi/4 (high) and 0.5.
         // The extreme values for u create terms during error correction that are infinite.
         final long xuLo = Long.MIN_VALUE + (1024 << 10);
         final long xuHi = Long.MAX_VALUE - (1023 << 10);
         // Call sampler with smallest possible w that is not 0. This creates a finite z
         // but an infinite d due to the use of alpha -> 0.
         final long x = 3L;
         final long[] longs = {xuLo, x, xuHi, x, 0, x};
 
         assertUWSequence(new double[] {
             -PI_4 + 1024 * DU * PI_4, SMALL_W,
             PI_4 - 1024 * DU * PI_4, SMALL_W,
             0.0, SMALL_W
         }, longs);
 
         // alpha must be small to create infinite d and beta with low skew
         // to allow the direction switch in phi/2 to create opposite directions.
         // If the skew is too large then the skew dominates the direction.
         // When u=0.5 then f=0 and a standard sum of (1 + eps * d) * f + tau * d
         // with d=inf would cause inf * 0 = NaN.
         final double alpha = 0.03;
         final StableSampler sampler = StableSampler.of(createRngWithSequence(longs), alpha, beta);
         final double x1 = sampler.sample();
         final double x2 = sampler.sample();
         final double x3 = sampler.sample();
         // Expect the limit of the support (the direction is controlled by extreme phi)
         final double max = Double.POSITIVE_INFINITY;
         Assertions.assertEquals(-max, x1);
         Assertions.assertEquals(max, x2);
         // Expect the sampler to avoid inf * 0
         Assertions.assertNotEquals(Double.NaN, x3);
         // When f=0 the sample should be in the middle (beta=0) or skewed in the direction of beta
         if (beta == 0) {
             // In the middle
             Assertions.assertEquals(0.0, x3);
         } else {
             // At the support limit
             Assertions.assertEquals(Math.copySign(max, beta), x3);
         }
     }
 
     /**
      * Test the {@code alpha=1} CMS sampler when the random generator outputs create
      * deviates that cause the value {@code phi/2} to be at the extreme limits.
      */
     @Test
     void testAlpha1SamplesWithExtremePhi() {
         // The numerator is:
         // 1 + 2 * phiby2 * tau
         // tau = beta / pi/2 when alpha=1
         //     = +/-2 / pi when alpha=1, beta = +/-1
         // This should not create zero if phi/2 is not pi/4.
         // Test the limits of phi/2 to check samples are finite.
 
         // Add a long to create an ordinary value for w of 1.0.
         // u -> -pi/4
         final long[] longs1 = {Long.MIN_VALUE + (1 << 10), 2703662416942444033L};
         assertUWSequence(new double[] {
             -PI_4 + PI_4 * DU, 1.0,
         }, longs1);
         final StableSampler sampler1 = StableSampler.of(createRngWithSequence(longs1), 1.0, 1.0);
         final double x1 = sampler1.sample();
         Assertions.assertTrue(Double.isFinite(x1), "Sampler did not recover");
 
         // u -> pi/4
         final long[] longs2 = {Long.MAX_VALUE, 2703662416942444033L};
         assertUWSequence(new double[] {
             PI_4 - PI_4 * DU, 1.0,
         }, longs2);
         final StableSampler sampler2 = StableSampler.of(createRngWithSequence(longs2), 1.0, -1.0);
         final double x2 = sampler2.sample();
         Assertions.assertTrue(Double.isFinite(x2), "Sampler did not recover");
 
         // Sample should be a reflection
         Assertions.assertEquals(x1, -x2);
     }
 
     /**
      * Test the support of the distribution when {@code gamma = 1} and
      * {@code delta = 0}. A non-infinite support applies when {@code alpha < 0} and
      * {@code |beta| = 1}.
      */
     @Test
     void testSupport() {
         testSupport(1.0, 0.0);
     }
 
     /**
      * Test the support of the distribution when {@code gamma != 1} and
      * {@code delta != 0}. A non-infinite support applies when {@code alpha < 0} and
      * {@code |beta| = 1}.
      */
     @Test
     void testSupportWithTransformation() {
         // This tests extreme values which should not create NaN results
         for (final double gamma : new double[] {0.78, 1.23, Double.MAX_VALUE, Double.MIN_VALUE}) {
             for (final double delta : new double[] {0.43, 12.34, Double.MAX_VALUE}) {
                 testSupport(gamma, delta);
                 testSupport(gamma, -delta);
             }
         }
     }
 
     /**
      * Test the support of the distribution. This applies when {@code alpha < 0} and
      * {@code |beta| = 1}.
      *
      * @param alpha Stability parameter. Must be in range {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in range {@code [-1, 1]}.
      * @param gamma Scale parameter. Must be strictly positive and finite.
      * @param delta Location parameter. Must be finite.
      */
     private static void testSupport(double gamma, double delta) {
         // When alpha is small (<=0.1) the computation becomes limited by floating-point precision.
         final double[] alphas = {2.0, 1.5, 1.0, Math.nextDown(1), 0.99, 0.75, 0.5, 0.25, 0.1, 0.01};
         for (final double alpha : alphas) {
             testSupport(alpha, 1, gamma, delta);
             testSupport(alpha, -1, gamma, delta);
         }
     }
     /**
      * Test the support of the distribution. This applies when {@code alpha < 0} and
      * {@code |beta| = 1}.
      *
      * @param alpha Stability parameter. Must be in range {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in range {@code [-1, 1]}.
      * @param gamma Scale parameter. Must be strictly positive and finite.
      * @param delta Location parameter. Must be finite.
      */
     private static void testSupport(double alpha, double beta, double gamma, double delta) {
         // This is the inclusive bounds (no infinite values)
         final double[] support = getSupport(alpha, beta);
         // Do not scale the max value. It acts as an effective infinity.
         double lower;
         if (support[0] == -Double.MAX_VALUE) {
             lower = Double.NEGATIVE_INFINITY;
         } else {
             lower = support[0] * gamma + delta;
         }
         double upper;
         if (support[1] == Double.MAX_VALUE) {
             upper = Double.POSITIVE_INFINITY;
         } else {
             upper = support[1] * gamma + delta;
         }
         // Create an RNG that will generate extreme values:
         // Here we use 4 recursions into the tail of the exponential. The large exponential
         // deviate is approximately 30.3.
         final long[] longs = new long[] {
             // Note: Add a long of Long.MIN_VALUE to test the sampler ignores this value.
             // Hits edge case for generation of phi/4 in (-pi/4, pi/4)
             Long.MIN_VALUE,
 
             // phi/2 -> -pi/4, w=0
             Long.MIN_VALUE + (1 << 10), 0,
             // phi/2 -> -pi/4, w=large
             Long.MIN_VALUE + (1 << 10), -1, -1, -1, -1, -1, -1, -1, -1, 0,
             // phi/2 -> pi/4, w=0
             Long.MAX_VALUE, 0,
             // phi/2 -> pi/4, w=large
             Long.MAX_VALUE, -1, -1, -1, -1, -1, -1, -1, -1, 0,
             // phi/2=0, w=0
             0, 0,
             // phi/2=0, w=inf
             0, -1, -1, -1, -1, -1, -1, -1, -1, 0,
 
             // Add non extreme exponential deviate to test only extreme u
             // phi/2 -> -pi/4, w=1
             Long.MIN_VALUE + (1 << 10), 2703662416942444033L,
             // phi/2 -> pi/4, w=1
             Long.MAX_VALUE, 2703662416942444033L,
             // phi/2=0, w=1
             0, 2703662416942444033L,
 
             // Add non extreme uniform deviate to test only extreme w
             // phi/2=pi/5, w=0
             Long.MIN_VALUE >> 1, 0,
             // phi/2=pi/5, w=large
             Long.MIN_VALUE >> 1, -1, -1, -1, -1, -1, -1, -1, -1, 0,
             // phi/2=pi/5, w=0
             Long.MIN_VALUE >>> 1, 0,
             // phi/2=pi/5, w=large
             Long.MIN_VALUE >>> 1, -1, -1, -1, -1, -1, -1, -1, -1, 0,
         };
 
         // Validate series
         final double phiby2low = -PI_4 + PI_4 * DU;
         final double phiby2high = PI_4 - PI_4 * DU;
         assertUWSequence(new double[] {
             phiby2low, 0,
             phiby2low, LARGE_W,
             phiby2high, 0,
             phiby2high, LARGE_W,
             0, 0,
             0, LARGE_W,
             phiby2low, 1.0,
             phiby2high, 1.0,
             0, 1.0,
             -PI_4 / 2, 0,
             -PI_4 / 2, LARGE_W,
             PI_4 / 2, 0,
             PI_4 / 2, LARGE_W,
         }, longs);
 
         final StableSampler sampler = StableSampler.of(
             createRngWithSequence(longs), alpha, beta, gamma, delta);
         for (int i = 0; i < 100; i++) {
             final double x = sampler.sample();
             if (!(lower <= x && x <= upper)) {
                 Assertions.fail(String.format("Invalid sample. alpha=%s, beta=%s, gamma=%s, delta=%s [%s, %s] x=%s",
                     alpha, beta, gamma, delta, lower, upper, x));
             }
         }
     }
 
     /**
      * Gets the support of the distribution. This returns the inclusive bounds. So exclusive
      * infinity is computed as the maximum finite value. Compute the value {@code mu} which is the
      * mean of the distribution when {@code alpha > 1}.
      *
      * <pre>
      * x in [mu, +inf)    if alpha < 1, beta = 1
      * x in (-inf, mu]    if alpha < 1, beta = -1
      * x in (-inf, -inf)  otherwise
      * </pre>
      *
      * @param alpha the alpha
      * @param beta the beta
      * @return the support ({lower, upper, mu})
      */
     private static double[] getSupport(double alpha, double beta) {
         // Convert alpha as used by the sampler
         double eps = 1 - alpha;
         double meps1 = 1 - eps;
 
         // Since pi is approximate the symmetry is lost by wrapping.
         // Keep within the domain using (2-alpha).
         double mu;
         if (alpha > 1) {
             mu = beta * Math.tan((2 - meps1) * PI_2);
         } else {
             // Special case where tan(pi/4) is not 1 (it is Math.nextDown(1.0)).
             // This is needed when testing the Levy case during sampling.
             if (alpha == 0.5) {
                 mu = -beta;
             } else {
                 mu = -beta * Math.tan(meps1 * PI_2);
             }
         }
 
         // Standard support
         double lower = -Double.MAX_VALUE;
         double upper = Double.MAX_VALUE;
         if (meps1 < 1) {
             if (beta == 1) {
                 // alpha < 0, beta = 1
                 lower = mu;
             } else if (beta == -1) {
                 // alpha < 0, beta = -1
                 upper = mu;
             }
         }
         return new double[] {lower, upper, mu};
     }
 
     /**
      * Assumption test:
      * Test the random deviates u and w can be generated by manipulating the RNG.
      */
     @Test
     void testRandomDeviatesUandW() {
         // Extremes of the uniform deviate generated using the same method as the sampler
         final double d = DU * PI_4;
         // Test in (-pi/4, pi/4)
         Assertions.assertNotEquals(-PI_4, getU(createRngWithSequence(Long.MIN_VALUE)));
         Assertions.assertEquals(-PI_4 + d, getU(createRngWithSequence(Long.MIN_VALUE + (1 << 10))));
         Assertions.assertEquals(-PI_4 / 2, getU(createRngWithSequence(Long.MIN_VALUE >> 1)));
         Assertions.assertEquals(-d, getU(createRngWithSequence(-1)));
         Assertions.assertEquals(0.0, getU(createRngWithSequence(0)));
         Assertions.assertEquals(d, getU(createRngWithSequence(1 << 10)));
         Assertions.assertEquals(PI_4 / 2, getU(createRngWithSequence(Long.MIN_VALUE >>> 1)));
         Assertions.assertEquals(PI_4 - d, getU(createRngWithSequence(Long.MAX_VALUE)));
 
         // Extremes of the exponential sampler
         Assertions.assertEquals(0, ZigguratSampler.Exponential.of(
                 createRngWithSequence(0L)).sample());
         Assertions.assertEquals(SMALL_W, ZigguratSampler.Exponential.of(
                 createRngWithSequence(3)).sample());
         Assertions.assertEquals(0.5, ZigguratSampler.Exponential.of(
                 createRngWithSequence(1446480648965178882L)).sample());
         Assertions.assertEquals(1.0, ZigguratSampler.Exponential.of(
             createRngWithSequence(2703662416942444033L)).sample());
         Assertions.assertEquals(2.5, ZigguratSampler.Exponential.of(
                 createRngWithSequence(6092639261715210240L)).sample());
         Assertions.assertEquals(5.0, ZigguratSampler.Exponential.of(
             createRngWithSequence(-6261465550279131136L)).sample());
         Assertions.assertEquals(TAIL_W, ZigguratSampler.Exponential.of(
                 createRngWithSequence(-1, -1, 0)).sample());
         Assertions.assertEquals(3 * TAIL_W, ZigguratSampler.Exponential.of(
                 createRngWithSequence(-1, -1, -1, -1, -1, -1, 0)).sample(), 1e-14);
     }
 
     /**
      * Gets a uniform random variable in {@code (-pi/4, pi/4)}.
      *
      * <p>Copied from the StableSampler for testing. In the main sampler the variable u
      * is named either {@code phi} in {@code (-pi/2, pi/2)}, or
      * {@code phiby2} in {@code (-pi/4, pi/4)}. Here we test phiby2 for the CMS algorithm.
      *
      * @return u
      */
     private static double getU(UniformRandomProvider rng) {
         final double x = getU(rng.nextLong());
         if (x == -PI_4) {
             return getU(rng);
         }
         return x;
     }
 
     /**
      * Gets a uniform random variable in {@code [-pi/4, pi/4)} from a long value.
      *
      * <p>Copied from the StableSampler for testing. In the main sampler the variable u
      * is named either {@code phi} in {@code (-pi/2, pi/2)}, or
      * {@code phiby2} in {@code (-pi/4, pi/4)}. Here we test phiby2 for the CMS algorithm.
      *
      * <p>Examples of different output where {@code d} is the gap between values of {@code phi/2}
      * and is equal to {@code pi * 2^-55 = pi/4 * 2^-53}:
      *
      * <pre>
      * Long.MIN_VALUE                  -pi/4
      * Long.MIN_VALUE + (1 << 10)      -pi/4 + d
      * Long.MIN_VALUE >> 1             -pi/5
      * -1                              -d
      * 0                               0.0
      * 1 << 10                         d
      * Long.MIN_VALUE >>> 1            pi/5
      * Long.MAX_VALUE                  pi/4 - d
      * </pre>
      *
      * @return u
      */
     private static double getU(long x) {
         return (x >> 10) * PI_4_SCALED;
     }
 
     /**
      * Creates a RNG that will return the provided output for the next double and
      * next long functions. When the sequence is complete a valid random output
      * continues.
      *
      * <p>The sampler generates (in order):
      * <ol>
      * <li>{@code phi/2} in {@code (-pi/4, pi/4)} using long values from the RNG.
      * <li>{@code w} using the {@link ZigguratSampler.Exponential}.
      * This uses a long values from the RNG.
      * </ol>
      *
      * <p>Careful control of the sequence can generate any value for {@code w} and {@code phi/2}.
      * The sampler creates a uniform deviate first, then an exponential deviate second.
      * Examples of different output where {@code d} is the gap between values of {@code phi/2}
      * and is equal to {@code pi * 2^-55 = pi/4 * 2^-53}:
      *
      * <pre>
      * longs                           phi/2               w
      * Long.MIN_VALUE                  try again [1]
      * Long.MIN_VALUE + (1 << 10)      -pi/4 + d
      * Long.MIN_VALUE >> 1             -pi/5
      * -1                              -d
      * 0                               0.0                 0
      * 1 << 10                         d
      * Long.MIN_VALUE >>> 1            pi/5
      * Long.MAX_VALUE                  pi/4 - d
      * 3                                                   6.564735882096453E-19
      * 1446480648965178882L                                0.5
      * 2703662416942444033L                                1.0
      * 6092639261715210240L                                2.5
      * -6261465550279131136L                               5.0
      * -1, -1, 0                                           7.569274694148063
      * -1L * 2n, 0                                         n * 7.569274694148063  [2]
      * </pre>
      *
      * <ol>
      * <li>When phi/2=-pi/4 the method will ignore the value and obtain another long value.
      * <li>To create a large value for the exponential sampler requires recursion. Each input
      * of 2 * -1L will add 7.569274694148063 to the total. A long of zero will stop recursion.
      * </ol>
      *
      * @param longs the initial sequence of longs
      * @return the uniform random provider
      */
     private static UniformRandomProvider createRngWithSequence(final long... longs) {
         // Note:
         // The StableSampler uniform deviate is generated from a long.
         // It is ignored if zero, a value of 1 << 11 generates the smallest value (2^-53).
         //
         // The ZigguratSampler.Exponential uses a single long value >98% of the time.
         // To create a certain value x the input y can be obtained by reversing the
         // computation of the corresponding precomputed factor X. The lowest 8 bits of y
         // choose the index i into X so must be set as the lowest bits.
         // The long is shifted right 1 before multiplying by X so this must be reversed.
         //
         // To find y to obtain the sample x use:
         // double[] X = { /* from ZigguratSampler.Exponential */ }
         // double x = 1.0; // or any other value < 7.5
         // for (int i = 0; i < X.length; i++) {
         //      // Add back the index to the lowest 8 bits.
         //      // This will work if the number is so big that the lower bits
         //      // are zerod when casting the 53-bit mantissa to a long.
         //      long y = ((long) (x / X[i]) << 1) + i;
         //      if ((y >>> 1) * X[i] == x) {
         //          // Found y!
         //      }
         // }
 
         // Start with a valid RNG.
         // This is required for nextDouble() since invoking super.nextDouble() when
         // the sequence has expired will call nextLong() and may use the intended
         // sequence of longs.
         final UniformRandomProvider rng = RandomAssert.seededRNG();
 
         // A RNG with the provided output
         return new SplitMix64(0L) {
             private int l;
 
             @Override
             public long nextLong() {
                 if (l == longs.length) {
                     return rng.nextLong();
                 }
                 return longs[l++];
             }
         };
     }
 
     /**
      * Assert the sequence of output from a uniform deviate and exponential deviate
      * created using the same method as the sampler.
      *
      * <p>The RNG is created using
      * {@link #createRngWithSequence(long[])}. See the method javadoc for
      * examples of how to generate different deviate values.
      *
      * @param expected the expected output (u1, w1, u2, w2, u3, w3, ...)
      * @param longs the initial sequence of longs
      */
     private static void assertUWSequence(double[] expected, long[] longs) {
         final UniformRandomProvider rng = createRngWithSequence(longs);
 
         // Validate series
         final SharedStateContinuousSampler exp = ZigguratSampler.Exponential.of(rng);
         for (int i = 0; i < expected.length; i += 2) {
             final int j = i / 2;
             Assertions.assertEquals(expected[i], getU(rng), () -> j + ": Incorrect u");
             if (i + 1 < expected.length) {
                 Assertions.assertEquals(expected[i + 1], exp.sample(), () -> j + ": Incorrect w");
             }
         }
     }
 
     /**
      * Test the sampler output is a continuous function of {@code alpha} and {@code beta}.
      * This test verifies the switch to the dedicated {@code alpha=1} or {@code beta=0}
      * samplers computes a continuous function of the parameters.
      */
     @Test
     void testSamplerOutputIsContinuousFunction() {
         // Test alpha passing through 1 when beta!=0 (switch to an alpha=1 sampler)
         for (final double beta : new double[] {0.5, 0.2, 0.1, 0.001}) {
             testSamplerOutputIsContinuousFunction(1 + 8096 * DU, beta, 1.0, beta, 1 - 8096 * DU, beta, 0);
             testSamplerOutputIsContinuousFunction(1 + 1024 * DU, beta, 1.0, beta, 1 - 1024 * DU, beta, 1);
             // Not perfect when alpha -> 1
             testSamplerOutputIsContinuousFunction(1 + 128 * DU, beta, 1.0, beta, 1 - 128 * DU, beta, 1);
             testSamplerOutputIsContinuousFunction(1 + 16 * DU, beta, 1.0, beta, 1 - 16 * DU, beta, 4);
             // This works with ulp=0. Either this is a lucky random seed or because the approach
             // to 1 creates equal output.
             testSamplerOutputIsContinuousFunction(1 + DU, beta, 1.0, beta, 1 - DU, beta, 0);
         }
         // Test beta passing through 0 when alpha!=1 (switch to a beta=0 sampler)
         for (final double alpha : new double[] {1.5, 1.2, 1.1, 1.001}) {
             testSamplerOutputIsContinuousFunction(alpha, 8096 * DU, alpha, 0, alpha, -8096 * DU, 0);
             testSamplerOutputIsContinuousFunction(alpha, 1024 * DU, alpha, 0, alpha, -1024 * DU, 0);
             testSamplerOutputIsContinuousFunction(alpha, 128 * DU, alpha, 0, alpha, -128 * DU, 1);
             // Not perfect when beta is very small
             testSamplerOutputIsContinuousFunction(alpha, 16 * DU, alpha, 0, alpha, -16 * DU, 64);
             testSamplerOutputIsContinuousFunction(alpha, DU, alpha, 0, alpha, -DU, 4);
         }
 
         // Note: No test for transition to the Cauchy case (alpha=1, beta=0).
         // Requires a RNG that discards output that would be used to create a exponential
         // deviate. Just create one each time a request for nextLong is performed and
         // ensure nextLong >>> 11 is not zero.
 
         // When the parameters create a special case sampler this will not work.
         // alpha passing through 0.5 when beta=1 (Levy sampler)
         // alpha -> 2 as the sampler (Gaussian sampler).
     }
 
     /**
      * Test sampler output is a continuous function of {@code alpha} and
      * {@code beta}. Create 3 samplers with the same RNG and test the middle sampler
      * computes a value between the upper and lower sampler.
      *
      * @param alpha1 lower sampler alpha
      * @param beta1 lower sampler beta
      * @param alpha2 middle sampler alpha
      * @param beta2 middle sampler beta
      * @param alpha3 upper sampler alpha
      * @param beta3 upper sampler beta
      */
     private static void testSamplerOutputIsContinuousFunction(double alpha1, double beta1,
                                                               double alpha2, double beta2,
                                                               double alpha3, double beta3,
                                                               int ulp) {
         final long seed = 0x62738468L;
         final UniformRandomProvider rng1 = RandomSource.XO_RO_SHI_RO_128_PP.create(seed);
         final UniformRandomProvider rng2 = RandomSource.XO_RO_SHI_RO_128_PP.create(seed);
         final UniformRandomProvider rng3 = RandomSource.XO_RO_SHI_RO_128_PP.create(seed);
         final StableSampler sampler1 = StableSampler.of(rng1, alpha1, beta1);
         final StableSampler sampler2 = StableSampler.of(rng2, alpha2, beta2);
         final StableSampler sampler3 = StableSampler.of(rng3, alpha3, beta3);
         final Supplier<String> msg = () -> String.format("alpha=%s, beta=%s", alpha2, beta2);
         for (int i = 0; i < 1000; i++) {
             final double x1 = sampler1.sample();
             final double x2 = sampler2.sample();
             final double x3 = sampler3.sample();
             // x2 should be in between x1 and x3
             if (x3 > x1) {
                 if (x2 > x3) {
                     // Should be the same
                     Assertions.assertEquals(x3, x2, ulp * Math.ulp(x3), msg);
                 } else if (x2 < x1) {
                     Assertions.assertEquals(x1, x2, ulp * Math.ulp(x1), msg);
                 }
             } else if (x3 < x1) {
                 if (x2 < x3) {
                     // Should be the same
                     Assertions.assertEquals(x3, x2, ulp * Math.ulp(x3), msg);
                 } else if (x2 > x1) {
                     Assertions.assertEquals(x1, x2, ulp * Math.ulp(x1), msg);
                 }
             }
         }
     }
 
     /**
      * Test the SharedStateSampler implementation for each case using a different implementation.
      */
     @Test
     void testSharedStateSampler() {
         // Gaussian case
         testSharedStateSampler(2.0, 0.0);
         // Cauchy case
         testSharedStateSampler(1.0, 0.0);
         // Levy case
         testSharedStateSampler(0.5, 1.0);
         // Hit code coverage of alpha=0.5 (Levy case) but beta != 1
         testSharedStateSampler(0.5, 0.1);
         // Beta 0 (symmetric) case
         testSharedStateSampler(1.3, 0.0);
         // Alpha 1 case
         testSharedStateSampler(1.0, 0.23);
         // Alpha close to 1
         testSharedStateSampler(Math.nextUp(1.0), 0.23);
         // General case
         testSharedStateSampler(1.3, 0.1);
         // Small alpha cases
         testSharedStateSampler(1e-5, 0.1);
         testSharedStateSampler(1e-5, 0.0);
         // Large alpha case.
         // This hits code coverage for computing tau from (1-alpha) -> -1
         testSharedStateSampler(1.99, 0.1);
     }
 
     /**
      * Test the SharedStateSampler implementation. This tests with and without the
      * {@code gamma} and {@code delta} parameters.
      *
      * @param alpha Alpha.
      * @param beta Beta.
      */
     private static void testSharedStateSampler(double alpha, double beta) {
         final UniformRandomProvider rng1 = RandomAssert.seededRNG();
         final UniformRandomProvider rng2 = RandomAssert.seededRNG();
         StableSampler sampler1 = StableSampler.of(rng1, alpha, beta);
         StableSampler sampler2 = sampler1.withUniformRandomProvider(rng2);
         RandomAssert.assertProduceSameSequence(sampler1, sampler2);
         // Test shifted
         sampler1 = StableSampler.of(rng1, alpha, beta, 1.3, 13.2);
         sampler2 = sampler1.withUniformRandomProvider(rng2);
         RandomAssert.assertProduceSameSequence(sampler1, sampler2);
     }
 
     /**
      * Test the implementation of the transformed sampler (scaled and translated).
      */
     @Test
     void testTransformedSampler() {
         // Gaussian case
         // The Gaussian case has its own scaling where the StdDev is gamma * sqrt(2).
         // (N(x) * sqrt(2)) * gamma != N(x) * (sqrt(2) * gamma)
         // Test with a delta
         testTransformedSampler(2.0, 0.0, 1);
         // Cauchy case
         testTransformedSampler(1.0, 0.0);
         // Levy case
         testTransformedSampler(0.5, 1.0);
         // Symmetric case
         testTransformedSampler(1.3, 0.0);
         // Alpha 1 case
         testTransformedSampler(1.0, 0.23);
         // Alpha close to 1
         testTransformedSampler(Math.nextUp(1.0), 0.23);
         // General case
         testTransformedSampler(1.3, 0.1);
         // Small alpha case
         testTransformedSampler(1e-5, 0.1);
         // Large alpha case.
         // This hits the case for computing tau from (1-alpha) -> -1.
         testTransformedSampler(1.99, 0.1);
     }
 
     /**
      * Test the implementation of the transformed sampler (scaled and translated).
      * The transformed output must match exactly.
      *
      * @param alpha Alpha.
      * @param beta Beta.
      */
     private static void testTransformedSampler(double alpha, double beta) {
         testTransformedSampler(alpha, beta, 0);
     }
 
     /**
      * Test the implementation of the transformed sampler (scaled and translated).
      * The transformed output must match within the provided ULP.
      *
      * @param alpha Alpha.
      * @param beta Beta.
      * @param ulp Allowed ULP difference.
      */
     private static void testTransformedSampler(double alpha, double beta, int ulp) {
         final UniformRandomProvider[] rngs = RandomAssert.createRNG(2);
         final UniformRandomProvider rng1 = rngs[0];
         final UniformRandomProvider rng2 = rngs[1];
         final double gamma = 3.4;
         final double delta = -17.3;
         final StableSampler sampler1 = StableSampler.of(rng1, alpha, beta, gamma, delta);
         final ContinuousSampler sampler2 = createTransformedSampler(rng2, alpha, beta, gamma, delta);
         if (ulp == 0) {
             RandomAssert.assertProduceSameSequence(sampler1, sampler2);
         } else {
             for (int i = 0; i < 10; i++) {
                 final double x1 = sampler1.sample();
                 final double x2 = sampler2.sample();
                 Assertions.assertEquals(x1, x2, ulp * Math.ulp(x1));
             }
         }
     }
 
     /**
      * Create a transformed sampler from a normalized sampler scaled and translated by
      * gamma and delta.
      *
      * @param rng Source of randomness.
      * @param alpha Alpha.
      * @param beta Beta.
      * @param gamma Gamma.
      * @param delta Delta.
      * @return the transformed sampler
      */
     private static ContinuousSampler createTransformedSampler(UniformRandomProvider rng,
                                                               double alpha, double beta,
                                                               final double gamma, final double delta) {
         final StableSampler delegate = StableSampler.of(rng, alpha, beta);
         return new ContinuousSampler() {
             @Override
             public double sample() {
                 return gamma * delegate.sample() + delta;
             }
         };
     }
 
     /**
      * Test symmetry when when u and beta are mirrored around 0.5 and 0 respectively.
      */
     @Test
     void testSymmetry() {
         final byte[] seed = RandomSource.KISS.createSeed();
         for (final double alpha : new double[] {1e-4, 0.78, 1, 1.23}) {
             for (final double beta : new double[] {-0.43, 0.23}) {
                 for (final double gamma : new double[] {0.78, 1, 1.23}) {
                     for (final double delta : new double[] {-0.43, 0, 0.23}) {
                         // The sampler generates u then w.
                         // If u is not -pi/4 then only a single long is used.
                         // This can be reversed around 0 by reversing the upper 54-bits.
                         // w will use 1 long only for fast lookup and then additional longs
                         // for edge of the ziggurat sampling. Fast look-up is always used
                         // when the lowest 8-bits create a value below 252.
 
                         // Use the same random source for two samplers.
                         final UniformRandomProvider rng1 = RandomSource.KISS.create(seed);
                         final UniformRandomProvider rng2 = RandomSource.KISS.create(seed);
 
                         // RNG which will not return 0 for every other long.
                         final UniformRandomProvider forward = new SplitMix64(0) {
                             private int i;
                             @Override
                             public long nextLong() {
                                 // Manipulate alternate longs
                                 if ((i++ & 0x1) == 0) {
                                     // This must not be Long.MIN_VALUE.
                                     // So set the lowest bit of the upper 54-bits.
                                     final long x = rng1.nextLong() >>> 10 | 1L;
                                     // Shift back
                                     return x << 10;
                                 }
                                 // For the exponential sample ensure the lowest 8-bits are < 252.
                                 long x;
                                 do {
                                     x = rng1.nextLong();
                                 } while ((x & 0xff) >= 252);
                                 return x;
                             }
                         };
 
                         // RNG which will not return 0 for every other long but this long is reversed.
                         final UniformRandomProvider reverse = new SplitMix64(0) {
                             private final long upper = 1L << 54;
                             private int i;
                             @Override
                             public long nextLong() {
                                 // Manipulate alternate longs
                                 if ((i++ & 0x1) == 0) {
                                     // This must not be Long.MIN_VALUE.
                                     // So set the lowest bit of the upper 54-bits.
                                     final long x = rng2.nextLong() >>> 10 | 1L;
                                     // Reverse then shift back
                                     return (upper - x) << 10;
                                 }
                                 // For the exponential sample ensure the lowest 8-bits are < 252.
                                 long x;
                                 do {
                                     x = rng2.nextLong();
                                 } while ((x & 0xff) >= 252);
                                 return x;
                             }
                         };
 
                         final StableSampler s1 = StableSampler.of(forward, alpha, beta, gamma, delta);
                         // Since mirroring applies before the shift of delta this must be negated too
                         final StableSampler s2 = StableSampler.of(reverse, alpha, -beta, gamma, -delta);
                         for (int i = 0; i < 100; i++) {
                             Assertions.assertEquals(s1.sample(), -s2.sample());
                         }
                     }
                 }
             }
         }
     }
 
     /**
      * Test symmetry for the Levy case ({@code alpha = 0.5} and {@code beta = 1}.
      */
     @Test
     void testSymmetryLevy() {
         final double alpha = 0.5;
         final double beta = 1.0;
         final byte[] seed = RandomSource.KISS.createSeed();
         final UniformRandomProvider rng1 = RandomSource.KISS.create(seed);
         final UniformRandomProvider rng2 = RandomSource.KISS.create(seed);
         for (final double gamma : new double[] {0.78, 1, 1.23}) {
             for (final double delta : new double[] {-0.43, 0, 0.23}) {
                 final StableSampler s1 = StableSampler.of(rng1, alpha, beta, gamma, delta);
                 // Since mirroring applies before the shift of delta this must be negated too
                 final StableSampler s2 = StableSampler.of(rng2, alpha, -beta, gamma, -delta);
                 for (int i = 0; i < 100; i++) {
                     Assertions.assertEquals(s1.sample(), -s2.sample());
                 }
             }
         }
     }
 
     /**
      * Test the toString method for cases not hit in the rest of the test suite.
      * This test asserts the toString method always contains the string 'stable'
      * even for parameters that create the Gaussian, Cauchy or Levy cases.
      */
     @Test
     void testToString() {
         final UniformRandomProvider rng = RandomAssert.seededRNG();
         for (final double[] p : new double[][] {
             {1.3, 0.1},
             {2.0, 0.0},
             {1.0, 0.0},
             {0.5, 1.0},
             {1e-5, 0},
             {1e-5, 0.1},
             {0.7, 0.1, 3.0, 4.5},
         }) {
             StableSampler sampler;
             if (p.length == 2) {
                 sampler = StableSampler.of(rng, p[0], p[1]);
             } else {
                 sampler = StableSampler.of(rng, p[0], p[1], p[2], p[3]);
             }
             final String s = sampler.toString().toLowerCase();
             Assertions.assertTrue(s.contains("stable"));
         }
     }
 
     /**
      * Demonstrate the CMS sampler matches the Weron sampler when {@code alpha != 1}.
      * This shows the two are equivalent; they should match as the formulas are rearrangements.
      * Avoid testing as {@code alpha -> 1} as the Weron sampler loses bits of precision in
      * the output sample.
      *
      * <p>Note: Uses direct instantiation via the package-private constructors. This avoids
      * the factory method constructor to directly select the implementation. Constructor
      * parameters are not validated.
      */
     @Test
     void testImplementationsMatch() {
         // Avoid extreme samples. Do this by manipulating the output of nextLong.
         // Generation of the random deviate u uses the top 54-bits of the long.
         // Unset a high bit to ensure getU cannot approach pi/4.
         // Set a low bit to ensure getU cannot approach -pi/4.
         final long unsetHighBit = ~(1L << 54);
         final long setLowBit = 1L << 53;
         final double hi = getU(Long.MAX_VALUE & unsetHighBit);
         final double lo = getU(Long.MIN_VALUE | setLowBit);
         // The limits are roughly pi/4 and -pi/4
         Assertions.assertEquals(PI_4, hi, 2e-3);
         Assertions.assertEquals(-PI_4, lo, 2e-3);
         Assertions.assertEquals(0.0, lo + hi, 1e-3);
 
         // Setting a bit ensure the exponential sampler cannot be zero
         final UniformRandomProvider rng = createRngWithSequence(setLowBit);
         final double w = ZigguratSampler.Exponential.of(rng).sample();
         Assertions.assertNotEquals(0.0, w);
         // This is the actual value; it is small but not extreme.
         Assertions.assertEquals(0.0036959349092519837, w);
 
         final RandomSource source = RandomSource.XO_RO_SHI_RO_128_SS;
         final long seed = 0x83762b3daf1c43L;
         final UniformRandomProvider rng1 = new SplitMix64(0L) {
             private UniformRandomProvider delegate = source.create(seed);
             @Override
             public long next() {
                 final long x = delegate.nextLong();
                 return (x & unsetHighBit) | setLowBit;
             }
         };
         final UniformRandomProvider rng2 = new SplitMix64(0L) {
             private UniformRandomProvider delegate = source.create(seed);
             @Override
             public long next() {
                 final long x = delegate.nextLong();
                 return (x & unsetHighBit) | setLowBit;
             }
         };
 
         // Not too close to alpha=1
         final double[] alphas = {0.3, 0.5, 1.2, 1.5};
         final double[] betas = {-0.5, -0.3, -0.1, 0};
 
         final double relative = 1e-5;
         final double absolute = 1e-10;
 
         for (final double alpha : alphas) {
             for (final double beta : betas) {
                 final Supplier<String> msg = () -> String.format("alpha=%s, beta=%s", alpha, beta);
                 // WARNING:
                 // Created by direct access to package-private constructor.
                 // This is for testing only as these do not validate the parameters.
                 StableSampler s1;
                 StableSampler s2;
                 if (beta == 0) {
                     s1 = new Beta0CMSStableSampler(rng1, alpha);
                     s2 = new Beta0WeronStableSampler(rng2, alpha);
                 } else {
                     s1 = new CMSStableSampler(rng1, alpha, beta);
                     s2 = new WeronStableSampler(rng2, alpha, beta);
                 }
                 for (int i = 0; i < 1000; i++) {
                     final double x = s1.sample();
                     final double y = s2.sample();
                     Assertions.assertEquals(x, y, Math.max(absolute, Math.abs(x) * relative), msg);
                 }
             }
         }
     }
 
     /**
      * Demonstrate the general CMS sampler matches the {@code beta = 0} sampler.
      * The {@code beta = 0} sampler implements the same algorithm with cancelled terms removed.
      *
      * <p>Note: Uses direct instantiation via the package-private constructors. This avoids
      * the factory method constructor to directly select the implementation. Constructor
      * parameters are not validated.
      */
     @Test
     void testSpecializedBeta0CMSImplementation() {
         final RandomSource source = RandomSource.XO_RO_SHI_RO_128_SS;
         // Should be robust to any seed
         final byte[] seed = source.createSeed();
         final UniformRandomProvider rng1 = source.create(seed);
         final UniformRandomProvider rng2 = source.create(seed);
 
         final double[] alphas = {0.3, 0.5, 1.2, 1.5};
         for (final double alpha : alphas) {
             // WARNING:
             // Created by direct access to package-private constructor.
             // This is for testing only as these do not validate the parameters.
             final StableSampler sampler1 = new CMSStableSampler(rng1, alpha, 0.0);
             final StableSampler sampler2 = new Beta0CMSStableSampler(rng2, alpha);
             RandomAssert.assertProduceSameSequence(sampler1, sampler2);
         }
     }
 
     /**
      * Demonstrate the general Weron sampler matches the {@code beta = 0} sampler.
      * The {@code beta = 0} sampler implements the same algorithm with cancelled terms removed.
      *
      * <p>Note: Uses direct instantiation via the package-private constructors. This avoids
      * the factory method constructor to directly select the implementation. Constructor
      * parameters are not validated.
      */
     @Test
     void testSpecializedBeta0WeronImplementation() {
         final RandomSource source = RandomSource.XO_RO_SHI_RO_128_SS;
         // Should be robust to any seed
         final byte[] seed = source.createSeed();
         final UniformRandomProvider rng1 = source.create(seed);
         final UniformRandomProvider rng2 = source.create(seed);
 
         final double[] alphas = {0.3, 0.5, 1.2, 1.5};
         for (final double alpha : alphas) {
             // WARNING:
             // Created by direct access to package-private constructor.
             // This is for testing only as these do not validate the parameters.
             final StableSampler sampler1 = new WeronStableSampler(rng1, alpha, 0.0);
             final StableSampler sampler2 = new Beta0WeronStableSampler(rng2, alpha);
             RandomAssert.assertProduceSameSequence(sampler1, sampler2);
         }
     }
 
     /**
      * Test the Weron sampler when the term t1 is zero in the numerator.
      * This hits an edge case where sin(alpha * phi + atan(-zeta)) is zero.
      *
      * @see #testSinAlphaPhiMinusAtanZeta()
      */
     @Test
     void testWeronImplementationEdgeCase() {
         double alpha = 0.25;
         // Solved in testSinAlphaPhiMinusAtanZeta()
         double beta = -0.48021693505171;
         // Require phi = PI_4.
         // This is the equivalent of phi/2 = pi/5
         final long x = Long.MIN_VALUE >>> 1;
         final long[] longs = new long[] {
             // phi/2=pi/5, w=0
             x, 0,
             // phi/2=pi/5, w=large
             x, -1, -1, -1, -1, -1, -1, -1, -1, 0,
             // phi/2=pi/5, w=1
             x, 2703662416942444033L,
         };
 
         // Validate series
         assertUWSequence(new double[] {
             PI_4 / 2, 0,
             PI_4 / 2, LARGE_W,
             PI_4 / 2, 1.0,
         }, longs);
 
         final double zeta = -beta * Math.tan(alpha * PI_2);
         Assertions.assertEquals(0.0, alpha * PI_4 + Math.atan(-zeta));
 
         final UniformRandomProvider rng = createRngWithSequence(longs);
         final StableSampler sampler = new WeronStableSampler(rng, alpha, beta);
         // zeta is the offset used to shift the 1-parameterization to the
         // 0-parameterization. This is returned when other terms multiply to zero.
         Assertions.assertEquals(zeta, sampler.sample());
         Assertions.assertEquals(zeta, sampler.sample());
         Assertions.assertEquals(zeta, sampler.sample());
     }
 }