Unreachable code in ef_j0.c differs from glibc

Unreachable code in ef_j0.c differs from glibc

["Keith Packard"]

[-- Attachment #1: Type: text/plain, Size: 1786 bytes --]


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

[-- Attachment #2: signature.asc --]
[-- Type: application/pgp-signature, Size: 832 bytes --]

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
joseph@codesourcery.com