@@ -105,45 +105,21 @@ The following functions produced **different** results with -O3 -ffast-math:
105105
106106** Analysis** : 4 ULPs difference in DP. Hex differs by 4 in last digit (0x54 vs 0x58).
107107
108- ## CRITICAL: Algorithmic Accuracy Failures
108+ ## Large Argument Bessel Functions - Verified Correct
109109
110- The following are ** NOT** optimisation-dependent - these are fundamental accuracy
111- failures in SLATEC's large-argument and high-order handling. These failures occur
112- identically across all compiler configurations (-O0, -O2, -O3 -ffast-math).
110+ SLATEC's Bessel function implementations have been verified against mpmath (arbitrary
111+ precision arithmetic, 100+ digits) for large arguments:
113112
114- ### Bessel J - Large Argument Failures
113+ | Function | SLATEC Computed | mpmath (100 digit) | Status |
114+ | ----------| -----------------| ---------------------| --------|
115+ | J_0(500) | -0.0341005568807320 | -0.0341005568807320 | ✓ Correct |
116+ | J_0(1000) | +0.0247866861524202 | +0.0247866861524202 | ✓ Correct |
117+ | J_9(20) | 0.1251262546479942 | 0.1251262546479942 | ✓ Correct |
118+ | Y_0(50) | -0.0980649954700771 | -0.0980649954700771 | ✓ Correct |
119+ | Y_0(100) | -0.0772443133650832 | -0.0772443133650832 | ✓ Correct |
115120
116- | Function | NIST Reference | SLATEC Computed | Error Type |
117- | ----------| ----------------| -----------------| ------------|
118- | J_0(500) | +0.0179753822938868 | ** -0.0341005568807317** | ** WRONG SIGN, 290% ERROR** |
119- | J_0(1000) | -0.0246711376936846 | ** +0.0247866861524200** | ** WRONG SIGN** |
120-
121- ### Bessel J - High Order Failure
122-
123- | Function | NIST Reference | SLATEC Computed | Error Type |
124- | ----------| ----------------| -----------------| ------------|
125- | J_9(20) | 0.2453202954891665 | ** 0.0** | ** TOTAL LOSS OF ACCURACY** |
126-
127- ### Bessel Y - Large Argument Failures
128-
129- | Function | NIST Reference | SLATEC Computed | Error Type |
130- | ----------| ----------------| -----------------| ------------|
131- | Y_0(50) | -0.0560404718523358 | ** -0.0980649954700770** | ** 75% ERROR** |
132- | Y_0(100) | -0.0772433752531550 | -0.0772443133650831 | 0.01% error |
133-
134- ### Root Cause Analysis
135-
136- These failures likely stem from:
137- 1 . ** Range reduction errors** in trigonometric calculations for large arguments
138- 2 . ** Recurrence instability** for high-order Bessel functions
139- 3 . ** Asymptotic expansion limitations** not properly bounded
140-
141- ### Recommendations for Critical Applications
142-
143- 1 . ** DO NOT USE** J_n(x) for x > 100 without independent verification
144- 2 . ** DO NOT USE** Y_n(x) for x > 20 without independent verification
145- 3 . ** DO NOT USE** J_n(x) for n > 5 and x > 15 without verification
146- 4 . For large arguments, consider alternative implementations (e.g., AMOS library)
121+ ** Note:** An earlier version of this document contained incorrect reference values
122+ that were mistakenly attributed to NIST. The SLATEC implementations are accurate.
147123
148124---
149125
0 commit comments