BUG: amos/amos.h: fix translation bugs in binu and rati#138
Conversation
|
Thanks for digging so deep into this difficult code. Tests should go into this folder. You could mimic the approach from the other tests in that folder. With a little help from the LLM of your choice, it should hopefully not be too difficult (at least that was my experience as someone who never used ctest before). Do you have any specific questions? |
| if (nw > 0) { | ||
| return -1; | ||
| nz = -1; | ||
| if (nw == -2) { | ||
| nz = -2; | ||
| } | ||
| return nz; |
There was a problem hiding this comment.
The if condition in 2792 is unreachable as first we check for nw>0 and then for nw==-2. Was this in the original Fortran source?
There was a problem hiding this comment.
Thanks, of course you are right, I reverted this change.
Yes, this actually was this way in the Fortran source: https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/special/amos/zbinu.f#L99
|
Thank you for the review! I added a test which fails on main and passes here. Mainly, I was worried about the size of the reference data. Here, I could reduce it to an array of length 15. However, in one of my other PRs, the length would be larger than 200 and I am not sure, if this is nice to put into the code. When this is merged, I will update #92 and #93 and add tests to the same file. |
| if ((az <= 2.) || (az * az * 0.25 <= (dfnu + 1.0))) { | ||
| /* GOTO 10 */ | ||
| nw = seri(z, fnu, kode, n, cy, tol, elim, alim); | ||
| nw = seri(z, fnu, kode, nn, cy, tol, elim, alim); |
There was a problem hiding this comment.
| // MILLER ALGORITHM NORMALIZED BY THE SERIES | ||
| // | ||
| nw = mlri(z, fnu, kode, n, cy, tol); | ||
| nw = mlri(z, fnu, kode, nn, cy, tol); |
There was a problem hiding this comment.
| // ASYMPTOTIC EXPANSION FOR LARGE Z | ||
| // | ||
| nw = asyi(z, fnu, kode, n, cy, rl, tol, elim, alim); | ||
| nw = asyi(z, fnu, kode, nn, cy, rl, tol, elim, alim); |
There was a problem hiding this comment.
| if (az <= rl) { | ||
| /* 70 */ | ||
| nw = mlri(z, fnu, kode, n, cy, tol); | ||
| nw = mlri(z, fnu, kode, nn, cy, tol); |
There was a problem hiding this comment.
|
|
||
| if (p1 == 0.) { | ||
| p1 = std::complex<double>(0, 0); | ||
| p1 = std::complex<double>(tol, tol); |
There was a problem hiding this comment.
|
|
||
| for (i = 2; i < (n + 1); i++) { | ||
| pt = cdfnu + t1 * rz * cy[k]; | ||
| pt = cdfnu + t1 * rz + cy[k]; |
There was a problem hiding this comment.
|
The translation errors look genuine. I am not so convinced about the outcome though. When I run your test data against scipy's main branch, the maximal relative error I get is in the order of 10 eps. import numpy as np
from mpmath import mp
from scipy.special import jv
mp.dps = 1000
n = 15
z = 70 - 0.7j
for n in range(1, n + 1):
reference = complex(mp.besselj(n, mp.mpc(z)))
scipy_result = jv(n, z)
rel_error = abs(reference - scipy_result) / abs(reference)
print(f"n={n}, mpmath.j{n}(z)={reference}, scipy.special.jv({n}, z)={scipy_result}, relative error={rel_error}")Do your changes decrease the relative error even more? |
|
No, I expect no changes for scipy from this PR. Instead, this fixes a feature of the amos library which at the moment cannot be accessed through scipy (which is probably the reason why these bugs were not caught during the translation process). In the test file you can find the parameter |
steppi
left a comment
steppi
left a comment
There was a problem hiding this comment.
Looks good to me. Thanks @JaRoSchm and thanks @dschmitz89 for the review
* BUG: amos/amos.h: fix translation bugs in binu and rati * remove unreachable code * add test
3 participants
Reference issue
What does this implement/fix?
Hi, I detected these errors in the same way as in #92. This time I used
besjwithz = 70.26595687262548 - 0.7026595687262548jand n = 50 in https://github.com/JaRoSchm/pyamos/blob/main/compare_xsf_zbessel.py for finding the bugs. This was simply incorrectly translated from FORTRAN, see https://github.com/scipy/scipy/blob/4edfcaa3ce8a387450b6efce968572def71be089/scipy/special/amos/. Regarding testing the same question/comment from #92 applies.Additional information
AI Generation Disclosure
No AI tools used