Probit function in Java

Probit.java

/**
 * Provides a {@link #probit(double)} method to convert from the (0.0,1.0) domain to a Gaussian range.
 * Uses an algorithm by Peter John Acklam, as implemented by Sherali Karimov.
 * Free to use, but you should credit Acklam and Karimov appropriately.
 */
public final class Probit {
    private static final double LOW = 0.02425;
    private static final double HIGH = 1.0 - LOW;

    // Coefficients in rational approximations.
    private static final double
            A0 = -3.969683028665376e+01,
            A1 =  2.209460984245205e+02, 
            A2 = -2.759285104469687e+02,
            A3 =  1.383577518672690e+02, 
            A4 = -3.066479806614716e+01,
            A5 =  2.506628277459239e+00;

    private static final double
            B0 = -5.447609879822406e+01,
            B1 =  1.615858368580409e+02, 
            B2 = -1.556989798598866e+02,
            B3 =  6.680131188771972e+01, 
            B4 = -1.328068155288572e+01;

    private static final double
            C0 = -7.784894002430293e-03, 
            C1 = -3.223964580411365e-01,
            C2 = -2.400758277161838e+00,
            C3 = -2.549732539343734e+00, 
            C4 =  4.374664141464968e+00,
            C5 =  2.938163982698783e+00;

    private static final double
            D0 = 7.784695709041462e-03,
            D1 = 3.224671290700398e-01,
            D2 = 2.445134137142996e+00,
            D3 = 3.754408661907416e+00;

    /**
     * A way of taking a double in the (0.0, 1.0) range and mapping it to a Gaussian or normal distribution, so high
     * inputs correspond to high outputs, and similarly for the low range. This is centered on 0.0 and its standard
     * deviation seems to be 1.0 (the same as {@link Random#nextGaussian()}). It is slightly faster (about a 10%
     * increase in throughput, with comparable overhead for both methods) than the approach used by
     * {@link Random#nextGaussian()}, even when both use a faster RNG algorithm and when both values produced by
     * nextGaussian() are used.
     * <br>
     * This uses an algorithm by Peter John Acklam, as implemented by Sherali Karimov.
     * <a href="https://web.archive.org/web/20150910002142/http://home.online.no/~pjacklam/notes/invnorm/impl/karimov/StatUtil.java">Source</a>.
     * <a href="https://web.archive.org/web/20151030215612/http://home.online.no/~pjacklam/notes/invnorm/">Information on the algorithm</a>.
     * @param d should be between 0 and 1, exclusive, but other values are tolerated (they return infinite results)
     * @return a normal-distributed double centered on 0.0
     */
    @SuppressWarnings("divzero") // This can legitimately return infinite doubles, which it produces with zero division.
    public static double probit(final double d) {
        if (d <= 0 || d >= 1) {
            return (d - 0.5) / 0.0;
        }
        // Rational approximation for lower region:
        else if (d < LOW) {
            final double q = Math.sqrt(-2 * Math.log(d));
            return (((((C0 * q + C1) * q + C2) * q + C3) * q + C4) * q + C5) / ((((D0 * q + D1) * q + D2) * q + D3) * q + 1);
        }
        // Rational approximation for upper region:
        else if (HIGH < d) {
            final double q = Math.sqrt(-2 * Math.log(1 - d));
            return -(((((C0 * q + C1) * q + C2) * q + C3) * q + C4) * q + C5) / ((((D0 * q + D1) * q + D2) * q + D3) * q + 1);
        }
        // Rational approximation for central region:
        else {
            final double q = d - 0.5;
            final double r = q * q;
            return (((((A0 * r + A1) * r + A2) * r + A3) * r + A4) * r + A5) * q / (((((B0 * r + B1) * r + B2) * r + B3) * r + B4) * r + 1);
        }
    }
}