BUG: fix y1 for large inputs - #231
Merged
Merged
Conversation
y1 formed the asymptotic phase as the double-precision subtraction `x - 3*pi/4`, which injects an absolute phase error of order ulp(x)/2 and turns it into relative error in the result. Accuracy degrades from ~1e-13 at x=1e4 to a wrong answer by a factor of 7.7 at x=1e20. Applies the same trigonometric rewrite already used for j0 (scipy#131), y0 (scipy#140) and j1 (scipy#143), which y1 was not carried over to: sin(x - 3pi/4) = -(sin x + cos x)/sqrt(2) cos(x - 3pi/4) = (sin x - cos x)/sqrt(2) so p*sin(x-3pi/4) + w*q*cos(x-3pi/4) = [(w*q-p)*sin x - (p+w*q)*cos x]/sqrt(2), with the 1/sqrt(2) folded into SQRT2OPI -> SQRT1OPI exactly as y0 does. sin and cos are then evaluated on the unmodified x, so no phase precision is lost. Adds the missing y1 case to the large-input tests, which covered j0, y0 and j1 but not y1.
Contributor
|
Thanks @FireflySentinel ! |
dschmitz89
approved these changes
Jul 31, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Reference issue
Follow-up to j0 (#131), y0 (#140) and j1 (#143).
y1is the last one in that group still forming the phase asx - 3pi/4.What does this implement/fix?
xn = x - THPIO4rounds to the nearest double, so it costs aboutulp(x)/2ofabsolute phase, which
sin/costurn into relative error.I implemented the same rewrite as y0 (#140), so the trig is evaluated on the unmodified
x. I also adds themissing
y1case totest_bessel_functions_large_inputs.cpp.Additional information
Relative error vs mpmath at 1000 dps:
j0,y0andj1are at ~1e-16 on the same points, so this isy1only.Additionally,
yv/yninherit it through the recurrence. For example,yv(200, 1e7)is off by 4.1e-10 before this change.pixi run testspasses and the new test fails on main.AI Generation Disclosure
I used Claude Code to draft the fix and test.