Unreachable code in ef_j0.c differs from glibc

[keithp@keithp.com ("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 --]