Newlib and glibc both inherited a bunch of libm code from the venerable
SunPro C library, and it's still possible to diff the sources and see
how little they have diverged...
In trying to discover why newlib's bessel functions for 32-bit floats
are less accurate than glibc's, I compared the source code and found
some strange differences. While none of these explain the precision
issue, it seems like there are bugs here to fix.
(- is newlib, + is glibc)
(line 81, in j0f)
- if(ix>0x80000000) z = (invsqrtpi*cc)/__ieee754_sqrtf(x);
+ if(ix>0x48000000) z = (invsqrtpi*cc)/sqrtf(x);
(line 167 in y0f
- if(ix>0x80000000) z = (invsqrtpi*ss)/__ieee754_sqrtf(x);
+ if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x);
Line 81 and line 167 look like bugs to me -- 'ix' is the float with the
sign bit masked off, so those conditions can never be true.
openlibm has a different value that has a comment: