From mboxrd@z Thu Jan 1 00:00:00 1970 Return-Path: Received: from esa4.mentor.iphmx.com (esa4.mentor.iphmx.com [68.232.137.252]) by sourceware.org (Postfix) with ESMTPS id 3019E3857C50 for ; Thu, 10 Sep 2020 15:27:02 +0000 (GMT) DMARC-Filter: OpenDMARC Filter v1.3.2 sourceware.org 3019E3857C50 Authentication-Results: sourceware.org; dmarc=none (p=none dis=none) header.from=codesourcery.com Authentication-Results: sourceware.org; spf=pass smtp.mailfrom=joseph_myers@mentor.com IronPort-SDR: 3pDsp2PnbSl1OwCHyMPktb877TTHSDQ+VWpQ3YInqm61rPl4xbCTg6M3+XcFG+ypwOIBnIkC1Y fUMEnFXwpWAL0dE5/IaeZceyQk6cwfuvMq4L9zsO9CeC0nQFxrigv+prriNVtAvA4ccFfI4y2V tr8gkpxNJHr9ajV5Hcv4AOphEDGQXdl+Xz/OEXR8zgUXHjwFoK5ujQNgXfal6/Ug3XfC+B/y71 swvy9vG1/YIIRrvMUnJThRMnsaOWZV3hxmyX8xmAKPfdTljUCT+cZiBT249G8/+WkotRpHZMxI 7mA= X-IronPort-AV: E=Sophos;i="5.76,413,1592899200"; d="scan'208";a="52894885" Received: from orw-gwy-02-in.mentorg.com ([192.94.38.167]) by esa4.mentor.iphmx.com with ESMTP; 10 Sep 2020 07:27:00 -0800 IronPort-SDR: 8AaNKvKN+RAn1CxyPN/XhVDJCBAV3amJAHI8PqjSTt5+0kB1IvuzbsIUs769xVefegjK463oRX IMvcZOPBUmH0siPV2H9wBIRT/8lYvZHa8MJSyb7WUhedu1qm2QB6UxgKEIvBbaMaaCnKWXdRKq ge8p2/j/gU1txnU3Mvh4tM9iPVWVcaGkbKq3+qTQKX9EM+2sNb4xK3ZiXLKcSDTdB5n1E3xfik JhphNAc95pH6JyqXkCgy2pKK3VmYRFSskOuUQ5mNwp+00FiVYC6XQdzM66gi11Pev151njmSiY Bcc= Date: Thu, 10 Sep 2020 15:26:55 +0000 From: Joseph Myers X-X-Sender: jsm28@digraph.polyomino.org.uk To: Ruinland ChuanTzu Tsai CC: Subject: Re: In libm, sin(qNaN) doesn't expect FE_INVALID ? In-Reply-To: <20200910134901.GA14099@APC301.andestech.com> Message-ID: References: <20200903123442.GA5266@APC301.andestech.com> <20200908110214.GA3365@APC301.andestech.com> <20200910134901.GA14099@APC301.andestech.com> User-Agent: Alpine 2.21 (DEB 202 2017-01-01) MIME-Version: 1.0 X-Originating-IP: [137.202.0.90] X-ClientProxiedBy: svr-ies-mbx-02.mgc.mentorg.com (139.181.222.2) To svr-ies-mbx-01.mgc.mentorg.com (139.181.222.1) X-Spam-Status: No, score=-3127.2 required=5.0 tests=BAYES_00, HEADER_FROM_DIFFERENT_DOMAINS, KAM_DMARC_STATUS, SPF_HELO_PASS, SPF_PASS, TXREP autolearn=no autolearn_force=no version=3.4.2 X-Spam-Checker-Version: SpamAssassin 3.4.2 (2018-09-13) on server2.sourceware.org Content-Type: text/plain; charset="utf-8" Content-Transfer-Encoding: 8BIT X-Content-Filtered-By: Mailman/MimeDel 2.1.29 X-BeenThere: libc-alpha@sourceware.org X-Mailman-Version: 2.1.29 Precedence: list List-Id: Libc-alpha mailing list List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Thu, 10 Sep 2020 15:27:04 -0000 On Thu, 10 Sep 2020, Ruinland ChuanTzu Tsai wrote: > I'm a little bit confused by the implementation of ULPDIFF() inside > `math/libm-test-support.c` which is : > > ``` > #define ULPDIFF(given, expected) \ (FUNC(fabs) ((given) - (expected)) / ulp (expected) > ``` > > and it looks _not_ really the same as the formula inside glibc's docu- > mentation [1] : > ` |d.d...d - (z / 2^e)| / 2^(p - 1) ` > > ( For a number z with the representation d.d…d·2^e and p is the number > of bits in the mantissa of the floating-point number representation. ) > > The denominator part of these two seems to have different meaning ? It looks like that formula from the manual should actually be multiplying by 2^(p-1), not dividing, to get an actual figure in ulps. Note that there are at least two different measures of errors in ulps. The one used in glibc is that we take an ideal correctly rounded result, take the absolute value of the difference between that and the result returned by the function, and divide that by a unit in the last place of the correctly rounded result. This gives an error that is almost always an integer number of ulps (it can be a non-integer if the result returned has a lower exponent than the correctly rounded result). A correctly rounded result has a 0 ulps error by this definition (but that's not sufficient for being correctly rounded; correct rounding also requires the correct sign of 0 and correct exceptions). Another version sometimes seen in the literature defines ulps not for a correctly rounded result but for the infinite-precision mathematical result. When that's given as "the absolute value of the difference between the two floating-point numbers closest to x, one of which may equal x", note that if x is a power of 2 (and exceeds the magnitude of the least normal value), or rounds away from zero to a power of 2, then this gives a definition of ulp that's half the one used by glibc (and thus an error that's twice that of the glibc definition). Then the error in a function is determined by comparing the rounded value to the infinite-precision value, in terms of ulps of the infinite-precision value. With this definition, a correctly rounded result has error at most 0.5 ulps in round-to-nearest mode and less than 1 ulp in other modes (but again, that's not sufficient for being correctly rounded). > Besides this issue, I would like to know that is there any written > policy for loosening or tightening the ULPs for mathematic functions ? Only the functions bound to IEEE operations (sqrt, fma, etc.) are expected to be correctly rounded. For others, people have typically found performance can be improved without introducing large errors. My guess is that most functions could be made to achieve 1ulp errors in round-to-nearest and 2ulp in other modes (whichever definition is used) without making performance worse. > And if someone is introducing a new platform to glibc, are there any > rules to regulate ? e.g. "ccosh" mustn't have a ulp more than ...... The general rule for new platforms is to avoid having architecture-specific function implementations that aren't actually needed, and to improve performance by improving the generic C implementations instead; see . Architecture-specific versions of functions such as fma that are fully bound to IEEE operations may make sense, where there are relevant hardware instructions. Architecture-specific versions of transcendental functions are almost surely a bad idea. Once you're using the architecture-independent implementations, you should have the same ulps as for most other platforms (modulo minor differences arising from compiler choices in whether to contract operations, if one of the platforms has fused multiply-add instructions). The only architecture-specific implementation of ccosh is for m68k (the alpha version is just dealing with compatibility for past ABI changes). The m68k version really ought to go away because of its use of fsincos (see bug 13742 regarding use of fsincos on m68k, and note that emulators may well not accurately reflect hardware inaccuracy there). -- Joseph S. Myers joseph@codesourcery.com