Repository: scipy/scipy
Pull request: BUG: special: erfinv(x<<1) loses precision
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
http://scipy.github.io/devdocs/dev/contributor/development_workflow.html#checklist-before-submitting-a-pr

Also, please name and describe your PR as you would write a
commit message:
https://docs.scipy.org/doc/numpy/dev/development_workflow.html#writing-the-commit-message

Note that we are a team of volunteers; we appreciate your
patience during the review process.

Again, thanks for contributing!
-->

#### Reference issue
<!--Example: Closes gh-WXYZ.-->

closes gh-12758

#### What does this implement/fix?
<!--Please explain your changes.-->

Avoid the loss of precision for small arguments by using the Taylor series for small x.

#### Additional information
<!--Any additional information you think is important.-->

At small arguments, `erf(x) = 2/\sqrt{\pi} (x - x^3 / 3 + ...)`, so we only keep the linear term for `|x| < 1e-7` and then the inverse function `erfinv(x) = x * \sqrt{\pi} / 2`. 

Problem statement:
BUG: special: erfinv(y) is inaccurate for small abs(y).
`erfinv(y)` is implemented as `ndtri(0.5 * (y+1)) * NPY_SQRT1_2` (where `NPY_SQRT1_2` is `1/sqrt(2)`).  This limits the precision that can be obtained for small `abs(y)`, because for any `y` such that `abs(y)` is less than ~1e-16, `y + 1` is `1`.  So, for example, `erfinv(1e-20)` returns 0, but the correct result is `8.86226925452758e-21`:
```
In [45]: erfinv(1e-20)  # incorrect result                                                          
Out[45]: 0.0

In [46]: erf(8.86226925452758e-21)                                              
Out[46]: 1e-20
```
Octave and Boost return the correct value; we should too.

Octave:
```
octave:3> format long
octave:4> erfinv(1e-20)
ans =    8.86226925452758e-21
```
Boost
```
#include <iostream>
#include <boost/math/special_functions/erf.hpp>

using boost::math::erf_inv;
using namespace std;


int main()
{
    double p = 1e-20;
    double x = erf_inv(p);
    cout << scientific << setprecision(15) << x << endl;
    // 8.862269254527580e-21
}
```

Hints from issue discussion:

Commit message: TST: special: roundtrip erf(erfinv(x<<1))

Patch to review:
diff --git a/scipy/special/cephes/erfinv.c b/scipy/special/cephes/erfinv.c
index 4bdbc5c725aa..28f180f7d59f 100644
--- a/scipy/special/cephes/erfinv.c
+++ b/scipy/special/cephes/erfinv.c
@@ -17,6 +17,17 @@ double erfinv(double y) {
     const double domain_lb = -1;
     const double domain_ub = 1;
 
+    const double thresh = 1e-7;
+
+    /* 
+     * For small arguments, use the Taylor expansion
+     * erf(y) = 2/\sqrt{\pi} (y - y^3 / 3 + O(y^5)),    y\to 0 
+     * where we only retain the linear term.
+     * Otherwise, y + 1 loses precision for |y| << 1.
+     */
+    if ((-thresh < y) && (y < thresh)){
+        return y / M_2_SQRTPI;
+    } 
     if ((domain_lb < y) && (y < domain_ub)) {
         return ndtri(0.5 * (y+1)) * NPY_SQRT1_2;
     }
diff --git a/scipy/special/tests/test_erfinv.py b/scipy/special/tests/test_erfinv.py
index 61ad74315841..5792d320db6f 100644
--- a/scipy/special/tests/test_erfinv.py
+++ b/scipy/special/tests/test_erfinv.py
@@ -57,3 +57,19 @@ def test_literal_values(self):
     )
     def test_domain_bounds(self, f, x, y):
         assert_equal(f(x), y)
+
+    def test_erfinv_asympt(self):
+        # regression test for gh-12758: erfinv(x) loses precision at small x
+        # This test checks erfinv against the Taylor expansion:
+        # erf(x) = 2/\sqrt{\pi} (x - x^3 / 3 + O(x^5)),    x\to 0 
+        # where we only retain the linear term.
+        x = np.array([1e-20, 1e-15, 1e-14, 1e-10, 1e-8, 0.9e-7, 1.1e-7, 1e-6])
+        assert_allclose(sc.erfinv(x),
+                        np.sqrt(np.pi)/2 * x,
+                        rtol=1e-10)
+
+        # also test the roundtrip consistency
+        assert_allclose(sc.erf(sc.erfinv(x)),
+                        x,
+                        rtol=1e-10)
+


Code review comments:
diff_hunk: @@ -57,3 +57,19 @@ def test_literal_values(self):
     )
     def test_domain_bounds(self, f, x, y):
         assert_equal(f(x), y)
+
+    def test_erfinv_asympt(self):
+        # regression test for gh-12758: erfinv(x) loses precision at small x
+        # This test checks erfinv against the Taylor expansion:
+        # erf(x) = 2/\sqrt{\pi} (x - x^3 / 3 + O(x^5)),    x\to 0 
+        # where we only retain the linear term.
+        x = np.array([1e-20, 1e-15, 1e-14, 1e-10, 1e-8, 0.9e-7, 1.1e-7, 1e-6])
+        assert_allclose(sc.erfinv(x),
+                        np.sqrt(np.pi)/2 * x,
path: scipy/special/tests/test_erfinv.py
text: @user1:
```suggestion
        expected = np.array([8.86226925452758e-21,
                             8.862269254527581e-16,
                             8.86226925452758e-15,
                             8.862269254527581e-11,
                             8.86226925452758e-09,
                             7.97604232907484e-08,
                             9.74849617998037e-08,
                             8.8622692545299e-07])
        assert_allclose(sc.erfinv(x), expected,
```

@author:
Oops, sorry, I pushed a commit before seeing the suggested change. I can roll back the commit and accept it, if you want me to.

@author:
Oops, no, I cannot accept it anymore. Grrr. 

@author:
FTR, local rebase to remove the commit + force-push restores the 'commit suggestion' button. Done.