scipy: BUG: special: ive(v, z) returns nan for sufficiently large z

@ilayn First, I would like to apologize for not having worked on this issue for some time because of other duties. Motivated by your comment, I looked into the behavior on the imaginary axis of z if the gotos in lines 200 and 201 of zbesi.f are commented out. While on the real axis of z, everything seems to look fine, one can find problematic cases on the imaginary axis. Just as an example, I used the following code

from scipy import special
import mpmath as mp

z = complex(0, 10**20)
nu = 10**10
x1 = special.ive(nu, z)
x2 = mp.besseli(nu, z)
print(f"{z = }")
print(f"{nu = }")
print(f"Scipy:  {x1}")
print(f"mpmath: {x2}")
print(f"relative error: {abs(x1-x2)/abs(x1)}")

resulting in the following output:

z = 1e+20j
nu = 10000000000
Scipy:  (4.399561217393439e-11+2.886868016534155e-17j)
mpmath: (4.39956558504559e-11 + 1.36895178376941e-51j)
relative error: 1.19000355941554e-6

Note that in the case of purely imaginary arguments z, there is no exponential scaling in special.ive. Even though the results are close, a relative error of the order of 10⁻⁶ seems too large to me to be acceptable.

Without this precision problem on the imaginary axis, I would have suggested to remove lines 1673-1676 in your new _amos.c, but in view of the above example, I am pretty reluctant.

@steppi and @WarrenWeckesser Since @dschmitz89 suggested earlier that I contact you concerning a possible modification to zbesi.f: Do you think that the precision in the above example is still good enough in view of the fact that one could significantly extend the range of arguments?

@dschmitz89 : Ok, I will take a look when I find some spare time. Actually, I am more familiar with Fortran, but since we are talking about small modifications, C is fine for me as well.

Returning NaN for ν² > 2|z| is probably not a good idea because beyond the curve ν² = 2|z|, another approach to determine the Bessel function is used and valid results may be produced there. The problem is that the checks are done very early, i.e. in zbesi.f. Only later, in zbinu.f it is decided which algorithm to use. It is there where the call to zasyi.f happens if z is sufficiently large while ν² < 2|z|. A possible solution could be to move the checks to zbinu.f and then to modify the checks for the region where zasyi.f is applied. However, that would go against the logic of the present code.

The problem with precision in the region ν² < 2|z| could be related to the following problem. In the asymptotic expression for the Bessel function, numerators like 4ν²-1 appear. However, for large ν, 4ν²-1 and 4ν² might not differ numerically anymore. This is expected to happen at ν given by 0.5*sqrt(2**53), i.e. around 4.7e7. In https://user-images.githubusercontent.com/1056701/263541623-5d0b255d-0f93-450e-9556-870dcb3ebc32.png the problems happen even at smaller ν, but the error still seems tolerable. For a further analysis, I could implement the asymptotic formula in Python (for experimentation, not for production, so speed is not an issue here) and compare with the same implementation in mpmath. The comparison might point towards problems with precision.

According to Fig. 4 in D.E. Amos, SAND-83-0643, the large-ν expansion you are referring to is used in a different parameter regime, where ν² is much larger than 2|z|, i.e. in the upper right region of Fig. 4. From the example initially discussed. i.e. ive(1, 1e10), I assumed so far, that we are mainly interested in the region on the upper left, i.e. for ν²<2|z|. So far, I did not take a look into the code related to the parameter region where ν² ⪢ 2|z|.

Here are some first data on the relative error for ive relative to the result of mpmath.besseli (of course multiplied with the required exponential). The following plot is outdated. The linear increase of the error is an artefact of insufficient mpmath precision. For more details, see discussion further down. The asymptotic formula for large z, which we are discussing here and which is used in zasyi.f, cannot converge for ν² > 2|z|. That is why the four lines calculated for the indicated values of z end at different values of ν. For z=10¹⁰, the full range has a very good precision. However, for much larger values of z, the precision degrades as ν comes close to its maximally allowed value. Whether this results from an issue with the asymptotic expansion or from the implementation and the finite word size for floats and integers is not clear to me yet.

In the meantime, I have looked a bit more into the code of zasyi.f and can confirm that it implements the asymptotic expansion of the modified Bessel function even though in some places it might not be obvious at first sight. What I have not checked in detail so far are the error and convergence conditions. Only in two places a conversion from a float (FNU in one case and 2*RL in the other) to an integer is done via a single precision float, so there might be a relation to the condition in zbesi.f leading to IERR=4, but I am not sure about it.

I agree that an experimental approach commenting out the following lines in zbesi.f https://github.com/scipy/scipy/blob/e837f21f753e376514520b82053d8f49f290f046/scipy/special/amos/zbesi.f#L200-L201 and comparing with mpmath might give a better idea where the algorithm will fail and I hope that I will find some time in the next few days to work on it.