/*
    * 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 org.apache.commons.rng.UniformRandomProvider;
  
  /**
   * Samples from a stable distribution.
   *
   * <p>Several different parameterizations exist for the stable distribution.
   * This sampler uses the 0-parameterization distribution described in Nolan (2020) "Univariate Stable
   * Distributions: Models for Heavy Tailed Data". Springer Series in Operations Research and
   * Financial Engineering. Springer. Sections 1.7 and 3.3.3.
   *
   * <p>The random variable \( X \) has
   * the stable distribution \( S(\alpha, \beta, \gamma, \delta; 0) \) if its characteristic
   * function is given by:
   *
   * <p>\[ E(e^{iuX}) = \begin{cases} \exp \left (- \gamma^\alpha |u|^\alpha \left [1 - i \beta (\tan \frac{\pi \alpha}{2})(\text{sgn}(u)) \right ] + i \delta u \right ) &amp; \alpha \neq 1 \\
   * \exp \left (- \gamma |u| \left [1 + i \beta \frac{2}{\pi} (\text{sgn}(u)) \log |u| \right ] + i \delta u \right ) &amp; \alpha = 1 \end{cases} \]
   *
   * <p>The function is continuous with respect to all the parameters; the parameters \( \alpha \)
   * and \( \beta \) determine the shape and the parameters \( \gamma \) and \( \delta \) determine
   * the scale and location. The support of the distribution is:
   *
   * <p>\[ \text{support} f(x|\alpha,\beta,\gamma,\delta; 0) = \begin{cases} [\delta - \gamma \tan \frac{\pi \alpha}{2}, \infty) &amp; \alpha \lt 1\ and\ \beta = 1 \\
   * (-\infty, \delta + \gamma \tan \frac{\pi \alpha}{2}] &amp; \alpha \lt 1\ and\ \beta = -1 \\
   * (-\infty, \infty) &amp; otherwise \end{cases} \]
   *
   * <p>The implementation uses the Chambers-Mallows-Stuck (CMS) method as described in:
   * <ul>
   *  <li>Chambers, Mallows &amp; Stuck (1976) "A Method for Simulating Stable Random Variables".
   *      Journal of the American Statistical Association. 71 (354): 340–344.
   *  <li>Weron (1996) "On the Chambers-Mallows-Stuck method for simulating skewed stable
   *      random variables". Statistics &amp; Probability Letters. 28 (2): 165–171.
   * </ul>
   *
   * @see <a href="https://en.wikipedia.org/wiki/Stable_distribution">Stable distribution (Wikipedia)</a>
   * @see <a href="https://link.springer.com/book/10.1007/978-3-030-52915-4">Nolan (2020) Univariate Stable Distributions</a>
   * @see <a href="https://doi.org/10.1080%2F01621459.1976.10480344">Chambers et al (1976) JOASA 71: 340-344</a>
   * @see <a href="https://doi.org/10.1016%2F0167-7152%2895%2900113-1">Weron (1996).
   * Statistics &amp; Probability Letters. 28 (2): 165–171.</a>
   * @since 1.4
   */
  public abstract class StableSampler implements SharedStateContinuousSampler {
      /** pi / 2. */
      private static final double PI_2 = Math.PI / 2;
      /** The alpha value for the Gaussian case. */
      private static final double ALPHA_GAUSSIAN = 2;
      /** The alpha value for the Cauchy case. */
      private static final double ALPHA_CAUCHY = 1;
      /** The alpha value for the Levy case. */
      private static final double ALPHA_LEVY = 0.5;
      /** The alpha value for the {@code alpha -> 0} to switch to using the Weron formula.
       * Note that small alpha requires robust correction of infinite samples. */
      private static final double ALPHA_SMALL = 0.02;
      /** The beta value for the Levy case. */
      private static final double BETA_LEVY = 1.0;
      /** The gamma value for the normalized case. */
      private static final double GAMMA_1 = 1.0;
      /** The delta value for the normalized case. */
      private static final double DELTA_0 = 0.0;
      /** The tau value for zero. When tau is zero, this is effectively {@code beta = 0}. */
      private static final double TAU_ZERO = 0.0;
      /**
       * The lower support for the distribution.
       * This is the lower bound of {@code (-inf, +inf)}
       * If the sample is not within this bound ({@code lower < x}) then it is either
       * infinite or NaN and the result should be checked.
       */
      private static final double LOWER = Double.NEGATIVE_INFINITY;
      /**
       * The upper support for the distribution.
       * This is the upper bound of {@code (-inf, +inf)}.
       * If the sample is not within this bound ({@code x < upper}) then it is either
       * infinite or NaN and the result should be checked.
       */
      private static final double UPPER = Double.POSITIVE_INFINITY;
  
      /** Underlying source of randomness. */
      private final UniformRandomProvider rng;
  
      // Implementation notes
      //
      // The Chambers-Mallows-Stuck (CMS) method uses a uniform deviate u in (0, 1) and an
     // exponential deviate w to compute a stable deviate. Chambers et al (1976) published
     // a formula for alpha = 1 and alpha != 1. The function is discontinuous at alpha = 1
     // and to address this a trigonmoic rearrangement was provided using half angles that
     // is continuous with respect to alpha. The original discontinuous formulas were proven
     // in Weron (1996). The CMS rearrangement creates a deviate in the 0-parameterization
     // defined by Nolan (2020); the original discontinuous functions create a deviate in the
     // 1-parameterization defined by Nolan. A shift can be used to convert one parameterisation
     // to the other. The shift is the magnitude of the zeta term from the 1-parameterisation.
     // The following table shows how the zeta term -> inf when alpha -> 1 for
     // different beta (hence the discontinuity in the function):
     //
     // Zeta
     //             Beta
     // Alpha       1.0         0.5         0.25        0.1         0.0
     // 0.001       0.001571    0.0007854   0.0003927   0.0001571   0.0
     // 0.01        0.01571     0.007855    0.003927    0.001571    0.0
     // 0.05        0.07870     0.03935     0.01968     0.007870    0.0
     // 0.01        0.01571     0.007855    0.003927    0.001571    0.0
     // 0.1         0.1584      0.07919     0.03960     0.01584     0.0
     // 0.5         1.000       0.5000      0.2500      0.1000      0.0
     // 0.9         6.314       3.157       1.578       0.6314      0.0
     // 0.95        12.71       6.353       3.177       1.271       0.0
     // 0.99        63.66       31.83       15.91       6.366       0.0
     // 0.995       127.3       63.66       31.83       12.73       0.0
     // 0.999       636.6       318.3       159.2       63.66       0.0
     // 0.9995      1273        636.6       318.3       127.3       0.0
     // 0.9999      6366        3183        1592        636.6       0.0
     // 1.0         1.633E+16   8.166E+15   4.083E+15   1.633E+15   0.0
     //
     // For numerical simulation the 0-parameterization is favoured as it is continuous
     // with respect to all the parameters. When approaching alpha = 1 the large magnitude
     // of the zeta term used to shift the 1-parameterization results in cancellation and the
     // number of bits of the output sample is effected. This sampler uses the CMS method with
     // the continuous function as the base for the implementation. However it is not suitable
     // for all values of alpha and beta.
     //
     // The method computes a value log(z) with z in the interval (0, inf). When z is 0 or infinite
     // the computation can return invalid results. The open bound for the deviate u avoids
     // generating an extreme value that results in cancellation, z=0 and an invalid expression.
     // However due to floating point error this can occur
     // when u is close to 0 or 1, and beta is -1 or 1. Thus it is not enough to create
     // u by avoiding 0 or 1 and further checks are required.
     // The division by the deviate w also results in an invalid expression as the term z becomes
     // infinite as w -> 0. It should be noted that such events are extremely rare
     // (frequency in the 1 in 10^15), or will not occur at all depending on the parameters alpha
     // and beta.
     //
     // When alpha -> 0 then the distribution is extremely long tailed and the expression
     // using log(z) often computes infinity. Certain parameters can create NaN due to
     // 0 / 0, 0 * inf, or inf - inf. Thus the implementation must check the final result
     // and perform a correction if required, or generate another sample.
     // Correcting the original CMS formula has many edge cases depending on parameters. The
     // alternative formula provided by Weron is easier to correct when infinite values are
     // created. This correction is made easier by knowing that u is not 0 or 1 as certain
     // conditions on the intermediate terms can be eliminated. The implementation
     // thus generates u in the open interval (0,1) but leaves w unchecked and potentially 0.
     // The sample is generated and the result tested against the expected support. This detects
     // any NaN and infinite values. Incorrect samples due to the inability to compute log(z)
     // (extremely rare) and samples where alpha -> 0 has resulted in an infinite expression
     // for the value d are corrected using the Weron formula and returned within the support.
     //
     // The CMS algorithm is continuous for the parameters. However when alpha=1 or beta=0
     // many terms cancel and these cases are handled with specialised implementations.
     // The beta=0 case implements the same CMS algorithm with certain terms eliminated.
     // Correction uses the alternative Weron formula. When alpha=1 the CMS algorithm can
     // be corrected from infinite cases due to assumptions on the intermediate terms.
     //
     // The following table show the failure frequency (result not finite or, when beta=+/-1,
     // within the support) for the CMS algorithm computed using 2^30 random deviates.
     //
     // CMS failure rate
     //             Beta
     // Alpha       1.0         0.5         0.25        0.1         0.0
     // 1.999       0.0         0.0         0.0         0.0         0.0
     // 1.99        0.0         0.0         0.0         0.0         0.0
     // 1.9         0.0         0.0         0.0         0.0         0.0
     // 1.5         0.0         0.0         0.0         0.0         0.0
     // 1.1         0.0         0.0         0.0         0.0         0.0
     // 1.0         0.0         0.0         0.0         0.0         0.0
     // 0.9         0.0         0.0         0.0         0.0         0.0
     // 0.5         0.0         0.0         0.0         0.0         0.0
     // 0.25        0.0         0.0         0.0         0.0         0.0
     // 0.1         0.0         0.0         0.0         0.0         0.0
     // 0.05        0.0003458   0.0         0.0         0.0         0.0
     // 0.02        0.009028    6.938E-7    7.180E-7    7.320E-7    6.873E-7
     // 0.01        0.004878    0.0008555   0.0008553   0.0008554   0.0008570
     // 0.005       0.1519      0.02896     0.02897     0.02897     0.02897
     // 0.001       0.6038      0.3903      0.3903      0.3903      0.3903
     //
     // The sampler switches to using the error checked Weron implementation when alpha < 0.02.
     // Unit tests demonstrate the two samplers (CMS or Weron) product the same result within
     // a tolerance. The switch point is based on a consistent failure rate above 1 in a million.
     // At this point zeta is small and cancellation leading to loss of bits in the sample is
     // minimal.
     //
     // In common use the sampler will not have a measurable failure rate. The output will
     // be continuous as alpha -> 1 and beta -> 0. The evaluated function produces symmetric
     // samples when u and beta are mirrored around 0.5 and 0 respectively. To achieve this
     // the computation of certain parameters has been changed from the original implementation
     // to avoid evaluating Math.tan outside the interval (-pi/2, pi/2).
     //
     // Note: Chambers et al (1976) use an approximation to tan(x) / x in the RSTAB routine.
     // A JMH performance test is available in the RNG examples module comparing Math.tan
     // with various approximations. The functions are faster than Math.tan(x) / x.
     // This implementation uses a higher accuracy approximation than the original RSTAB
     // implementation; it has a mean ULP difference to Math.tan of less than 1 and has
     // a noticeable performance gain.
 
     /**
      * Base class for implementations of a stable distribution that requires an exponential
      * random deviate.
      */
     private abstract static class BaseStableSampler extends StableSampler {
         /** pi/2 scaled by 2^-53. */
         private static final double PI_2_SCALED = 0x1.0p-54 * Math.PI;
         /** pi/4 scaled by 2^-53. */
         private static final double PI_4_SCALED = 0x1.0p-55 * Math.PI;
         /** -pi / 2. */
         private static final double NEG_PI_2 = -Math.PI / 2;
         /** -pi / 4. */
         private static final double NEG_PI_4 = -Math.PI / 4;
 
         /** The exponential sampler. */
         private final ContinuousSampler expSampler;
 
         /**
          * @param rng Underlying source of randomness
          */
         BaseStableSampler(UniformRandomProvider rng) {
             super(rng);
             expSampler = ZigguratSampler.Exponential.of(rng);
         }
 
         /**
          * Gets a random value for the omega parameter ({@code w}).
          * This is an exponential random variable with mean 1.
          *
          * <p>Warning: For simplicity this does not check the variate is not 0.
          * The calling CMS algorithm should detect and handle incorrect samples as a result
          * of this unlikely edge case.
          *
          * @return omega
          */
         double getOmega() {
             // Note: Ideally this should not have a value of 0 as the CMS algorithm divides
             // by w and it creates infinity. This can result in NaN output.
             // Under certain parameterizations non-zero small w also creates NaN output.
             // Thus output should be checked regardless.
             return expSampler.sample();
         }
 
         /**
          * Gets a random value for the phi parameter.
          * This is a uniform random variable in {@code (-pi/2, pi/2)}.
          *
          * @return phi
          */
         double getPhi() {
             // See getPhiBy2 for method details.
             final double x = (nextLong() >> 10) * PI_2_SCALED;
             // Deliberate floating-point equality check
             if (x == NEG_PI_2) {
                 return getPhi();
             }
             return x;
         }
 
         /**
          * Gets a random value for the phi parameter divided by 2.
          * This is a uniform random variable in {@code (-pi/4, pi/4)}.
          *
          * <p>Note: Ideally this should not have a value of -pi/4 or pi/4 as the CMS algorithm
          * can generate infinite values when the phi/2 uniform deviate is +/-pi/4. This
          * can result in NaN output. Under certain parameterizations phi/2 close to the limits
          * also create NaN output. Thus output should be checked regardless. Avoiding
          * the extreme values simplifies the number of checks that are required.
          *
          * @return phi / 2
          */
         double getPhiBy2() {
             // As per o.a.c.rng.core.utils.NumberFactory.makeDouble(long) but using a
             // signed shift of 10 in place of an unsigned shift of 11. With a factor of 2^-53
             // this would produce a double in [-1, 1).
             // Here the multiplication factor incorporates pi/4 to avoid a separate
             // multiplication.
             final double x = (nextLong() >> 10) * PI_4_SCALED;
             // Deliberate floating-point equality check
             if (x == NEG_PI_4) {
                 // Sample again using recursion.
                 // A stack overflow due to a broken RNG will eventually occur
                 // rather than the alternative which is an infinite loop
                 // while x == -pi/4.
                 return getPhiBy2();
             }
             return x;
         }
     }
 
     /**
      * Class for implementations of a stable distribution transformed by scale and location.
      */
     private static final class TransformedStableSampler extends StableSampler {
         /** Underlying normalized stable sampler. */
         private final StableSampler sampler;
         /** The scale parameter. */
         private final double gamma;
         /** The location parameter. */
         private final double delta;
 
         /**
          * @param sampler Normalized stable sampler.
          * @param gamma Scale parameter. Must be strictly positive.
          * @param delta Location parameter.
          */
         TransformedStableSampler(StableSampler sampler, double gamma, double delta) {
             // No RNG required
             super(null);
             this.sampler = sampler;
             this.gamma = gamma;
             this.delta = delta;
         }
 
         @Override
         public double sample() {
             return gamma * sampler.sample() + delta;
         }
 
         @Override
         public StableSampler withUniformRandomProvider(UniformRandomProvider rng) {
             return new TransformedStableSampler(sampler.withUniformRandomProvider(rng),
                                                 gamma, delta);
         }
 
         @Override
         public String toString() {
             // Avoid a null pointer from the unset RNG instance in the parent class
             return sampler.toString();
         }
     }
 
     /**
      * Implement the {@code alpha = 2} stable distribution case (Gaussian distribution).
      */
     private static final class GaussianStableSampler extends StableSampler {
         /** sqrt(2). */
         private static final double ROOT_2 = Math.sqrt(2);
 
         /** Underlying normalized Gaussian sampler. */
         private final NormalizedGaussianSampler sampler;
         /** The standard deviation. */
         private final double stdDev;
         /** The mean. */
         private final double mean;
 
         /**
          * @param rng Underlying source of randomness
          * @param gamma Scale parameter. Must be strictly positive.
          * @param delta Location parameter.
          */
         GaussianStableSampler(UniformRandomProvider rng, double gamma, double delta) {
             super(rng);
             this.sampler = ZigguratSampler.NormalizedGaussian.of(rng);
             // A standardized stable sampler with alpha=2 has variance 2.
             // Set the standard deviation as sqrt(2) * scale.
             // Avoid this being infinity to avoid inf * 0 in the sample
             this.stdDev = Math.min(Double.MAX_VALUE, ROOT_2 * gamma);
             this.mean = delta;
         }
 
         /**
          * @param rng Underlying source of randomness
          * @param source Source to copy.
          */
         GaussianStableSampler(UniformRandomProvider rng, GaussianStableSampler source) {
             super(rng);
             this.sampler = ZigguratSampler.NormalizedGaussian.of(rng);
             this.stdDev = source.stdDev;
             this.mean = source.mean;
         }
 
         @Override
         public double sample() {
             return stdDev * sampler.sample() + mean;
         }
 
         @Override
         public GaussianStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
             return new GaussianStableSampler(rng, this);
         }
     }
 
     /**
      * Implement the {@code alpha = 1} and {@code beta = 0} stable distribution case
      * (Cauchy distribution).
      */
     private static final class CauchyStableSampler extends BaseStableSampler {
         /** The scale parameter. */
         private final double gamma;
         /** The location parameter. */
         private final double delta;
 
         /**
          * @param rng Underlying source of randomness
          * @param gamma Scale parameter. Must be strictly positive.
          * @param delta Location parameter.
          */
         CauchyStableSampler(UniformRandomProvider rng, double gamma, double delta) {
             super(rng);
             this.gamma = gamma;
             this.delta = delta;
         }
 
         /**
          * @param rng Underlying source of randomness
          * @param source Source to copy.
          */
         CauchyStableSampler(UniformRandomProvider rng, CauchyStableSampler source) {
             super(rng);
             this.gamma = source.gamma;
             this.delta = source.delta;
         }
 
         @Override
         public double sample() {
             // Note:
             // The CMS beta=0 with alpha=1 sampler reduces to:
             // S = 2 * a / a2, with a = tan(x), a2 = 1 - a^2, x = phi/2
             // This is a double angle identity for tan:
             // 2 * tan(x) / (1 - tan^2(x)) = tan(2x)
             // Here we use the double angle identity for consistency with the other samplers.
             final double phiby2 = getPhiBy2();
             final double a = phiby2 * SpecialMath.tan2(phiby2);
             final double a2 = 1 - a * a;
             final double x = 2 * a / a2;
             return gamma * x + delta;
         }
 
         @Override
         public CauchyStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
             return new CauchyStableSampler(rng, this);
         }
     }
 
     /**
      * Implement the {@code alpha = 0.5} and {@code beta = 1} stable distribution case
      * (Levy distribution).
      *
      * Note: This sampler can be used to output the symmetric case when
      * {@code beta = -1} by negating {@code gamma}.
      */
     private static final class LevyStableSampler extends StableSampler {
         /** Underlying normalized Gaussian sampler. */
         private final NormalizedGaussianSampler sampler;
         /** The scale parameter. */
         private final double gamma;
         /** The location parameter. */
         private final double delta;
 
         /**
          * @param rng Underlying source of randomness
          * @param gamma Scale parameter. Must be strictly positive.
          * @param delta Location parameter.
          */
         LevyStableSampler(UniformRandomProvider rng, double gamma, double delta) {
             super(rng);
             this.sampler = ZigguratSampler.NormalizedGaussian.of(rng);
             this.gamma = gamma;
             this.delta = delta;
         }
 
         /**
          * @param rng Underlying source of randomness
          * @param source Source to copy.
          */
         LevyStableSampler(UniformRandomProvider rng, LevyStableSampler source) {
             super(rng);
             this.sampler = ZigguratSampler.NormalizedGaussian.of(rng);
             this.gamma = source.gamma;
             this.delta = source.delta;
         }
 
         @Override
         public double sample() {
             // Levy(Z) = 1 / N(0,1)^2, where N(0,1) is a standard normalized variate
             final double norm = sampler.sample();
             // Here we must transform from the 1-parameterization to the 0-parameterization.
             // This is a shift of -beta * tan(pi * alpha / 2) = -1 when alpha=0.5, beta=1.
             final double z = (1.0 / (norm * norm)) - 1.0;
             // In the 0-parameterization the scale and location are a linear transform.
             return gamma * z + delta;
         }
 
         @Override
         public LevyStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
             return new LevyStableSampler(rng, this);
         }
     }
 
     /**
      * Implement the generic stable distribution case: {@code alpha < 2} and
      * {@code beta != 0}. This routine assumes {@code alpha != 1}.
      *
      * <p>Implements the Chambers-Mallows-Stuck (CMS) method using the
      * formula provided in Weron (1996) "On the Chambers-Mallows-Stuck method for
      * simulating skewed stable random variables" Statistics &amp; Probability
      * Letters. 28 (2): 165–171. This method is easier to correct from infinite and
      * NaN results by boxing intermediate infinite values.
      *
      * <p>The formula produces a stable deviate from the 1-parameterization that is
      * discontinuous at {@code alpha=1}. A shift is used to create the 0-parameterization.
      * This shift is very large as {@code alpha -> 1} and the output loses bits of precision
      * in the deviate due to cancellation. It is not recommended to use this sampler when
      * {@code alpha -> 1} except for edge case correction.
      *
      * <p>This produces non-NaN output for all parameters alpha, beta, u and w with
      * the correct orientation for extremes of the distribution support.
      * The formulas used are symmetric with regard to beta and u.
      *
      * @see <a href="https://doi.org/10.1016%2F0167-7152%2895%2900113-1">Weron, R
      * (1996). Statistics &amp; Probability Letters. 28 (2): 165–171.</a>
      */
     static class WeronStableSampler extends BaseStableSampler {
         /** Epsilon (1 - alpha). */
         protected final double eps;
         /** Alpha (1 - eps). */
         protected final double meps1;
         /** Cache of expression value used in generation. */
         protected final double zeta;
         /** Cache of expression value used in generation. */
         protected final double atanZeta;
         /** Cache of expression value used in generation. */
         protected final double scale;
         /** 1 / alpha = 1 / (1 - eps). */
         protected final double inv1mEps;
         /** (1 / alpha) - 1 = eps / (1 - eps). */
         protected final double epsDiv1mEps;
         /** The inclusive lower support for the distribution. */
         protected final double lower;
         /** The inclusive upper support for the distribution. */
         protected final double upper;
 
         /**
          * @param rng Underlying source of randomness
          * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
          * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
          */
         WeronStableSampler(UniformRandomProvider rng, double alpha, double beta) {
             super(rng);
             eps = 1 - alpha;
             // When alpha < 0.5, 1 - eps == alpha is not always true as the reverse is not exact.
             // Here we store 1 - eps in place of alpha. Thus eps + (1 - eps) = 1.
             meps1 = 1 - eps;
 
             // Compute pre-factors for the Weron formula used during error correction.
             if (meps1 > 1) {
                 // Avoid calling tan outside the domain limit [-pi/2, pi/2].
                 zeta = beta * Math.tan((2 - meps1) * PI_2);
             } else {
                 zeta = -beta * Math.tan(meps1 * PI_2);
             }
 
             // Do not store xi = Math.atan(-zeta) / meps1 due to floating-point division errors.
             // Directly store Math.atan(-zeta).
             atanZeta = Math.atan(-zeta);
             scale = Math.pow(1 + zeta * zeta, 0.5 / meps1);
             // Note: These terms are used interchangeably in formulas
             //    1         1
             // -------  = -----
             // (1-eps)    alpha
             inv1mEps = 1.0 / meps1;
             //    1             eps     (1-alpha)     1
             // -------  - 1 = ------- = --------- = ----- - 1
             // (1-eps)        (1-eps)     alpha     alpha
             epsDiv1mEps = inv1mEps - 1;
 
 
             // Compute the support. This applies when alpha < 1 and beta = +/-1
             if (alpha < 1 && Math.abs(beta) == 1) {
                 if (beta == 1) {
                     // alpha < 0, beta = 1
                     lower = zeta;
                     upper = UPPER;
                 } else {
                     // alpha < 0, beta = -1
                     lower = LOWER;
                     upper = zeta;
                 }
             } else {
                 lower = LOWER;
                 upper = UPPER;
             }
         }
 
         /**
          * @param rng Underlying source of randomness
          * @param source Source to copy.
          */
         WeronStableSampler(UniformRandomProvider rng, WeronStableSampler source) {
             super(rng);
             this.eps = source.eps;
             this.meps1 = source.meps1;
             this.zeta = source.zeta;
             this.atanZeta = source.atanZeta;
             this.scale = source.scale;
             this.inv1mEps = source.inv1mEps;
             this.epsDiv1mEps = source.epsDiv1mEps;
             this.lower = source.lower;
             this.upper = source.upper;
         }
 
         @Override
         public double sample() {
             final double phi = getPhi();
             final double w = getOmega();
             return createSample(phi, w);
         }
 
         /**
          * Create the sample. This routine is robust to edge cases and returns a deviate
          * at the extremes of the support. It correctly handles {@code alpha -> 0} when
          * the sample is increasingly likely to be +/- infinity.
          *
          * @param phi Uniform deviate in {@code (-pi/2, pi/2)}
          * @param w Exponential deviate
          * @return x
          */
         protected double createSample(double phi, double w) {
             // Here we use the formula provided by Weron for the 1-parameterization.
             // 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:
             // xi = atan(-zeta) / alpha
             // 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)
 
             // 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), inv1mEps);
 
             // 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 phi=(-pi/2, pi/2) and direct
             // use of arctan(-zeta).
             final double t2 = Math.pow(Math.cos(eps * phi - atanZeta) / w, epsDiv1mEps);
             if (t2 == 0) {
                 return zeta;
             }
 
             final double x = t1 * t2 * scale + zeta;
 
             // Check the bounds. Applies when alpha < 1 and beta = +/-1.
             if (x <= lower) {
                 return lower;
             }
             return x < upper ? x : upper;
         }
 
         @Override
         public WeronStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
             return new WeronStableSampler(rng, this);
         }
     }
 
     /**
      * Implement the generic stable distribution case: {@code alpha < 2} and
      * {@code beta != 0}. This routine assumes {@code alpha != 1}.
      *
      * <p>Implements the Chambers-Mallows-Stuck (CMS) method from Chambers, et al
      * (1976) A Method for Simulating Stable Random Variables. Journal of the
      * American Statistical Association Vol. 71, No. 354, pp. 340-344.
      *
      * <p>The formula produces a stable deviate from the 0-parameterization that is
      * continuous at {@code alpha=1}.
      *
      * <p>This is an implementation of the Fortran routine RSTAB. In the event the
      * computation fails then an alternative computation is performed using the
      * formula provided in Weron (1996) "On the Chambers-Mallows-Stuck method for
      * simulating skewed stable random variables" Statistics &amp; Probability
      * Letters. 28 (2): 165–171. This method is easier to correct from infinite and
      * NaN results. The error correction path is extremely unlikely to occur during
      * use unless {@code alpha -> 0}. In general use it requires the random deviates
      * w or u are extreme. See the unit tests for conditions that create them.
      *
      * <p>This produces non-NaN output for all parameters alpha, beta, u and w with
      * the correct orientation for extremes of the distribution support.
      * The formulas used are symmetric with regard to beta and u.
      */
     static class CMSStableSampler extends WeronStableSampler {
         /** 1/2. */
         private static final double HALF = 0.5;
         /** Cache of expression value used in generation. */
         private final double tau;
 
         /**
          * @param rng Underlying source of randomness
          * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
          * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
          */
         CMSStableSampler(UniformRandomProvider rng, double alpha, double beta) {
             super(rng, alpha, beta);
 
             // Compute the RSTAB pre-factor.
             tau = getTau(alpha, beta);
         }
 
         /**
          * @param rng Underlying source of randomness
          * @param source Source to copy.
          */
         CMSStableSampler(UniformRandomProvider rng, CMSStableSampler source) {
             super(rng, source);
             this.tau = source.tau;
         }
 
         /**
          * Gets tau. This is a factor used in the CMS algorithm. If this is zero then
          * a special case of {@code beta -> 0} has occurred.
          *
          * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
          * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
          * @return tau
          */
         static double getTau(double alpha, double beta) {
             final double eps = 1 - alpha;
             final double meps1 = 1 - eps;
             // Compute RSTAB prefactor
             final double tau;
 
             // tau is symmetric around alpha=1
             // tau -> beta / pi/2 as alpha -> 1
             // tau -> 0 as alpha -> 2 or 0
             // Avoid calling tan as the value approaches the domain limit [-pi/2, pi/2].
             if (Math.abs(eps) < HALF) {
                 // 0.5 < alpha < 1.5. Note: This works when eps=0 as tan(0) / 0 == 1.
                 tau = beta / (SpecialMath.tan2(eps * PI_2) * PI_2);
             } else {
                 // alpha >= 1.5 or alpha <= 0.5.
                 // Do not call tan with alpha > 1 as it wraps in the domain [-pi/2, pi/2].
                 // Since pi is approximate the symmetry is lost by wrapping.
                 // Keep within the domain using (2-alpha).
                 if (meps1 > 1) {
                     tau = beta * PI_2 * eps * (2 - meps1) * -SpecialMath.tan2((2 - meps1) * PI_2);
                 } else {
                     tau = beta * PI_2 * eps * meps1 * SpecialMath.tan2(meps1 * PI_2);
                 }
             }
 
             return tau;
         }
 
         @Override
         public double sample() {
             final double phiby2 = getPhiBy2();
             final double w = getOmega();
 
             // Compute as per the RSTAB routine.
 
             // 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
             // The method is implemented as per the RSTAB routine with the exceptions:
             // 1. The function tan2(x) is implemented with a higher precision approximation.
             // 2. The sample is tested against the expected distribution support.
             // Infinite intermediate terms that create infinite or NaN are corrected by
             // switching the formula and handling infinite terms.
 
             // Compute some tangents
             // a in (-1, 1)
             // bb in [1, 4/pi)
             // b in (-1, 1)
             final double a = phiby2 * SpecialMath.tan2(phiby2);
             final double bb = SpecialMath.tan2(eps * phiby2);
             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.
             // numerator=0 is not possible *in theory* when the uniform deviate generating phi
             // is in the open interval (0, 1). In practice it is possible to obtain <=0 due
             // to round-off error, typically when beta -> +/-1 and phiby2 -> -/+pi/4.
             // This can happen for any alpha.
             final double z = a2p * (b2 + 2 * phiby2 * bb * tau) / (w * a2 * b2p);
 
             // Compute the exponential-type expression
             // Note: z may be infinite, typically when w->0 and a2->0.
             // This can produce NaN under certain parameterizations due to multiplication by 0.
             final double alogz = Math.log(z);
             final double d = SpecialMath.d2(epsDiv1mEps * alogz) * (alogz * inv1mEps);
 
             // Pre-compute the multiplication factor.
             // The numerator may be zero. The denominator is not zero as a2 is bounded to
             // above zero when the uniform deviate that generates phiby2 is not 0 or 1.
             // The min value of a2 is 2^-52. Assume f cannot be infinite as the numerator
             // is computed with a in (-1, 1); b in (-1, 1); phiby2 in (-pi/4, pi/4); tau in
             // [-2/pi, 2/pi]; bb in [1, 4/pi); a2 in (0, 1] limiting the numerator magnitude.
             final double f = (2 * ((a - b) * (1 + a * b) - phiby2 * tau * bb * (b * a2 - 2 * a))) /
                     (a2 * b2p);
 
             // Compute the stable deviate:
             final double x = (1 + eps * d) * f + tau * d;
 
             // Test the support
             if (lower < x && x < upper) {
                 return x;
             }
 
             // Error correction path:
             // x is at the bounds, infinite or NaN (created by 0 / 0,  0 * inf, or inf - inf).
             // This is caused by extreme parameterizations of alpha or beta, or extreme values
             // from the random deviates.
             // Alternatively alpha < 1 and beta = +/-1 and the sample x is at the edge or
             // outside the support due to floating point error.
 
             // Here we use the formula provided by Weron which is easier to correct
             // when deviates are extreme or alpha -> 0. The formula is not continuous
             // as alpha -> 1 without a shift which reduces the precision of the sample;
             // for rare edge case correction this has minimal effect on sampler output.
             return createSample(phiby2 * 2, w);
         }
 
         @Override
         public CMSStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
             return new CMSStableSampler(rng, this);
         }
     }
 
     /**
      * Implement the stable distribution case: {@code alpha == 1} and {@code beta != 0}.
      *
      * <p>Implements the same algorithm as the {@link CMSStableSampler} with
      * the {@code alpha} assumed to be 1.
      *
      * <p>This sampler specifically requires that {@code beta / (pi/2) != 0}; otherwise
      * the parameters equal {@code alpha=1, beta=0} as the Cauchy distribution case.
      */
     static class Alpha1CMSStableSampler extends BaseStableSampler {
         /** Cache of expression value used in generation. */
         private final double tau;
 
         /**
          * @param rng Underlying source of randomness
          * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
          */
         Alpha1CMSStableSampler(UniformRandomProvider rng, double beta) {
             super(rng);
             tau = beta / PI_2;
         }
 
         /**
          * @param rng Underlying source of randomness
          * @param source Source to copy.
          */
         Alpha1CMSStableSampler(UniformRandomProvider rng, Alpha1CMSStableSampler source) {
             super(rng);
             this.tau = source.tau;
         }
 
         @Override
         public double sample() {
             final double phiby2 = getPhiBy2();
             final double w = getOmega();
 
             // Compute some tangents
             final double a = phiby2 * SpecialMath.tan2(phiby2);
 
             // Compute some necessary subexpressions
             final double da = a * a;
             final double a2 = 1 - da;
             final double a2p = 1 + da;
 
             // Compute coefficient.
 
             // numerator=0 is not possible when the uniform deviate generating phi
             // is in the open interval (0, 1) and alpha=1.
             final double z = a2p * (1 + 2 * phiby2 * tau) / (w * a2);
 
             // Compute the exponential-type expression
             // Note: z may be infinite, typically when w->0 and a2->0.
             // This can produce NaN under certain parameterizations due to multiplication by 0.
             // When alpha=1 the expression
             // d = d2((eps / (1-eps)) * alogz) * (alogz / (1-eps)) is eliminated to 1 * log(z)
             final double d = Math.log(z);
 
             // Pre-compute the multiplication factor.
             final double f = (2 * (a - phiby2 * tau * (-2 * a))) / a2;
 
             // Compute the stable deviate:
             // This does not require correction as f is finite (as per the alpha != 1 case),
             // tau is non-zero and only d can be infinite due to an extreme w -> 0.
             return f + tau * d;
         }
 
         @Override
         public Alpha1CMSStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
             return new Alpha1CMSStableSampler(rng, this);
         }
     }
 
     /**
      * Implement the generic stable distribution case: {@code alpha < 2} and {@code beta == 0}.
      *
      * <p>Implements the same algorithm as the {@link WeronStableSampler} with
      * the {@code beta} assumed to be 0.
      *
      * <p>This routine assumes {@code alpha != 1}; {@code alpha=1, beta=0} is the Cauchy
      * distribution case.
      */
     static class Beta0WeronStableSampler extends BaseStableSampler {
         /** Epsilon (1 - alpha). */
         protected final double eps;
         /** Epsilon (1 - alpha). */
         protected final double meps1;
         /** 1 / alpha = 1 / (1 - eps). */
         protected final double inv1mEps;
         /** (1 / alpha) - 1 = eps / (1 - eps). */
         protected final double epsDiv1mEps;
 
         /**
          * @param rng Underlying source of randomness
          * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
          */
         Beta0WeronStableSampler(UniformRandomProvider rng, double alpha) {
             super(rng);
             eps = 1 - alpha;
             meps1 = 1 - eps;
             inv1mEps = 1.0 / meps1;
             epsDiv1mEps = inv1mEps - 1;
         }
 
         /**
          * @param rng Underlying source of randomness
          * @param source Source to copy.
          */
         Beta0WeronStableSampler(UniformRandomProvider rng, Beta0WeronStableSampler source) {
             super(rng);
             this.eps = source.eps;
             this.meps1 = source.meps1;
             this.inv1mEps = source.inv1mEps;
             this.epsDiv1mEps = source.epsDiv1mEps;
         }
 
         @Override
         public double sample() {
             final double phi = getPhi();
             final double w = getOmega();
             return createSample(phi, w);
         }
 
         /**
          * Create the sample. This routine is robust to edge cases and returns a deviate
          * at the extremes of the support. It correctly handles {@code alpha -> 0} when
          * the sample is increasingly likely to be +/- infinity.
          *
          * @param phi Uniform deviate in {@code (-pi/2, pi/2)}
          * @param w Exponential deviate
          * @return x
          */
         protected double createSample(double phi, double w) {
             // As per the Weron sampler with beta=0 and terms eliminated.
             // Note that if alpha=1 this reduces to sin(phi) / cos(phi) => Cauchy case.
 
             // Compute terms.
             // Either term can be infinite or 0. Certain parameters compute 0 * inf.
             // Detect zeros and return as 0.
 
             // Note sin(alpha * phi) is only ever zero when phi=0. No value of alpha
             // multiplied by small phi can create zero due to the limited
             // precision of alpha imposed by alpha = 1 - (1-alpha). At this point cos(phi) = 1.
             // Thus 0/0 cannot occur.
             final double t1 = Math.sin(meps1 * phi) / Math.pow(Math.cos(phi), inv1mEps);
             if (t1 == 0) {
                 return 0;
             }
             final double t2 = Math.pow(Math.cos(eps * phi) / w, epsDiv1mEps);
             if (t2 == 0) {
                 return 0;
             }
             return t1 * t2;
         }
 
         @Override
         public Beta0WeronStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
             return new Beta0WeronStableSampler(rng, this);
         }
     }
 
     /**
      * Implement the generic stable distribution case: {@code alpha < 2} and {@code beta == 0}.
      *
      * <p>Implements the same algorithm as the {@link CMSStableSampler} with
      * the {@code beta} assumed to be 0.
      *
      * <p>This routine assumes {@code alpha != 1}; {@code alpha=1, beta=0} is the Cauchy
      * distribution case.
      */
     static class Beta0CMSStableSampler extends Beta0WeronStableSampler {
         /**
          * @param rng Underlying source of randomness
          * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
          */
         Beta0CMSStableSampler(UniformRandomProvider rng, double alpha) {
             super(rng, alpha);
         }
 
         /**
          * @param rng Underlying source of randomness
          * @param source Source to copy.
          */
         Beta0CMSStableSampler(UniformRandomProvider rng, Beta0CMSStableSampler source) {
             super(rng, source);
         }
 
         @Override
         public double sample() {
             final double phiby2 = getPhiBy2();
             final double w = getOmega();
 
             // Compute some tangents
             final double a = phiby2 * SpecialMath.tan2(phiby2);
             final double b = eps * phiby2 * SpecialMath.tan2(eps * phiby2);
             // Compute some necessary subexpressions
             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;
             // Compute coefficient.
             final double z = a2p * b2 / (w * a2 * b2p);
 
             // Compute the exponential-type expression
             final double alogz = Math.log(z);
             final double d = SpecialMath.d2(epsDiv1mEps * alogz) * (alogz * inv1mEps);
 
             // Pre-compute the multiplication factor.
             // The numerator may be zero. The denominator is not zero as a2 is bounded to
             // above zero when the uniform deviate that generates phiby2 is not 0 or 1.
             final double f = (2 * ((a - b) * (1 + a * b))) / (a2 * b2p);
 
             // Compute the stable deviate:
             final double x = (1 + eps * d) * f;
 
             // Test the support
             if (LOWER < x && x < UPPER) {
                 return x;
             }
 
             // Error correction path.
             // Here we use the formula provided by Weron which is easier to correct
             // when deviates are extreme or alpha -> 0.
             return createSample(phiby2 * 2, w);
         }
 
         @Override
         public Beta0CMSStableSampler withUniformRandomProvider(UniformRandomProvider rng) {
             return new Beta0CMSStableSampler(rng, this);
         }
     }
 
     /**
      * Implement special math functions required by the CMS algorithm.
      */
     static final class SpecialMath {
         /** pi/4. */
         private static final double PI_4 = Math.PI / 4;
         /** 4/pi. */
         private static final double FOUR_PI = 4 / Math.PI;
         /** tan2 product constant. */
         private static final double P0 = -0.5712939549476836914932149599e10;
         /** tan2 product constant. */
         private static final double P1 = 0.4946855977542506692946040594e9;
         /** tan2 product constant. */
         private static final double P2 = -0.9429037070546336747758930844e7;
         /** tan2 product constant. */
         private static final double P3 = 0.5282725819868891894772108334e5;
         /** tan2 product constant. */
         private static final double P4 = -0.6983913274721550913090621370e2;
         /** tan2 quotient constant. */
         private static final double Q0 = -0.7273940551075393257142652672e10;
         /** tan2 quotient constant. */
         private static final double Q1 = 0.2125497341858248436051062591e10;
         /** tan2 quotient constant. */
         private static final double Q2 = -0.8000791217568674135274814656e8;
         /** tan2 quotient constant. */
         private static final double Q3 = 0.8232855955751828560307269007e6;
         /** tan2 quotient constant. */
         private static final double Q4 = -0.2396576810261093558391373322e4;
         /**
          * The threshold to switch to using {@link Math#expm1(double)}. The following
          * table shows the mean (max) ULP difference between using expm1 and exp using
          * random +/-x with different exponents (n=2^30):
          *
          * <pre>
          * x        exp  positive x                 negative x
          * 64.0      6   0.10004021506756544  (2)   0.0                   (0)
          * 32.0      5   0.11177831795066595  (2)   0.0                   (0)
          * 16.0      4   0.0986650362610817   (2)   9.313225746154785E-10 (1)
          * 8.0       3   0.09863092936575413  (2)   4.9658119678497314E-6 (1)
          * 4.0       2   0.10015273280441761  (2)   4.547201097011566E-4  (1)
          * 2.0       1   0.14359260816127062  (2)   0.005623611621558666  (2)
          * 1.0       0   0.20160607434809208  (2)   0.03312791418284178   (2)
          * 0.5      -1   0.3993037799373269   (2)   0.28186883218586445   (2)
          * 0.25     -2   0.6307008266448975   (2)   0.5192863345146179    (2)
          * 0.125    -3   1.3862918205559254   (4)   1.386285437270999     (4)
          * 0.0625   -4   2.772640804760158    (8)   2.772612397558987     (8)
          * </pre>
          *
          * <p>The threshold of 0.5 has a mean ULP below 0.5 and max ULP of 2. The
          * transition is monotonic. Neither is true for the next threshold of 0.25.
          */
         private static final double SWITCH_TO_EXPM1 = 0.5;
 
         /** No instances. */
         private SpecialMath() {}
 
         /**
          * Evaluate {@code (exp(x) - 1) / x}. 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}.
          * <li>For {@code x=nan} this returns {@code nan}.
          * </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}.
          */
         static double d2(double x) {
             // Here expm1 is only used when use of expm1 and exp consistently
             // compute different results by more than 0.5 ULP.
             if (Math.abs(x) < SWITCH_TO_EXPM1) {
                 // Deliberate comparison to floating-point zero
                 if (x == 0) {
                     // Avoid 0 / 0 error
                     return 1.0;
                 }
                 return Math.expm1(x) / x;
             }
             // No use of expm1. Accuracy as x moves away from 0 is not required as the result
             // is divided by x and the accuracy of the final result is within a few ULP.
             if (x < Double.POSITIVE_INFINITY) {
                 return (Math.exp(x) - 1) / x;
             }
             // Upper bound (or NaN)
             return x;
         }
 
         /**
          * Evaluate {@code tan(x) / x}.
          *
          * <p>For {@code x} in the range {@code [0, pi/4]} this returns
          * a value in the range {@code [1, 4 / pi]}.
          *
          * <p>The following properties are desirable for the CMS algorithm:
          *
          * <ul>
          * <li>For {@code x=0} this returns {@code 1}.
          * <li>For {@code x=pi/4} this returns {@code 4/pi}.
          * <li>For {@code x=pi/4} this multiplied by {@code x} returns {@code 1}.
          * </ul>
          *
          * <p>This method is called by the CMS algorithm when {@code x < pi/4}.
          * In this case the method is almost as accurate as {@code Math.tan(x) / x}, does
          * not require checking for {@code x=0} and is faster.
          *
          * @param x the x
          * @return {@code tan(x) / x}.
          */
         static double tan2(double x) {
             if (Math.abs(x) > PI_4) {
                 // Reduction is not supported. Delegate to the JDK.
                 return Math.tan(x) / x;
             }
 
             // Testing with approximation 4283 from Hart et al, as used in the RSTAB
             // routine, showed the method was not accurate enough for use with
             // double computation. Hart et al state it has max relative error = 1e-10.66.
             // For tan(x) / x with x in [0, pi/4] values outside [1, 4/pi] were computed.
             // When testing verses Math.tan(x) the mean ULP difference is 93436.3.
 
             // Approximation 4288 from Hart et al (1968, P. 252).
             // Max relative error = 1e-26.68 (for tan(x)).
             // When testing verses Math.tan(x) the mean ULP difference is 0.590597.
 
             // The approximation is defined as:
             // tan(x*pi/4) = x * P(x^2) / Q(x^2)
             //   with P and Q polynomials of x squared.
             //
             // To create tan(x):
             // tan(x) = xi * P(xi^2) / Q(xi^2), xi = x * 4/pi
             // tan(x) / x = xi * P(xi^2) / Q(xi^2) / x
             // tan(x) / x = 4/pi * (P(xi^2) / Q(xi^2))
             //            = P(xi^2) / (pi/4 * Q(xi^2))
             // The later has a smaller mean ULP difference to Math.tan(x) / x.
             final double xi = x * FOUR_PI;
 
             // Use the power form with a reverse summation order to have smaller
             // magnitude terms first. Note: x < 1 so greater powers are smaller.
             // This has essentially the same accuracy as the nested form of the polynomials
             // for a marginal performance increase. See JMH examples for performance tests.
             final double x2 = xi * xi;
             final double x4 = x2 * x2;
             final double x6 = x4 * x2;
             final double x8 = x4 * x4;
             return          (x8 * P4 + x6 * P3 + x4 * P2 + x2 * P1 + P0) /
                     (PI_4 * (x8 * x2 + x8 * Q4 + x6 * Q3 + x4 * Q2 + x2 * Q1 + Q0));
         }
     }
 
     /**
      * @param rng Generator of uniformly distributed random numbers.
      */
     StableSampler(UniformRandomProvider rng) {
         this.rng = rng;
     }
 
     /**
      * Generate a sample from a stable distribution.
      *
      * <p>The distribution uses the 0-parameterization: S(alpha, beta, gamma, delta; 0).
      */
     @Override
     public abstract double sample();
 
     /** {@inheritDoc} */
     // Redeclare the signature to return a StableSampler not a SharedStateContinuousSampler
     @Override
     public abstract StableSampler withUniformRandomProvider(UniformRandomProvider rng);
 
     /**
      * Generates a {@code long} value.
      * Used by algorithm implementations without exposing access to the RNG.
      *
      * @return the next random value
      */
     long nextLong() {
         return rng.nextLong();
     }
 
     /** {@inheritDoc} */
     @Override
     public String toString() {
         // All variations use the same string representation, i.e. no changes
         // for the Gaussian, Levy or Cauchy case.
         return "Stable deviate [" + rng.toString() + "]";
     }
 
     /**
      * Creates a standardized sampler of a stable distribution with zero location and unit scale.
      *
      * <p>Special cases:
      *
      * <ul>
      * <li>{@code alpha=2} returns a Gaussian distribution sampler with
      *     {@code mean=0} and {@code variance=2} (Note: {@code beta} has no effect on the distribution).
      * <li>{@code alpha=1} and {@code beta=0} returns a Cauchy distribution sampler with
      *     {@code location=0} and {@code scale=1}.
      * <li>{@code alpha=0.5} and {@code beta=1} returns a Levy distribution sampler with
      *     {@code location=-1} and {@code scale=1}. This location shift is due to the
      *     0-parameterization of the stable distribution.
      * </ul>
      *
      * <p>Note: To allow the computation of the stable distribution the parameter alpha
      * is validated using {@code 1 - alpha} in the interval {@code [-1, 1)}.
      *
      * @param rng Generator of uniformly distributed random numbers.
      * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
      * @return the sampler
      * @throws IllegalArgumentException if {@code 1 - alpha < -1}; or {@code 1 - alpha >= 1};
      * or {@code beta < -1}; or {@code beta > 1}.
      */
     public static StableSampler of(UniformRandomProvider rng,
                                    double alpha,
                                    double beta) {
         validateParameters(alpha, beta);
         return create(rng, alpha, beta);
     }
 
     /**
      * Creates a sampler of a stable distribution. This applies a transformation to the
      * standardized sampler.
      *
      * <p>The random variable \( X \) has
      * the stable distribution \( S(\alpha, \beta, \gamma, \delta; 0) \) if:
      *
      * <p>\[ X = \gamma Z_0 + \delta \]
      *
      * <p>where \( Z_0 = S(\alpha, \beta; 0) \) is a standardized stable distribution.
      *
      * <p>Note: To allow the computation of the stable distribution the parameter alpha
      * is validated using {@code 1 - alpha} in the interval {@code [-1, 1)}.
      *
      * @param rng Generator of uniformly distributed random numbers.
      * @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 gamma Scale parameter. Must be strictly positive and finite.
      * @param delta Location parameter. Must be finite.
      * @return the sampler
      * @throws IllegalArgumentException if {@code 1 - alpha < -1}; or {@code 1 - alpha >= 1};
      * or {@code beta < -1}; or {@code beta > 1}; or {@code gamma <= 0}; or
      * {@code gamma} or {@code delta} are not finite.
      * @see #of(UniformRandomProvider, double, double)
      */
     public static StableSampler of(UniformRandomProvider rng,
                                    double alpha,
                                    double beta,
                                    double gamma,
                                    double delta) {
         validateParameters(alpha, beta, gamma, delta);
 
         // Choose the algorithm.
         // Reuse the special cases as they have transformation support.
 
         if (alpha == ALPHA_GAUSSIAN) {
             // Note: beta has no effect and is ignored.
             return new GaussianStableSampler(rng, gamma, delta);
         }
 
         // Note: As beta -> 0 the result cannot be computed differently to beta = 0.
         if (alpha == ALPHA_CAUCHY && CMSStableSampler.getTau(ALPHA_CAUCHY, beta) == TAU_ZERO) {
             return new CauchyStableSampler(rng, gamma, delta);
         }
 
         if (alpha == ALPHA_LEVY && Math.abs(beta) == BETA_LEVY) {
             // Support mirroring for negative beta by inverting the beta=1 Levy sample
             // using a negative gamma. Note: The delta is not mirrored as it is a shift
             // applied to the scaled and mirrored distribution.
             return new LevyStableSampler(rng, beta * gamma, delta);
         }
 
         // Standardized sampler
         final StableSampler sampler = create(rng, alpha, beta);
         // Transform
         return new TransformedStableSampler(sampler, gamma, delta);
     }
 
     /**
      * Creates a standardized sampler of a stable distribution with zero location and unit scale.
      *
      * @param rng Generator of uniformly distributed random numbers.
      * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
      * @return the sampler
      */
     private static StableSampler create(UniformRandomProvider rng,
                                         double alpha,
                                         double beta) {
         // Choose the algorithm.
         // The special case samplers have transformation support and use gamma=1.0, delta=0.0.
         // As alpha -> 0 the computation increasingly requires correction
         // of infinity to the distribution support.
 
         if (alpha == ALPHA_GAUSSIAN) {
             // Note: beta has no effect and is ignored.
             return new GaussianStableSampler(rng, GAMMA_1, DELTA_0);
         }
 
         // Note: As beta -> 0 the result cannot be computed differently to beta = 0.
         // This is based on the computation factor tau:
         final double tau = CMSStableSampler.getTau(alpha, beta);
 
         if (tau == TAU_ZERO) {
             // Symmetric case (beta skew parameter is effectively zero)
             if (alpha == ALPHA_CAUCHY) {
                 return new CauchyStableSampler(rng, GAMMA_1, DELTA_0);
             }
             if (alpha <= ALPHA_SMALL) {
                 // alpha -> 0 requires robust error correction
                 return new Beta0WeronStableSampler(rng, alpha);
             }
             return new Beta0CMSStableSampler(rng, alpha);
         }
 
         // Here beta is significant.
 
         if (alpha == 1) {
             return new Alpha1CMSStableSampler(rng, beta);
         }
 
         if (alpha == ALPHA_LEVY && Math.abs(beta) == BETA_LEVY) {
             // Support mirroring for negative beta by inverting the beta=1 Levy sample
             // using a negative gamma. Note: The delta is not mirrored as it is a shift
             // applied to the scaled and mirrored distribution.
             return new LevyStableSampler(rng, beta, DELTA_0);
         }
 
         if (alpha <= ALPHA_SMALL) {
             // alpha -> 0 requires robust error correction
             return new WeronStableSampler(rng, alpha, beta);
         }
 
         return new CMSStableSampler(rng, alpha, beta);
     }
 
     /**
      * Validate the parameters are in the correct range.
      *
      * @param alpha Stability parameter. Must be in the interval {@code (0, 2]}.
      * @param beta Skewness parameter. Must be in the interval {@code [-1, 1]}.
      * @throws IllegalArgumentException if {@code 1 - alpha < -1}; or {@code 1 - alpha >= 1};
      * or {@code beta < -1}; or {@code beta > 1}.
      */
     private static void validateParameters(double alpha, double beta) {
         // The epsilon (1-alpha) value must be in the interval [-1, 1).
         // Logic inversion will identify NaN
         final double eps = 1 - alpha;
         if (!(-1 <= eps && eps < 1)) {
             throw new IllegalArgumentException("alpha is not in the interval (0, 2]: " + alpha);
         }
         if (!(-1 <= beta && beta <= 1)) {
             throw new IllegalArgumentException("beta is not in the interval [-1, 1]: " + beta);
         }
     }
 
     /**
      * Validate the parameters are in the correct range.
      *
      * @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 gamma Scale parameter. Must be strictly positive and finite.
      * @param delta Location parameter. Must be finite.
      * @throws IllegalArgumentException if {@code 1 - alpha < -1}; or {@code 1 - alpha >= 1};
      * or {@code beta < -1}; or {@code beta > 1}; or {@code gamma <= 0}; or
      * {@code gamma} or {@code delta} are not finite.
      */
     private static void validateParameters(double alpha, double beta,
                                            double gamma, double delta) {
         validateParameters(alpha, beta);
 
         InternalUtils.requireStrictlyPositiveFinite(gamma, "gamma");
         InternalUtils.requireFinite(delta, "delta");
     }
 }