Repository: scipy/scipy
Pull request: BUG: special: Fix hyperu for b = 1, small x, small a
<!-- 
Thanks for contributing a pull request! Please ensure that
your PR satisfies the checklist before submitting:
https://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://scipy.github.io/devdocs/dev/contributor/development_workflow.html#writing-the-commit-message.
However, please only include an issue number in the description, not the title,
and please ensure that any code names containing underscores are enclosed in backticks.

Depending on your changes, you can skip CI operations and save time and energy: 
https://scipy.github.io/devdocs/dev/contributor/continuous_integration.html#skipping

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 #15650 
Ticks a box on https://github.com/scipy/scipy/issues/20223
#### What does this implement/fix?
<!--Please explain your changes.-->
This PR uses DLMF 13.3.7 as suggested by @mdhaber in https://github.com/scipy/scipy/issues/15650#issuecomment-2323741880 to improve accuracy of `hyperu` when `b == 1` and `x` and `a` are close to zero. On `main`, `hyperu` is returning garbage values for such input. I think the underlying issue is that Zhang and Jin weren't careful when implementing the hypergeometric functions, and sometimes used formula outside of the domains in which they are valid. @mdhaber proposed to use @mdhaber when `b == 1`, `|x| < 1` and `|a| < 1`, based on numerical experiments described in the linked comment above, and I agree with his findings. I've added test cases covering this region.

#### Additional information
<!--Any additional information you think is important.-->
`hyperu` involves wrappers on top of wrappers. There is a wrapper `hyperu` in `_hypergeometric.pxd` which calls a wrapper `hypU_wrap` in `xsf/specfun.h` which calls `chgu` in `xsf/specfun/specfun.h`. The lowest level is ilhan's direct translation of Zhang and Jin's Fortran code, the middle layer is based on the old `specfun_wrappers` which called functions from `specfun.f`, and the outer Cython wrapper has some fixes for special cases patched in. I chose to make @mdhaber's suggested fix in the outer Cython wrapper, which is the most convenient spot. 

Obviously, this wrappers on top of wrappers architecture is not ideal. The plan is to eventually rewrite `hyperu` along the lines of https://github.com/scipy/scipy/pull/20089, and lift acceptably accurate functions up out of `xsf/specfun/specfun.h` as was done for exponential integrals here: https://github.com/scipy/scipy/pull/21210/files#diff-5fbee505531722bdef130b290e08018cb6291bfb85a3a84810c155c28e496210.


Problem statement:
BUG: scipy.special.hyperu(a,b,x) gives wrong values for a close to 0 and b=1, x close to zero
### Describe your issue.

Hello,

I noticed that the special function hyperu(a,b,x) gives output values that are extremely off for some x close to 0, when a is close to 0 and b=1.
In the code example below, I compared the values with mpmath's hyperu. For a certain interval of input values x the output is off by several magnitudes (in this example the errors appear near x=0.01). I see no such errors for a's that are separated sufficiently from 0, e.g. |a|>0.5.



### Reproducing Code Example

```python
import numpy as np
import scipy
import matplotlib.pyplot as plt
import scipy.special
from mpmath import *

a = 0.1
x = np.linspace(0.001,0.05,1001)
y1 = np.zeros((1001,))
y2 = np.zeros((1001,))

for k in range(len(x)):
    y1[k] = scipy.special.hyperu(a,1,x[k])
    y2[k] = mp.hyperu(a,1,x[k])

print(scipy.special.hyperu(a,1,0.01))
print(mp.hyperu(a,1,0.01))

plt.plot(x,y1,'r')
plt.plot(x,y2,'b')
plt.show()
```


### Error message

```shell
no error message, just a miscalculation
```


### SciPy/NumPy/Python version information

('1.1.0', '1.16.2', sys.version_info(major=2, minor=7, micro=16, releaselevel='final', serial=0))

Hints from issue discussion:
upgraded versions:
1.7.3 1.21.5 sys.version_info(major=3, minor=7, micro=3, releaselevel='final', serial=0)

don't appear to give correct results either 
Hi @matthias-baur, thanks for the bug report. I haven't looked into the details yet, but we're aware that there are some problems with the hypergeometric function implementations in `scipy/special/specfun/specfun.f`. Numerical evaluation of hypergeometric functions involves handling a patchwork of special cases.  We're in the process of re-implementing SciPy's hypergeometric functions in Cython one special case at a time, falling back to the original Fortran implementations for cases that aren't handled yet. I've been working on this for the past 7 months, focusing on the Gaussian hypergeometric function `hyp2f1`, but it's a slow process due to limited spare time on my part. I plan to start working on fixing any problems that may exist in  `hyperu` and  hope to have some fixes ready for version 1.9.
Thanks for the response, that's good to hear! I hope this bug can be sorted out without too much hassle.
I looked into this a bit. Here's what the $\log_{10}$ of the error looks like as a function of $\log_{10}a$ and $\log_{10}z$ at $b=1$.

