Unreachable code in ef_j0.c differs from glibc

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

Unreachable code in ef_j0.c differs from glibc

Sourceware - newlib list mailing list

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:

        https://github.com/JuliaMath/openlibm/blob/master/src/e_j0f.c:

        if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */

Apple has the same value as newlib:

        https://opensource.apple.com/source/Libm/Libm-47.1/i386.subproj/e_j0f.c.auto.html

        if((unsigned long)ix>0x80000000ul) z = (invsqrtpi*cc)/sqrtf(x);

An older version from freebsd has the newlib value as well:

        web.mit.edu/freebsd/head/lib/msun/src/e_j0f.c

        if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);

(line 174 in y0f)
- if(ix<=0x32000000) { /* x < 2**-27 */
+ if(ix<=0x39800000) { /* x < 2**-13 */
     return(u00 + tpi*__ieee754_logf(x));
  }
  z = x*x;

Line 174 just seems like a missed optimization case; y0f for values <
2**-13 can be approximated with the shorter computation.

As to the accuracy issues, I haven't found the root cause yet.

--
-keith

signature.asc (847 bytes) Download Attachment
Reply | Threaded
Open this post in threaded view
|

Re: Unreachable code in ef_j0.c differs from glibc

Joseph Myers
On Tue, 17 Mar 2020, "Keith Packard" via Newlib wrote:

> (- 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);

That's 0x5c000000 in glibc after commit
ad180676b83dc1782d407dbff57dabbaab0c1f71.

--
Joseph S. Myers
[hidden email]