public double sample()

in commons-rng-sampling/src/main/java/org/apache/commons/rng/sampling/distribution/StableSampler.java [785:865]
23 lines of code
3 McCabe index (conditional complexity)

        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);
        }