<details>

```python3
import numpy as np
from mpmath import mp
from scipy import special
import matplotlib.pyplot as plt

def relative_error(res, ref):
    return np.where(res != ref, abs(res - ref) / abs(ref), 0)

@np.vectorize
def hyperu_ref(a, b, z):
    return float(mp.hyperu(a, b, z))

def plot(res, ref):
    err = relative_error(res, ref)
    logerr = np.log10(err)
    logerr[logerr < -15] = -15
    plt.contourf(np.log10(a).ravel(), np.log10(z).ravel(),
                 logerr, cmap='viridis')
    plt.xlabel("log10(a)")
    plt.ylabel("log10(z)")
    # plt.title(title)
    plt.colorbar()
    plt.show()

a = np.logspace(-16, 0, 30)
b = 1
z = np.logspace(-16, 0, 25)[:, np.newaxis]

res = special.hyperu(a, b, z)
ref = hyperu_ref(a, b, z)
plot(res, ref)
```

</details>

<img width="875" alt="image" src="https://github.com/user-attachments/assets/e51c7e1a-970e-4597-bbab-659153529c8f">

The plot looks similar when we replace $a$ with $-a$. Note that the relative error outside the problem region is extremely low, like `1e-14`.

$b = 1$ appears to have a special implementation because if we make it $b = 1 + 10^{-15}$, we get a totally different plot.

<img width="772" alt="image" src="https://github.com/user-attachments/assets/dfe328a2-3745-4b84-a111-5b4ae877f58b">

That error is not great either, but that's a different story.

