Class BoostErf


  • final class BoostErf
    extends java.lang.Object
    Implementation of the error function and its inverse.

    This code has been adapted from the Boost c++ implementation <boost/math/special_functions/erf.hpp>. The erf/erfc functions and their inverses are copyright John Maddock 2006 and subject to the Boost Software License.

    Additions made to support the erfcx function are original work under the Apache software license.

    See Also:
    Boost C++ Error functions
    • Field Summary

      Fields 
      Modifier and Type Field Description
      private static double COMPUTE_ERF
      Threshold for the erf implementation for |x| where the computation uses erf(x); otherwise erfc(x) is computed.
      private static double ERFCX_APPROX
      Threshold for the scaled complementary error function erfcx where the approximation (1 / sqrt(pi)) / x can be used.
      private static double ERFCX_NEG_X_MAX
      Threshold for the scaled complementary error function erfcx for negative x where 2 * exp(x*x) will overflow.
      private static double EXP_XX_1
      Threshold for the scaled complementary error function erfcx for x where exp(x*x) == 1; x <= t.
      private static double MULTIPLIER
      The multiplier used to split the double value into high and low parts.
      private static double ONE_OVER_ROOT_PI
      1 / sqrt(pi).
    • Constructor Summary

      Constructors 
      Modifier Constructor Description
      private BoostErf()
      Private constructor.
    • Method Summary

      All Methods Static Methods Concrete Methods 
      Modifier and Type Method Description
      (package private) static double erf​(double x)
      Returns the error function.
      (package private) static double erfc​(double x)
      Returns the complementary error function.
      (package private) static double erfcInv​(double z)
      Returns the inverse complementary error function.
      (package private) static double erfcx​(double x)
      Returns the scaled complementary error function.
      private static double erfImp​(double z, boolean invert, boolean scaled)
      53-bit implementation for the error function.
      (package private) static double erfInv​(double z)
      Returns the inverse error function.
      private static double erfInvImp​(double p, double q)
      Common implementation for inverse erf and erfc functions.
      (package private) static double expmxx​(double x)
      Compute exp(-x*x) with high accuracy.
      (package private) static double expxx​(double x)
      Compute exp(x*x) with high accuracy.
      private static double expxx​(double a, double b)
      Compute exp(a+b) with high accuracy assuming a+b = a.
      private static double highPartUnscaled​(double value)
      Implement Dekker's method to split a value into two parts.
      private static double squareLow​(double hx, double lx, double xx)
      Compute the low part of the double length number (z,zz) for the exact square of x using Dekker's mult12 algorithm.
      private static double squareLowUnscaled​(double x, double xx)
      Compute the low part of the double length number (z,zz) for the exact square of x using Dekker's mult12 algorithm.
      • Methods inherited from class java.lang.Object

        clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
    • Field Detail

      • MULTIPLIER

        private static final double MULTIPLIER
        The multiplier used to split the double value into high and low parts. From Dekker (1971): "The constant should be chosen equal to 2^(p - p/2) + 1, where p is the number of binary digits in the mantissa". Here p is 53 and the multiplier is 2^27 + 1.
        See Also:
        Constant Field Values
      • ONE_OVER_ROOT_PI

        private static final double ONE_OVER_ROOT_PI
        1 / sqrt(pi). Used for the scaled complementary error function erfcx.
        See Also:
        Constant Field Values
      • ERFCX_APPROX

        private static final double ERFCX_APPROX
        Threshold for the scaled complementary error function erfcx where the approximation (1 / sqrt(pi)) / x can be used.
        See Also:
        Constant Field Values
      • COMPUTE_ERF

        private static final double COMPUTE_ERF
        Threshold for the erf implementation for |x| where the computation uses erf(x); otherwise erfc(x) is computed. The final result is achieved by suitable application of symmetry.
        See Also:
        Constant Field Values
      • ERFCX_NEG_X_MAX

        private static final double ERFCX_NEG_X_MAX
        Threshold for the scaled complementary error function erfcx for negative x where 2 * exp(x*x) will overflow. Value is 26.62873571375149.
      • EXP_XX_1

        private static final double EXP_XX_1
        Threshold for the scaled complementary error function erfcx for x where exp(x*x) == 1; x <= t. Value is (1 + 5/16) * 2^-27 = 9.778887033462524E-9.

        Note: This is used for performance. If set to 0 then the result is computed using expm1(x*x) with the same final result.

        See Also:
        Constant Field Values
    • Constructor Detail

      • BoostErf

        private BoostErf()
        Private constructor.
    • Method Detail

      • erfc

        static double erfc​(double x)
        Returns the complementary error function.
        Parameters:
        x - the value.
        Returns:
        the complementary error function.
      • erf

        static double erf​(double x)
        Returns the error function.
        Parameters:
        x - the value.
        Returns:
        the error function.
      • erfImp

        private static double erfImp​(double z,
                                     boolean invert,
                                     boolean scaled)
        53-bit implementation for the error function.

        Note: The scaled flag only applies when z >= 0.5 and invert == true. This functionality is used to compute erfcx(z) for positive z.

        Parameters:
        z - Point to evaluate
        invert - true to invert the result (for the complementary error function)
        scaled - true to compute the scaled complementary error function
        Returns:
        the error function result
      • erfcx

        static double erfcx​(double x)
        Returns the scaled complementary error function.
         erfcx(x) = exp(x^2) * erfc(x)
         
        Parameters:
        x - the value.
        Returns:
        the scaled complementary error function.
      • erfcInv

        static double erfcInv​(double z)
        Returns the inverse complementary error function.
        Parameters:
        z - Value (in [0, 2]).
        Returns:
        t such that z = erfc(t)
      • erfInv

        static double erfInv​(double z)
        Returns the inverse error function.
        Parameters:
        z - Value (in [-1, 1]).
        Returns:
        t such that z = erf(t)
      • erfInvImp

        private static double erfInvImp​(double p,
                                        double q)
        Common implementation for inverse erf and erfc functions.
        Parameters:
        p - P-value
        q - Q-value (1-p)
        Returns:
        the inverse
      • expxx

        static double expxx​(double x)
        Compute exp(x*x) with high accuracy. This is performed using information in the round-off from x*x.

        This is accurate at large x to 1 ulp.

        At small x the accuracy cannot be improved over using exp(x*x). This occurs at x <= 1.

        Warning: This has no checks for overflow. The method is never called when x*x > log(MAX_VALUE/2).

        Parameters:
        x - Value
        Returns:
        exp(x*x)
      • expmxx

        static double expmxx​(double x)
        Compute exp(-x*x) with high accuracy. This is performed using information in the round-off from x*x.

        This is accurate at large x to 1 ulp until exp(-x*x) is close to sub-normal. For very small exp(-x*x) the adjustment is sub-normal and bits can be lost in the adjustment for a max observed error of < 2 ulp.

        At small x the accuracy cannot be improved over using exp(-x*x). This occurs at x <= 1.

        Parameters:
        x - Value
        Returns:
        exp(-x*x)
      • expxx

        private static double expxx​(double a,
                                    double b)
        Compute exp(a+b) with high accuracy assuming a+b = a.

        This is accurate at large positive a to 1 ulp. If a is negative and exp(a) is close to sub-normal a bit of precision may be lost when adjusting result as the adjustment is sub-normal (max observed error < 2 ulp). For the use case of multiplication of a number less than 1 by exp(-x*x), a = -x*x, the result will be sub-normal and the rounding error is lost.

        At small |a| the accuracy cannot be improved over using exp(a) as the round-off is too small to create terms that can adjust the standard result by more than 0.5 ulp. This occurs at |a| <= 1.

        Parameters:
        a - High bits of a split number
        b - Low bits of a split number
        Returns:
        exp(a+b)
      • squareLowUnscaled

        private static double squareLowUnscaled​(double x,
                                                double xx)
        Compute the low part of the double length number (z,zz) for the exact square of x using Dekker's mult12 algorithm. The standard precision product x*x must be provided. The number x is split into high and low parts using Dekker's algorithm.

        Warning: This method does not perform scaling in Dekker's split and large finite numbers can create NaN results.

        Parameters:
        x - Number to square
        xx - Standard precision product x*x
        Returns:
        the low part of the square double length number
      • highPartUnscaled

        private static double highPartUnscaled​(double value)
        Implement Dekker's method to split a value into two parts. Multiplying by (2^s + 1) creates a big value from which to derive the two split parts.
         c = (2^s + 1) * a
         a_big = c - a
         a_hi = c - a_big
         a_lo = a - a_hi
         a = a_hi + a_lo
         

        The multiplicand allows a p-bit value to be split into (p-s)-bit value a_hi and a non-overlapping (s-1)-bit value a_lo. Combined they have (p-1) bits of significand but the sign bit of a_lo contains a bit of information. The constant is chosen so that s is ceil(p/2) where the precision p for a double is 53-bits (1-bit of the mantissa is assumed to be 1 for a non sub-normal number) and s is 27.

        This conversion does not use scaling and the result of overflow is NaN. Overflow may occur when the exponent of the input value is above 996.

        Splitting a NaN or infinite value will return NaN.

        Parameters:
        value - Value.
        Returns:
        the high part of the value.
        See Also:
        Math.getExponent(double)
      • squareLow

        private static double squareLow​(double hx,
                                        double lx,
                                        double xx)
        Compute the low part of the double length number (z,zz) for the exact square of x using Dekker's mult12 algorithm. The standard precision product x*x must be provided. The number x should already be split into low and high parts.

        Note: This uses the high part of the result (z,zz) as x * x and not hx * hx + hx * lx + lx * hx as specified in Dekker's original paper. See Shewchuk (1997) for working examples.

        Parameters:
        hx - High part of factor.
        lx - Low part of factor.
        xx - Square of the factor.
        Returns:
        lx * ly - (((xy - hx * hy) - lx * hy) - hx * ly)
        See Also:
        Shewchuk (1997) Theorum 18