Erf.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      https://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. /*
  19.  * This is not the original file distributed by the Apache Software Foundation
  20.  * It has been modified by the Hipparchus project
  21.  */
  22. package org.hipparchus.special;

  24. import org.hipparchus.CalculusFieldElement;
  25. import org.hipparchus.Field;
  26. import org.hipparchus.util.FastMath;

  28. /**
  29.  * This is a utility class that provides computation methods related to the
  30.  * error functions.
  31.  *
  32.  */
  33. public class Erf {

  35.     /**
  36.      * The number {@code X_CRIT} is used by {@link #erf(double, double)} internally.
  37.      * This number solves {@code erf(x)=0.5} within 1ulp.
  38.      * More precisely, the current implementations of
  39.      * {@link #erf(double)} and {@link #erfc(double)} satisfy:<br>
  40.      * {@code erf(X_CRIT) < 0.5},<br>
  41.      * {@code erf(Math.nextUp(X_CRIT) > 0.5},<br>
  42.      * {@code erfc(X_CRIT) = 0.5}, and<br>
  43.      * {@code erfc(Math.nextUp(X_CRIT) < 0.5}
  44.      */
  45.     private static final double X_CRIT = 0.4769362762044697;

  47.     /**
  48.      * Default constructor.  Prohibit instantiation.
  49.      */
  50.     private Erf() {}

  52.     /**
  53.      * Returns the error function.
  54.      *
  55.      * \[
  56.      * \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{t=0}^x e^{-t^2}dt
  57.      * \]
  58.      *
  59.      * <p>This implementation computes erf(x) using the
  60.      * {@link Gamma#regularizedGammaP(double, double, double, int) regularized gamma function},
  61.      * following <a href="http://mathworld.wolfram.com/Erf.html"> Erf</a>, equation (3)</p>
  62.      *
  63.      * <p>The value returned is always between -1 and 1 (inclusive).
  64.      * If {@code abs(x) > 40}, then {@code erf(x)} is indistinguishable from
  65.      * either 1 or -1 as a double, so the appropriate extreme value is returned.
  66.      * </p>
  67.      *
  68.      * @param x the value.
  69.      * @return the error function erf(x)
  70.      * @throws org.hipparchus.exception.MathIllegalStateException
  71.      * if the algorithm fails to converge.
  72.      * @see Gamma#regularizedGammaP(double, double, double, int)
  73.      */
  74.     public static double erf(double x) {
  75.         if (FastMath.abs(x) > 40) {
  76.             return x > 0 ? 1 : -1;
  77.         }
  78.         final double ret = Gamma.regularizedGammaP(0.5, x * x, 1.0e-15, 10000);
  79.         return x < 0 ? -ret : ret;
  80.     }

  82.     /**
  83.      * Returns the error function.
  84.      *
  85.      * \[
  86.      * \mathrm{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{t=0}^x e^{-t^2}dt
  87.      * \]
  88.      *
  89.      * <p>This implementation computes erf(x) using the
  90.      * {@link Gamma#regularizedGammaP(double, double, double, int) regularized gamma function},
  91.      * following <a href="http://mathworld.wolfram.com/Erf.html"> Erf</a>, equation (3)</p>
  92.      *
  93.      * <p>The value returned is always between -1 and 1 (inclusive).
  94.      * If {@code abs(x) > 40}, then {@code erf(x)} is indistinguishable from
  95.      * either 1 or -1 as a double, so the appropriate extreme value is returned.
  96.      * </p>
  97.      *
  98.      * @param <T> type of the field elements
  99.      * @param x the value.
  100.      * @return the error function erf(x)
  101.      * @throws org.hipparchus.exception.MathIllegalStateException
  102.      * if the algorithm fails to converge.
  103.      * @see Gamma#regularizedGammaP(double, double, double, int)
  104.      */
  105.     public static <T extends CalculusFieldElement<T>> T erf(T x) {
  106.         final Field<T> field = x.getField();
  107.         final T one = field.getOne();

  109.         if (FastMath.abs(x.getReal()) > 40) {
  110.             return x.getReal() > 0 ? one : one.negate();
  111.         }
  112.         final T ret = Gamma.regularizedGammaP(one.newInstance(0.5), x.square(), 1.0e-15, 10000);
  113.         return x.getReal() < 0 ? ret.negate() : ret;
  114.     }


  117.     /**
  118.      * Returns the complementary error function.
  119.      *
  120.      * \[
  121.      * \mathrm{erfc}(x) =  \frac{2}{\sqrt{\pi}} \int_{t=x}^\infty e^{-t^2}dt = 1 - \mathrm{erf}
  122.      *
  123.      * <p>This implementation computes erfc(x) using the
  124.      * {@link Gamma#regularizedGammaQ(double, double, double, int) regularized gamma function},
  125.      * following <a href="http://mathworld.wolfram.com/Erf.html"> Erf</a>, equation (3).</p>
  126.      *
  127.      * <p>The value returned is always between 0 and 2 (inclusive).
  128.      * If {@code abs(x) > 40}, then {@code erf(x)} is indistinguishable from
  129.      * either 0 or 2 as a double, so the appropriate extreme value is returned.
  130.      * </p>
  131.      *
  132.      * @param x the value
  133.      * @return the complementary error function erfc(x)
  134.      * @throws org.hipparchus.exception.MathIllegalStateException
  135.      * if the algorithm fails to converge.
  136.      * @see Gamma#regularizedGammaQ(double, double, double, int)
  137.      */
  138.     public static double erfc(double x) {
  139.         if (FastMath.abs(x) > 40) {
  140.             return x > 0 ? 0 : 2;
  141.         }
  142.         final double ret = Gamma.regularizedGammaQ(0.5, x * x, 1.0e-15, 10000);
  143.         return x < 0 ? 2 - ret : ret;
  144.     }

  146.     /**
  147.      * Returns the complementary error function.
  148.      *
  149.      * \[
  150.      * erfc(x) = \frac{2}{\sqrt{\pi}} \int_x^\infty e^{-t^2}dt = 1 - erf(x)
  151.      * \]
  152.      *
  153.      * <p>This implementation computes erfc(x) using the
  154.      * {@link Gamma#regularizedGammaQ(double, double, double, int) regularized gamma function}, following <a
  155.      * href="http://mathworld.wolfram.com/Erf.html"> Erf</a>, equation (3).</p>
  156.      *
  157.      * <p>The value returned is always between 0 and 2 (inclusive).
  158.      * If {@code abs(x) > 40}, then {@code erf(x)} is indistinguishable from either 0 or 2 as a double, so the
  159.      * appropriate extreme value is returned. <b>This implies that the current implementation does not allow the use of
  160.      * {@link org.hipparchus.dfp.Dfp Dfp} with extended precision.</b>
  161.      * </p>
  162.      *
  163.      * @param x the value
  164.      * @param <T> type of the field elements
  165.      *
  166.      * @return the complementary error function erfc(x)
  167.      *
  168.      * @throws org.hipparchus.exception.MathIllegalStateException if the algorithm fails to converge.
  169.      * @see Gamma#regularizedGammaQ(double, double, double, int)
  170.      */
  171.     public static <T extends CalculusFieldElement<T>> T erfc(T x) {
  172.         final Field<T> field = x.getField();
  173.         final T        zero  = field.getZero();
  174.         final T        one   = field.getOne();

  176.         if (FastMath.abs(x.getReal()) > 40) {
  177.             return x.getReal() > 0 ? zero : one.newInstance(2.);
  178.         }
  179.         final T ret = Gamma.regularizedGammaQ(one.newInstance(0.5), x.square(), 1.0e-15, 10000);
  180.         return x.getReal() < 0 ? ret.negate().add(2.) : ret;
  181.     }

  183.     /**
  184.      * Returns the difference between erf(x1) and erf(x2).
  185.      * <p>
  186.      * The implementation uses either erf(double) or erfc(double)
  187.      * depending on which provides the most precise result.
  188.      *
  189.      * @param x1 the first value
  190.      * @param x2 the second value
  191.      * @return erf(x2) - erf(x1)
  192.      */
  193.     public static double erf(double x1, double x2) {
  194.         if(x1 > x2) {
  195.             return -erf(x2, x1);
  196.         }

  198.         return
  199.         x1 < -X_CRIT ?
  200.             x2 < 0.0 ?
  201.                 erfc(-x2) - erfc(-x1) :
  202.                 erf(x2) - erf(x1) :
  203.             x2 > X_CRIT && x1 > 0.0 ?
  204.                 erfc(x1) - erfc(x2) :
  205.                 erf(x2) - erf(x1);
  206.     }

  208.     /**
  209.      * Returns the difference between erf(x1) and erf(x2).
  210.      * <p>
  211.      * The implementation uses either erf(double) or erfc(double)
  212.      * depending on which provides the most precise result.
  213.      *
  214.      * @param x1 the first value
  215.      * @param x2 the second value
  216.      * @param <T> type of the field elements
  217.      *
  218.      * @return erf(x2) - erf(x1)
  219.      */
  220.     public static <T extends CalculusFieldElement<T>> T erf(T x1, T x2) {

  222.         if (x1.getReal() > x2.getReal()) {
  223.             return erf(x2, x1).negate();
  224.         }

  226.         return
  227.                 x1.getReal() < -X_CRIT ?
  228.                         x2.getReal() < 0.0 ?
  229.                                 erfc(x2.negate()).subtract(erfc(x1.negate())) :
  230.                                 erf(x2).subtract(erf(x1)) :
  231.                         x2.getReal() > X_CRIT && x1.getReal() > 0.0 ?
  232.                                 erfc(x1).subtract(erfc(x2)) :
  233.                                 erf(x2).subtract(erf(x1));
  234.     }

  236.     /**
  237.      * Returns the inverse erf.
  238.      * <p>
  239.      * This implementation is described in the paper:
  240.      * <a href="http://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf">Approximating
  241.      * the erfinv function</a> by Mike Giles, Oxford-Man Institute of Quantitative Finance,
  242.      * which was published in GPU Computing Gems, volume 2, 2010.
  243.      * The source code is available <a href="http://gpucomputing.net/?q=node/1828">here</a>.
  244.      * </p>
  245.      * @param x the value
  246.      * @return t such that x = erf(t)
  247.      */
  248.     public static double erfInv(final double x) {

  250.         // beware that the logarithm argument must be
  251.         // computed as (1.0 - x) * (1.0 + x),
  252.         // it must NOT be simplified as 1.0 - x * x as this
  253.         // would induce rounding errors near the boundaries +/-1
  254.         double w = - FastMath.log((1.0 - x) * (1.0 + x));
  255.         double p;

  257.         if (w < 6.25) {
  258.             w -= 3.125;
  259.             p =  -3.6444120640178196996e-21;
  260.             p =   -1.685059138182016589e-19 + p * w;
  261.             p =   1.2858480715256400167e-18 + p * w;
  262.             p =    1.115787767802518096e-17 + p * w;
  263.             p =   -1.333171662854620906e-16 + p * w;
  264.             p =   2.0972767875968561637e-17 + p * w;
  265.             p =   6.6376381343583238325e-15 + p * w;
  266.             p =  -4.0545662729752068639e-14 + p * w;
  267.             p =  -8.1519341976054721522e-14 + p * w;
  268.             p =   2.6335093153082322977e-12 + p * w;
  269.             p =  -1.2975133253453532498e-11 + p * w;
  270.             p =  -5.4154120542946279317e-11 + p * w;
  271.             p =    1.051212273321532285e-09 + p * w;
  272.             p =  -4.1126339803469836976e-09 + p * w;
  273.             p =  -2.9070369957882005086e-08 + p * w;
  274.             p =   4.2347877827932403518e-07 + p * w;
  275.             p =  -1.3654692000834678645e-06 + p * w;
  276.             p =  -1.3882523362786468719e-05 + p * w;
  277.             p =    0.0001867342080340571352 + p * w;
  278.             p =  -0.00074070253416626697512 + p * w;
  279.             p =   -0.0060336708714301490533 + p * w;
  280.             p =      0.24015818242558961693 + p * w;
  281.             p =       1.6536545626831027356 + p * w;
  282.         } else if (w < 16.0) {
  283.             w = FastMath.sqrt(w) - 3.25;
  284.             p =   2.2137376921775787049e-09;
  285.             p =   9.0756561938885390979e-08 + p * w;
  286.             p =  -2.7517406297064545428e-07 + p * w;
  287.             p =   1.8239629214389227755e-08 + p * w;
  288.             p =   1.5027403968909827627e-06 + p * w;
  289.             p =   -4.013867526981545969e-06 + p * w;
  290.             p =   2.9234449089955446044e-06 + p * w;
  291.             p =   1.2475304481671778723e-05 + p * w;
  292.             p =  -4.7318229009055733981e-05 + p * w;
  293.             p =   6.8284851459573175448e-05 + p * w;
  294.             p =   2.4031110387097893999e-05 + p * w;
  295.             p =   -0.0003550375203628474796 + p * w;
  296.             p =   0.00095328937973738049703 + p * w;
  297.             p =   -0.0016882755560235047313 + p * w;
  298.             p =    0.0024914420961078508066 + p * w;
  299.             p =   -0.0037512085075692412107 + p * w;
  300.             p =     0.005370914553590063617 + p * w;
  301.             p =       1.0052589676941592334 + p * w;
  302.             p =       3.0838856104922207635 + p * w;
  303.         } else if (!Double.isInfinite(w)) {
  304.             w = FastMath.sqrt(w) - 5.0;
  305.             p =  -2.7109920616438573243e-11;
  306.             p =  -2.5556418169965252055e-10 + p * w;
  307.             p =   1.5076572693500548083e-09 + p * w;
  308.             p =  -3.7894654401267369937e-09 + p * w;
  309.             p =   7.6157012080783393804e-09 + p * w;
  310.             p =  -1.4960026627149240478e-08 + p * w;
  311.             p =   2.9147953450901080826e-08 + p * w;
  312.             p =  -6.7711997758452339498e-08 + p * w;
  313.             p =   2.2900482228026654717e-07 + p * w;
  314.             p =  -9.9298272942317002539e-07 + p * w;
  315.             p =   4.5260625972231537039e-06 + p * w;
  316.             p =  -1.9681778105531670567e-05 + p * w;
  317.             p =   7.5995277030017761139e-05 + p * w;
  318.             p =  -0.00021503011930044477347 + p * w;
  319.             p =  -0.00013871931833623122026 + p * w;
  320.             p =       1.0103004648645343977 + p * w;
  321.             p =       4.8499064014085844221 + p * w;
  322.         } else {
  323.             // this branch does not appears in the original code, it
  324.             // was added because the previous branch does not handle
  325.             // x = +/-1 correctly. In this case, w is positive infinity
  326.             // and as the first coefficient (-2.71e-11) is negative.
  327.             // Once the first multiplication is done, p becomes negative
  328.             // infinity and remains so throughout the polynomial evaluation.
  329.             // So the branch above incorrectly returns negative infinity
  330.             // instead of the correct positive infinity.
  331.             p = Double.POSITIVE_INFINITY;
  332.         }

  334.         return p * x;

  336.     }

  338.     /**
  339.      * Returns the inverse erf.
  340.      * <p>
  341.      * This implementation is described in the paper:
  342.      * <a href="http://people.maths.ox.ac.uk/gilesm/files/gems_erfinv.pdf">Approximating
  343.      * the erfinv function</a> by Mike Giles, Oxford-Man Institute of Quantitative Finance,
  344.      * which was published in GPU Computing Gems, volume 2, 2010.
  345.      * The source code is available <a href="http://gpucomputing.net/?q=node/1828">here</a>.
  346.      * </p>
  347.      * @param <T> type of the filed elements
  348.      * @param x the value
  349.      * @return t such that x = erf(t)
  350.      */
  351.     public static <T extends CalculusFieldElement<T>> T erfInv(final T x) {
  352.         final T one = x.getField().getOne();

  354.         // beware that the logarithm argument must be
  355.         // computed as (1.0 - x) * (1.0 + x),
  356.         // it must NOT be simplified as 1.0 - x * x as this
  357.         // would induce rounding errors near the boundaries +/-1
  358.         T w = one.subtract(x).multiply(one.add(x)).log().negate();
  359.         T p;

  361.         if (w.getReal() < 6.25) {
  362.             w = w.subtract(3.125);
  363.             p = one.newInstance(-3.6444120640178196996e-21);
  364.             p = p.multiply(w).add(-1.685059138182016589e-19);
  365.             p = p.multiply(w).add(1.2858480715256400167e-18);
  366.             p = p.multiply(w).add(1.115787767802518096e-17);
  367.             p = p.multiply(w).add(-1.333171662854620906e-16);
  368.             p = p.multiply(w).add(2.0972767875968561637e-17);
  369.             p = p.multiply(w).add(6.6376381343583238325e-15);
  370.             p = p.multiply(w).add(-4.0545662729752068639e-14);
  371.             p = p.multiply(w).add(-8.1519341976054721522e-14);
  372.             p = p.multiply(w).add(2.6335093153082322977e-12);
  373.             p = p.multiply(w).add(-1.2975133253453532498e-11);
  374.             p = p.multiply(w).add(-5.4154120542946279317e-11);
  375.             p = p.multiply(w).add(1.051212273321532285e-09);
  376.             p = p.multiply(w).add(-4.1126339803469836976e-09);
  377.             p = p.multiply(w).add(-2.9070369957882005086e-08);
  378.             p = p.multiply(w).add(4.2347877827932403518e-07);
  379.             p = p.multiply(w).add(-1.3654692000834678645e-06);
  380.             p = p.multiply(w).add(-1.3882523362786468719e-05);
  381.             p = p.multiply(w).add(0.0001867342080340571352);
  382.             p = p.multiply(w).add(-0.00074070253416626697512);
  383.             p = p.multiply(w).add(-0.0060336708714301490533);
  384.             p = p.multiply(w).add(0.24015818242558961693);
  385.             p = p.multiply(w).add(1.6536545626831027356);
  386.         }
  387.         else if (w.getReal() < 16.0) {
  388.             w = w.sqrt().subtract(3.25);
  389.             p = one.newInstance(2.2137376921775787049e-09);
  390.             p = p.multiply(w).add(9.0756561938885390979e-08);
  391.             p = p.multiply(w).add(-2.7517406297064545428e-07);
  392.             p = p.multiply(w).add(1.8239629214389227755e-08);
  393.             p = p.multiply(w).add(1.5027403968909827627e-06);
  394.             p = p.multiply(w).add(-4.013867526981545969e-06);
  395.             p = p.multiply(w).add(2.9234449089955446044e-06);
  396.             p = p.multiply(w).add(1.2475304481671778723e-05);
  397.             p = p.multiply(w).add(-4.7318229009055733981e-05);
  398.             p = p.multiply(w).add(6.8284851459573175448e-05);
  399.             p = p.multiply(w).add(2.4031110387097893999e-05);
  400.             p = p.multiply(w).add(-0.0003550375203628474796);
  401.             p = p.multiply(w).add(0.00095328937973738049703);
  402.             p = p.multiply(w).add(-0.0016882755560235047313);
  403.             p = p.multiply(w).add(0.0024914420961078508066);
  404.             p = p.multiply(w).add(-0.0037512085075692412107);
  405.             p = p.multiply(w).add(0.005370914553590063617);
  406.             p = p.multiply(w).add(1.0052589676941592334);
  407.             p = p.multiply(w).add(3.0838856104922207635);
  408.         }
  409.         else if (!w.isInfinite()) {
  410.             w = w.sqrt().subtract(5.0);
  411.             p = one.newInstance(-2.7109920616438573243e-11);
  412.             p = p.multiply(w).add(-2.5556418169965252055e-10);
  413.             p = p.multiply(w).add(1.5076572693500548083e-09);
  414.             p = p.multiply(w).add(-3.7894654401267369937e-09);
  415.             p = p.multiply(w).add(7.6157012080783393804e-09);
  416.             p = p.multiply(w).add(-1.4960026627149240478e-08);
  417.             p = p.multiply(w).add(2.9147953450901080826e-08);
  418.             p = p.multiply(w).add(-6.7711997758452339498e-08);
  419.             p = p.multiply(w).add(2.2900482228026654717e-07);
  420.             p = p.multiply(w).add(-9.9298272942317002539e-07);
  421.             p = p.multiply(w).add(4.5260625972231537039e-06);
  422.             p = p.multiply(w).add(-1.9681778105531670567e-05);
  423.             p = p.multiply(w).add(7.5995277030017761139e-05);
  424.             p = p.multiply(w).add(-0.00021503011930044477347);
  425.             p = p.multiply(w).add(-0.00013871931833623122026);
  426.             p = p.multiply(w).add(1.0103004648645343977);
  427.             p = p.multiply(w).add(4.8499064014085844221);
  428.         }
  429.         else {
  430.             // this branch does not appear in the original code, it
  431.             // was added because the previous branch does not handle
  432.             // x = +/-1 correctly. In this case, w is positive infinity
  433.             // and as the first coefficient (-2.71e-11) is negative.
  434.             // Once the first multiplication is done, p becomes negative
  435.             // infinity and remains so throughout the polynomial evaluation.
  436.             // So the branch above incorrectly returns negative infinity
  437.             // instead of the correct positive infinity.
  438.             p = one.multiply(Double.POSITIVE_INFINITY);
  439.         }

  441.         return p.multiply(x);

  443.     }

  445.     /**
  446.      * Returns the inverse erfc.
  447.      * @param x the value
  448.      * @return t such that x = erfc(t)
  449.      */
  450.     public static double erfcInv(final double x) {
  451.         return erfInv(1 - x);
  452.     }

  454.     /**
  455.      * Returns the inverse erfc.
  456.      * @param x the value
  457.      * @param <T> type of the field elements
  458.      * @return t such that x = erfc(t)
  459.      */
  460.     public static <T extends CalculusFieldElement<T>> T erfcInv(final T x) {
  461.         return erfInv(x.negate().add(1));
  462.     }

  464. }
Created with JaCoCo 0.8.12.202403310830