It would be ideal to fix `hyperu` everywhere, but that's a tall order, and it may take some time. It also might be nice to find the real problem and fix it, but that also hasn't happened in recent years. In the meantime, we could patch the problem for $b=1$ when $z < 1$ and $|a| < 1$ using [DLMF 13.3.7](https://dlmf.nist.gov/13.3).

```python3
def hyperu_b1(a, b, z):
    # Using DLMF Recurrence Relation 13.3.7
    a = a + 1
    return -(b - 2*a - z)*special.hyperu(a, b, z) - a*(a-b+1)*special.hyperu(a+1, b, z)
```

The problem goes away, and the error is on par with what it should be for the usual $b=1$ case.

<img width="809" alt="image" src="https://github.com/user-attachments/assets/217c97ac-4af1-45c8-b90b-a5ce55d215b7">

Thanks @mdhaber! Very helpful. I'll follow up this week with a PR with your fix.

Commit message: Add test cases for gh-15650

Patch to review:
diff --git a/scipy/special/_hypergeometric.pxd b/scipy/special/_hypergeometric.pxd
index d3b4b6ea7ca5..4eeafd4c7a68 100644
--- a/scipy/special/_hypergeometric.pxd
+++ b/scipy/special/_hypergeometric.pxd
@@ -26,5 +26,9 @@ cdef inline double hyperu(double a, double b, double x) noexcept nogil:
         else:
             # DLMF 13.2.14-15 and 13.2.19-21
             return cephes_poch_wrap(1.0 - b + a, -a)
+    if b == 1 and abs(a) < 1 and abs(x) < 1:
+        # DLMF 13.3.7. Fixes gh-15650
+        a = a + 1
+        return -(b - 2*a - x)*hypU_wrap(a, b, x) - a*(a - b + 1)*hypU_wrap(a + 1, b, x)
 
     return hypU_wrap(a, b, x)
diff --git a/scipy/special/tests/test_hypergeometric.py b/scipy/special/tests/test_hypergeometric.py
index 77c8454192d5..a7541e4689bd 100644
--- a/scipy/special/tests/test_hypergeometric.py
+++ b/scipy/special/tests/test_hypergeometric.py
@@ -23,6 +23,41 @@ def test_special_cases(self):
     def test_nan_inputs(self, a, b, x):
         assert np.isnan(sc.hyperu(a, b, x)) == np.any(np.isnan([a, b, x]))
 
+    @pytest.mark.parametrize(
+        'a,b,x,expected',
+        [(1e-15, 1.0, 1e-15, 1.0000000000000346),
+         (1e-15, 1.0, 4.6415888336127916e-11, 1.0000000000000238),
+         (1e-15, 1.0, 2.1544346900318865e-06, 1.000000000000013),
+         (1e-15, 1.0, 0.1, 1.0000000000000022),
+         (4.6415888336127916e-11, 1.0, 1e-15, 1.000000001603148),
+         (4.6415888336127916e-11, 1.0, 4.6415888336127916e-11, 1.0000000011043908),
+         (4.6415888336127916e-11, 1.0, 2.1544346900318865e-06, 1.0000000006056338),
+         (4.6415888336127916e-11, 1.0, 0.1, 1.0000000001068765),
+         (2.1544346900318865e-06, 1.0, 1e-15, 1.0000744116183264),
+         (2.1544346900318865e-06, 1.0, 4.6415888336127916e-11, 1.0000512613332657),
+         (2.1544346900318865e-06, 1.0, 2.1544346900318865e-06, 1.0000281110482052),
+         (2.1544346900318865e-06, 1.0, 0.1, 1.0000049607649177),
+         (0.1, 1.0, 1e-15, 4.6048315083246285),
+         (0.1, 1.0, 4.6415888336127916e-11, 3.475343054552083),
+         (0.1, 1.0, 2.1544346900318865e-06, 2.345854924985846),
+         (0.1, 1.0, 0.1, 1.2202292433349389)]
+    )
+    def test_gh_15650(self, a, b, x, expected):
+        # See https://github.com/scipy/scipy/issues/15650
+        #
+        # Reference values computed with mpmath
+        #
+        # import itertools as it
+        # import numpy as np
+        #
+        # from mpmath import mp
+        #
+        # cases = []
+        # for a, x in it.product(np.logspace(-15, -1, 4), repeat=2):
+        #     with mp.workdps(100):
+        #         cases.append((a, 1.0, x, float(mp.hyperu(a, 1.0, x))))
+        assert_allclose(sc.hyperu(a, b, x), expected, rtol=1e-13)
+
 
 class TestHyp1f1:
 


Code review comments:
diff_hunk: @@ -26,5 +26,9 @@ cdef inline double hyperu(double a, double b, double x) noexcept nogil:
         else:
             # DLMF 13.2.14-15 and 13.2.19-21
             return cephes_poch_wrap(1.0 - b + a, -a)
+    if b == 1 and abs(a) < 1 and abs(x) < 1:
+        # DLMF 13.3.7. Fixes gh-15650
+        a = a + 1
+        return -(b - 2*a - x)*hypU_wrap(a, b, x) - a*(a - b + 1)*hypU_wrap(a + 1, b, x)
path: scipy/special/_hypergeometric.pxd
text: @user1:
Re: https://github.com/scipy/scipy/pull/21495/files#r1740896703, can this be:
```suggestion
        return -(b - 2*a - x)*hyperu(a, b, x) - a*(a - b + 1)*hyperu(a + 1, b, x)
```
? 


@author:
I'm not sure how much impact it actually has on performance, but using `hypU_wrap` saves us from repeating the checks for `nan`, `x <  0`, `x == 0` on the previous lines. 

@author:
I see, if  `a < 0` then the recurrence still kicks us back into the problematic region. Yes, it needs to be `hyperu`.

@user1:
If 0.9 is the threshold, yes. But it doesn't have to be. We just had to pick the threshold more carefully. See below.

diff_hunk: @@ -26,5 +26,9 @@ cdef inline double hyperu(double a, double b, double x) noexcept nogil:
         else:
             # DLMF 13.2.14-15 and 13.2.19-21
             return cephes_poch_wrap(1.0 - b + a, -a)
+    if b == 1 and abs(a) < 1 and abs(x) < 1:
path: scipy/special/_hypergeometric.pxd
text: @user1:
It looked to me like if argument `x` is negative, you get NaN, which is why I didn't mention it. I suppose `abs(x)` is fine, though.

@author:
Just an oversight on my part. You're right; it should just be `x`.

diff_hunk: @@ -23,6 +23,41 @@ def test_special_cases(self):
     def test_nan_inputs(self, a, b, x):
         assert np.isnan(sc.hyperu(a, b, x)) == np.any(np.isnan([a, b, x]))
 
+    @pytest.mark.parametrize(
+        'a,b,x,expected',
+        [(1e-15, 1.0, 1e-15, 1.0000000000000346),
path: scipy/special/tests/test_hypergeometric.py
text: @user1:
I think negative `a` should be included in the tests. We will likely find that there is still a rough edge as `a` approaches `-1` because it will call `hyperu` with `a + 1`, which will be in the problem region. If this can recurse the call  to `hyperu`, I think it would work. 
(If that won't work, we might need to adjusting the criterion `abs(a) < 1` to something like `abs(a) < 0.9`, but it might leave a sliver of the problematic behavior unless we fine-tune the threshold), or we could use the recurrence relation in the other direction, which I'd need to test.)