BUG: Incorrect translation of AMOS Fortran code resulting in incorrect/NaN outputs #46
Description
This following line in the function kscl in amos.h is incorrectly translated from the original Fortran code:
cs = (exp(std::real(cs)) / tol) * (cos(std::imag(cs)) + sin(std::imag(cs) * std::complex<double>(0, 1)));
The original Fortran code (in zkscl.f) is:
STR = DEXP(CSR)/TOL
CSR = STR*DCOS(CSI)
CSI = STR*DSIN(CSI)
so really it should be translated as
cs = (exp(std::real(cs))/tol)*std::complex<double>(cos(std::imag(cs)), sin(std::imag(cs)));
This resulted in incorrect/NaN outputs for special.kv on certain edge case values. For example
>>> from scipy import special
>>> special.kv(0, 680-1000j)
np.complex128(nan+infj)
>>> special.kv(0, 680-680j)
np.complex128(-4.5537300329442675e-298-0.02992162857456993j)
when really the outputs should be 1.9017e-298+1.7134e-297j and -4.5537e-298-1.8787e-297 respectively.