From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from elaine.keithp.com (home.keithp.com [63.227.221.253]) by sourceware.org (Postfix) with ESMTPS id F280F38A1403 for ; Tue, 17 Mar 2020 20:01:32 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org F280F38A1403 Received: from localhost (localhost [127.0.0.1]) by elaine.keithp.com (Postfix) with ESMTP id 9FD573F2B53E for ; Tue, 17 Mar 2020 13:01:31 -0700 (PDT) X-Virus-Scanned: Debian amavisd-new at keithp.com Received: from elaine.keithp.com ([127.0.0.1]) by localhost (elaine.keithp.com [127.0.0.1]) (amavisd-new, port 10024) with LMTP id fZ5LqZlxqFG9 for ; Tue, 17 Mar 2020 13:01:31 -0700 (PDT) Received: from keithp.com (koto.keithp.com [10.0.0.2]) by elaine.keithp.com (Postfix) with ESMTPSA id 3C2813F2B53D for ; Tue, 17 Mar 2020 13:01:31 -0700 (PDT) Received: by keithp.com (Postfix, from userid 1000) id E93791582127; Tue, 17 Mar 2020 13:01:30 -0700 (PDT) From: keithp@keithp.com ("Keith Packard") To: newlib@sourceware.org Subject: Unreachable code in ef_j0.c differs from glibc Date: Tue, 17 Mar 2020 13:01:30 -0700 Message-ID: <874kumn42t.fsf@keithp.com> MIME-Version: 1.0 Content-Type: multipart/signed; boundary="=-=-="; micalg=pgp-sha256; protocol="application/pgp-signature" X-Spam-Status: No, score=-0.2 required=5.0 tests=DKIM_SIGNED, DKIM_VALID, DKIM_VALID_AU, DKIM_VALID_EF, RCVD_IN_DNSWL_NONE, SPF_HELO_NONE, SPF_PASS autolearn=ham autolearn_force=no version=3.4.2 X-Spam-Checker-Version: SpamAssassin 3.4.2 (2018-09-13) on server2.sourceware.org X-BeenThere: newlib@sourceware.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Newlib mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 17 Mar 2020 20:01:34 -0000 --=-=-= Content-Type: text/plain 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 --=-=-= Content-Type: application/pgp-signature; name="signature.asc" -----BEGIN PGP SIGNATURE----- iQIzBAEBCAAdFiEEw4O3eCVWE9/bQJ2R2yIaaQAAABEFAl5xLJoACgkQ2yIaaQAA ABFnMQ/8D25fpf9iXkOK+V8mTy+Ni3TsH/pxejKW9yjFkEpeUgizwW11BPPy/vGu UvLvy/HMe4XiosLuG/7nweCzGJ/9QZCmtNG0BFpNWqLu6QDbXNAFatJ2LFP8z5RM GyeZ+O4AJak52d+y6ZCTmkUdjXr7W2qv4qwkq14NxgAMG8XNzWhX7NzOltidSiVR CEBvJULyc5A5AZXzKHJGSBoCVtuJAOgtWBPshnfQZ5M7Bnb2x2F2j+DcidnQ4sFX ppXbkXpp0jbv3uBofIuWCwLMAAEbkPshA9uuhVtjWpqk5AR1ZlYfuIcz3epw3vYa 8cdP1JjSCJUbGPQSsdFW93S7rJ0yXfDh6GMJFyg5WkuUhkxgxnt/9SAb4ubBDq22 yp7bFMdwZgRSvwddguzKe8Vx6LNtQiai/FDpB52Yl3I2Fp+3Z9+Ydfd3d1WCRn9B KvBzDPZVnMGryImpooK1WL4/DHiGFdtvNkEgE7+vf8KVmMz8fUjzf9hMHtEyx2/K QzMT5mTJ1Cyk4DM+nRz+qBF7phNdwlMKLGnhtYPt5V0ABbmZrg7dCKwiAlR3Ma+p Iqx6yd0x9xAB7YFRyiQ99HakILdt6WBpOCt5cEgyeo+g221d/BW84Aq8quh/YlAK 5/fSO2BsbzexBlLxhbUpSVkF2+lUVq7XrDUDX1ww2Paz8bhSoEg= =0EMk -----END PGP SIGNATURE----- --=-=-=--