Test case example: // gcc main.c -std=c99 -lm <-m32> #include <stdio.h> #include <fenv.h> #include <stdlib.h> // abs function for int argument extern float j0f(float); extern double j0(double); #define i32(a) (((int*)&(a))[0]) int main() { float a,r1,r2; i32(a)=0x413caa1f; //argument fesetround(FE_TONEAREST); r1 = j0f(a); //actual result r2 = j0(a); //actual result 2 printf("inputs: a = (0x%08x)\t%e\n", i32(a), a); printf("actual = (0x%08x)\t%.13e\n", i32(r1), r1); printf("expected = (0x%08x)\t%.13e\n", i32(r2), r2); printf("error = %d ulp\n", abs(i32(r1) - i32(r2))); return 0; } Results x86_32: inputs: a = (0x413caa1f) 1.179153e+01 actual = (0xb471e828) -2.2529331999976e-07 expected = (0xb471d452) -2.2522115727952e-07 error = 5078 ulp Results x86_64: inputs: a = (0x413caa1f) 1.179153e+01 actual = (0xb471d9d6) -2.2524122300638e-07 expected = (0xb471d452) -2.2522115727952e-07 error = 1412 ulp
*** Bug 14468 has been marked as a duplicate of this bug. ***
For reference for anyone working on improving accuracy for Bessel functions (an issue for all the versions for all floating-point formats, not just x86 / x86_64): http://www.cl.cam.ac.uk/~jrh13/papers/bessel.ps.gz discusses how to achieve correct rounding for Bessel functions of any fixed constant order (jn/yn will be harder).
>> please give testcases that use hex floats New test case for J0f Bessel function. #include <stdio.h> #include <math.h> #include <fenv.h> extern float j0f(float); extern double j0(double); int main() { float a,r1; double r2; a=0x1.79543ep3; //argument fesetround(FE_TONEAREST); r1 = j0f(a); //actual result r2 = j0(a); //actual result 2 printf("inputs: a = %.6a\n", a); printf("actual = %.6a\n", r1); printf("expected = %.6a\n", r2); printf("error = %f ulp\n", fabs(r1 - r2) * exp2(23.0-logbf(r1))); return 0; } 64 bit result: inputs: a = 0x1.79543ep+3 actual = -0x1.e3b3acp-23 expected = -0x1.e3a8a3p-23 error = 1412.411834 ulp 32 bit result: inputs: a = 0x1.79543ep+3 actual = -0x1.e3d050p-23 expected = -0x1.e3a8a3p-23 error = 5078.411826 ulp
here is the worst case found by exhaustive search on all binary32 values with glib 2.31 on x86_64: #include <stdio.h> #include <math.h> #include <fenv.h> extern float j0f(float); extern double j0(double); int main() { float a,r1; double r2; a=0x1.4665d2p24; fesetround(FE_TONEAREST); r1 = j0f(a); //actual result r2 = j0(a); //actual result 2 printf("inputs: a = %.6a\n", a); printf("actual = %.6a\n", r1); printf("expected = %.6a\n", r2); printf("error = %f ulp\n", fabs(r1 - r2) * exp2(23.0-logbf(r2))); return 0; } inputs: a = 0x1.4665d2p+24 actual = 0x1.1c6742p-40 expected = 0x1.4a040dp-49 error = 4760682489.747789 ulp
the large error for huge operands can be avoided by removing the optimization "if(ix>0x48000000)" around line 63 of sysdeps/ieee754/flt-32/e_j0f.c. Then for x=0x1.4665d2p24 (x=21390802) we get an error of 239 ulps only (which corresponds to the cancellation between u*cc and v*ss in the expression u*cc-v*ss at line 66).
To fix the large errors properly, all of the j0 / j1 / y0 / y1 functions probably need to be reimplemented along the lines described in <http://www.cl.cam.ac.uk/~jrh13/papers/bessel.ps.gz> - expansions around each zero of the functions for small arguments, expansions based on phase and magnitude or similar functions for larger arguments (plus range reduction for phase) where the phase is either rigorously / exhaustively shown to be accurate enough (for float) or is heuristically accurate enough (for other floating-point formats, if exhaustive checks around the zeroes, until the phase approximation gets simple enough to use results on rational approximation of pi, are infeasible).
<https://sourceware.org/ml/libc-alpha/2013-03/msg00345.html> has my comments on these thresholds in various Bessel function implementations (note that in some cases they also serve to avoid spurious underflow / overflow exceptions, so removing such checks in one place would require adding them in another).
with the change proposed in comment #5, the maximal ulp error drops from 4760682496 to 6177902, and the number of binary32 inputs which give an error of 2 ulps or more drops from 310998452 to 269351612 (13% less).
I would expect removing that check to result in spurious overflow and underflow exceptions where pzerof and qzerof compute 1/(x*x), as they don't have any other checks to avoid such exceptions. Replacing the 0x48000000 checks for float by something like 0x5c000000 should be heuristically safe (although the functions could still do with rewrites along the lines described by Harrison).
The master branch has been updated by Joseph Myers <jsm28@sourceware.org>: https://sourceware.org/git/gitweb.cgi?p=glibc.git;h=ad180676b83dc1782d407dbff57dabbaab0c1f71 commit ad180676b83dc1782d407dbff57dabbaab0c1f71 Author: Joseph Myers <joseph@codesourcery.com> Date: Fri Feb 14 14:16:25 2020 +0000 Adjust thresholds in Bessel function implementations (bug 14469). A recent discussion in bug 14469 notes that a threshold in float Bessel function implementations, used to determine when to use a simpler implementation approach, results in substantially inaccurate results. As I discussed in <https://sourceware.org/ml/libc-alpha/2013-03/msg00345.html>, a heuristic argument suggests 2^(S+P) as the right order of magnitude for a suitable threshold, where S is the number of significand bits in the floating-point type and P is the number of significant bits in the representation of the floating-point type, and the float and ldbl-96 implementations use thresholds that are too small. Some threshold does need using, there or elsewhere in the implementation, to avoid spurious underflow and overflow for large arguments. This patch sets the thresholds in the affected implementations to more heuristically justifiable values. Results will still be inaccurate close to zeroes of the functions (thus this patch does *not* fix any of the bugs for Bessel function inaccuracy); fixing that would require a different implementation approach, likely along the lines described in <http://www.cl.cam.ac.uk/~jrh13/papers/bessel.ps.gz>. So the justification for a change such as this would be statistical rather than based on particular tests that had excessive errors and no longer do so (no doubt such tests could be found, but would probably be too fragile to add to the testsuite, as liable to give large errors again from very small implementation changes or even from compiler changes). See <https://sourceware.org/ml/libc-alpha/2020-02/msg00638.html> for such statistics of the resulting improvements for float functions. Tested (glibc testsuite) for x86_64.
with glibc-2.32, the maximum error is now 6177902 ulps on x86_64 (for x=3.153646966e+38)
fixed by commit 9acda61