From owner-freebsd-numerics@FreeBSD.ORG Sun Jun 2 20:50:20 2013 Return-Path: Delivered-To: numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by hub.freebsd.org (Postfix) with ESMTP id 5D25BF13 for ; Sun, 2 Jun 2013 20:50:20 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail12.syd.optusnet.com.au (mail12.syd.optusnet.com.au [211.29.132.193]) by mx1.freebsd.org (Postfix) with ESMTP id 022A5189E for ; Sun, 2 Jun 2013 20:50:19 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail12.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id r52KoCi1006680 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 3 Jun 2013 06:50:18 +1000 Date: Mon, 3 Jun 2013 06:50:12 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: numerics@freebsd.org Subject: interesting bug in cexpf() caused by extra exponent range Message-ID: <20130603061718.M7113@besplex.bde.org> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.0 cv=e4Ne0tV/ c=1 sm=1 a=w8zoROYCmMcA:10 a=kj9zAlcOel0A:10 a=PO7r1zJSAAAA:8 a=JzwRw_2MAAAA:8 a=x7WX8rXFEGMA:10 a=mX8LXYcyzczS5YDSj0cA:9 a=CjuIK1q_8ugA:10 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sun, 02 Jun 2013 20:50:20 -0000 I found that cexpf(192.066 + I*1.176e-38) is broken with clang iff clang uses the i387 (the default without -march). Only the imaginary part is detected as being broken by my tests. This is because the imaginary part is evaluated as expf(192.066) * sinf(1.176e-38). expf(192.066) is expected to overflow, but it is evaluated as huge*huge = 1e60 in extra precision and exponent range. clang but not gcc manages to keep the extras. sinf(1.176e-38) is tiny and multiplication reduces the value to below FLT_MAX so that it no longer overflows to Inf when the extra exponent range is killed. I think that it is correct and good to evaluate expf(192.066) with extra precision and exponent range and then use these extras (except C11 breaks the returning of the extras). However, 1e60 is not the correct evaluation, which is 2.589e83. With the correct evaluation, the result of the multiplication (3.047e45) would exceed FLT_MAX in extra precision/range so it would overflow to Inf in float precision. I think the bug is in cexpf(). Its thresholds depend on overflow occuring, but without a strict assignment there is no gurantee that overflow occurs: @ if (hx >= exp_ovfl && hx <= cexp_ovfl) { @ /* @ * x is between 88.7 and 192, so we must scale to avoid @ * overflow in expf(x). @ */ @ return (__ldexp_cexpf(z, 0)); @ } else { @ /* @ * Cases covered here: @ * - x < exp_ovfl and exp(x) won't overflow (common case) @ * - x > cexp_ovfl, so exp(x) * s overflows for all s > 0 Actually, it doesn't overflow for all s > 0, because exp(x) is garbage instead of Inf as expected. @ * - x = +-Inf (generated by exp()) @ * - x = NaN (spurious inexact exception from y) @ */ @ exp_x = expf(x); The garbage would be converted to Inf if the compiler were perfectly pessimized to C90/99/(11?) spec. There is the usual bizarreness that if this used exp(x) to intentionally get extra precision/range, then the assignment would work and and lose the extras. Here the extra precision would be useful in the non-overflowing case, but the algorithm depends on overflow occuring. It would have to do a little more classification to force a value of Inf above cexp_ovfl if it used exp() and assigned it to a double or long double so as to keep the extra precision. @ return (cpackf(exp_x * cosf(y), exp_x * sinf(y))); @ } @ } I think a strict assignment is needed to fix this. Similarly in cexp(). This is fragile in user code too. Naive code might use expf(190) * 1e-38 (no assignment). Then overflow shouldn't occur, but the result is garbage. Bruce From owner-freebsd-numerics@FreeBSD.ORG Mon Jun 3 01:11:39 2013 Return-Path: Delivered-To: numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by hub.freebsd.org (Postfix) with ESMTP id B903AF7 for ; Mon, 3 Jun 2013 01:11:39 +0000 (UTC) (envelope-from das@freebsd.org) Received: from zim.MIT.EDU (50-196-151-174-static.hfc.comcastbusiness.net [50.196.151.174]) by mx1.freebsd.org (Postfix) with ESMTP id 862EC1F19 for ; Mon, 3 Jun 2013 01:11:38 +0000 (UTC) Received: from zim.MIT.EDU (localhost [127.0.0.1]) by zim.MIT.EDU (8.14.7/8.14.2) with ESMTP id r531BV4u098818; Sun, 2 Jun 2013 18:11:31 -0700 (PDT) (envelope-from das@freebsd.org) Received: (from das@localhost) by zim.MIT.EDU (8.14.7/8.14.2/Submit) id r531BURt098817; Sun, 2 Jun 2013 18:11:30 -0700 (PDT) (envelope-from das@freebsd.org) Date: Sun, 2 Jun 2013 18:11:30 -0700 From: David Schultz To: Bruce Evans Subject: Re: interesting bug in cexpf() caused by extra exponent range Message-ID: <20130603011130.GA84261@zim.MIT.EDU> References: <20130603061718.M7113@besplex.bde.org> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20130603061718.M7113@besplex.bde.org> Cc: numerics@freebsd.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 03 Jun 2013 01:11:39 -0000 On Mon, Jun 03, 2013, Bruce Evans wrote: > I found that cexpf(192.066 + I*1.176e-38) is broken with clang iff > clang uses the i387 (the default without -march). Only the imaginary > part is detected as being broken by my tests. > > This is because the imaginary part is evaluated as expf(192.066) * > sinf(1.176e-38). expf(192.066) is expected to overflow, but it > is evaluated as huge*huge = 1e60 in extra precision and exponent > range. clang but not gcc manages to keep the extras. > sinf(1.176e-38) is tiny and multiplication reduces the value to > below FLT_MAX so that it no longer overflows to Inf when the > extra exponent range is killed. > > I think that it is correct and good to evaluate expf(192.066) with > extra precision and exponent range and then use these extras (except > C11 breaks the returning of the extras). However, 1e60 is not the > correct evaluation, which is 2.589e83. With the correct evaluation, > the result of the multiplication (3.047e45) would exceed FLT_MAX > in extra precision/range so it would overflow to Inf in float > precision. > > I think the bug is in cexpf(). Its thresholds depend on overflow > occuring, but without a strict assignment there is no gurantee that > overflow occurs: I think you have identified two potential bugs, but neither are in cexpf(): one is in the compiler, and the other is in expf(). First the compiler bug: cexpf() is careful to assign the result of expf() to a float variable, and the C standard requires that the compiler not retain any extra precision. Second, the expf() bug: If the compiler decides to return the result in a register that has more than float precision, expf() needs to be prepared to deal with that. It can return INFINITY, or it can return the correct answer (which will overflow to infinity when the extra precision is dropped). Returning 1e60 is bogus. I think just using a macro, e.g., RETURN(huge * huge), to ensure that the extra precision gets dropped, is the most straightforward solution. > This is fragile in user code too. Naive code might use expf(190) * 1e-38 > (no assignment). Then overflow shouldn't occur, but the result is garbage. Yes, I think this is why it needs to be fixed in expf(), not in callers of expf(). From owner-freebsd-numerics@FreeBSD.ORG Mon Jun 3 11:06:49 2013 Return-Path: Delivered-To: freebsd-numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.FreeBSD.org [8.8.178.115]) by hub.freebsd.org (Postfix) with ESMTP id C8FE35DB for ; Mon, 3 Jun 2013 11:06:49 +0000 (UTC) (envelope-from owner-bugmaster@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:1900:2254:206c::16:87]) by mx1.freebsd.org (Postfix) with ESMTP id A190613BB for ; Mon, 3 Jun 2013 11:06:49 +0000 (UTC) Received: from freefall.freebsd.org (localhost [127.0.0.1]) by freefall.freebsd.org (8.14.7/8.14.7) with ESMTP id r53B6nra015109 for ; Mon, 3 Jun 2013 11:06:49 GMT (envelope-from owner-bugmaster@FreeBSD.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.7/8.14.7/Submit) id r53B6nlv015107 for freebsd-numerics@FreeBSD.org; Mon, 3 Jun 2013 11:06:49 GMT (envelope-from owner-bugmaster@FreeBSD.org) Date: Mon, 3 Jun 2013 11:06:49 GMT Message-Id: <201306031106.r53B6nlv015107@freefall.freebsd.org> X-Authentication-Warning: freefall.freebsd.org: gnats set sender to owner-bugmaster@FreeBSD.org using -f From: FreeBSD bugmaster To: freebsd-numerics@FreeBSD.org Subject: Current problem reports assigned to freebsd-numerics@FreeBSD.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 03 Jun 2013 11:06:49 -0000 Note: to view an individual PR, use: http://www.freebsd.org/cgi/query-pr.cgi?pr=(number). The following is a listing of current problems submitted by FreeBSD users. These represent problem reports covering all versions including experimental development code and obsolete releases. S Tracker Resp. Description -------------------------------------------------------------------------------- o stand/175811 numerics libstdc++ needs complex support in order use C99 o bin/170206 numerics [msun] [patch] complex arcsinh, log, etc. 2 problems total. From owner-freebsd-numerics@FreeBSD.ORG Wed Jun 5 09:00:32 2013 Return-Path: Delivered-To: numerics@freebsd.org Received: from mx1.freebsd.org (mx1.FreeBSD.org [8.8.178.115]) by hub.freebsd.org (Postfix) with ESMTP id 56A9CE61 for ; Wed, 5 Jun 2013 09:00:32 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail17.syd.optusnet.com.au (mail17.syd.optusnet.com.au [211.29.132.198]) by mx1.freebsd.org (Postfix) with ESMTP id BB2BD1CDE for ; Wed, 5 Jun 2013 09:00:31 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail17.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id r5590Lms018977 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Wed, 5 Jun 2013 19:00:23 +1000 Date: Wed, 5 Jun 2013 19:00:21 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: numerics@freebsd.org Subject: isinf() on NaN args Message-ID: <20130605183354.V1244@besplex.bde.org> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.0 cv=e4Ne0tV/ c=1 sm=1 a=kj9zAlcOel0A:10 a=PO7r1zJSAAAA:8 a=JzwRw_2MAAAA:8 a=Gbq2Tcle6C4A:10 a=BKSoS5k6mVHXEOtCdAMA:9 a=CjuIK1q_8ugA:10 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 05 Jun 2013 09:00:32 -0000 Is isinf(x) (or even fpclassify(x)) permitted to raise FE_INVALID if x is a (quiet) NaN? It shouldn't, but I can't find where C99 requires this. I can only find where C99 requires FE_INVALID to not be raised for comparison macros. It should raise if x is a signaling NaN, but C99 doesn't specify signaling NaNs and the hardware behaviour is more varied and the software bugs are larger for signaling NaNs, so getting this wrong in many cases doesn't matter much. Bruce From owner-freebsd-numerics@FreeBSD.ORG Wed Jun 5 16:06:01 2013 Return-Path: Delivered-To: numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by hub.freebsd.org (Postfix) with ESMTP id 36B71488 for ; Wed, 5 Jun 2013 16:06:01 +0000 (UTC) (envelope-from das@CSAIL.MIT.EDU) Received: from zim.MIT.EDU (50-196-151-174-static.hfc.comcastbusiness.net [50.196.151.174]) by mx1.freebsd.org (Postfix) with ESMTP id 1D6B610FC for ; Wed, 5 Jun 2013 16:06:00 +0000 (UTC) Received: from zim.MIT.EDU (localhost [127.0.0.1]) by zim.MIT.EDU (8.14.7/8.14.2) with ESMTP id r55G5xZ7044481; Wed, 5 Jun 2013 09:05:59 -0700 (PDT) (envelope-from das@CSAIL.MIT.EDU) Received: (from das@localhost) by zim.MIT.EDU (8.14.7/8.14.2/Submit) id r55G5wYK044480; Wed, 5 Jun 2013 09:05:58 -0700 (PDT) (envelope-from das@CSAIL.MIT.EDU) Date: Wed, 5 Jun 2013 09:05:58 -0700 From: David Schultz To: Bruce Evans Subject: Re: isinf() on NaN args Message-ID: <20130605160558.GA44189@zim.MIT.EDU> References: <20130605183354.V1244@besplex.bde.org> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20130605183354.V1244@besplex.bde.org> Cc: numerics@freebsd.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 05 Jun 2013 16:06:01 -0000 On Wed, Jun 05, 2013, Bruce Evans wrote: > Is isinf(x) (or even fpclassify(x)) permitted to raise FE_INVALID if > x is a (quiet) NaN? It shouldn't, but I can't find where C99 requires > this. I can only find where C99 requires FE_INVALID to not be raised > for comparison macros. > > It should raise if x is a signaling NaN, but C99 doesn't specify > signaling NaNs and the hardware behaviour is more varied and the > software bugs are larger for signaling NaNs, so getting this wrong in > many cases doesn't matter much. The equivalent operation in IEEE-754R never raises an exception, even for signaling NaNs. If we ignore signaling NaNs for a moment, I think the general idea is to raise an invalid exception in two situations: 1. at the moment when a NaN is produced from non-NaN inputs (this is similar to the rule that overflow is raised only when the overflow first occurs) 2. for operations that don't return a floating-point result, whenever there's no sensible answer to return The second case covers things like conversion to integer, ilogb, and comparison operators. One could argue that classification macros like isinf() always have a sensible true or false answer, though, because you're asking a question about what kind of value you have. From owner-freebsd-numerics@FreeBSD.ORG Wed Jun 5 16:38:38 2013 Return-Path: Delivered-To: numerics@freebsd.org Received: from mx1.freebsd.org (mx1.FreeBSD.org [8.8.178.115]) by hub.freebsd.org (Postfix) with ESMTP id C1D0BCDA for ; Wed, 5 Jun 2013 16:38:38 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.95.76.21]) by mx1.freebsd.org (Postfix) with ESMTP id A48D91255 for ; Wed, 5 Jun 2013 16:38:38 +0000 (UTC) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.6/8.14.6) with ESMTP id r55GbDlT074217; Wed, 5 Jun 2013 09:37:13 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.6/8.14.6/Submit) id r55GbDmR074216; Wed, 5 Jun 2013 09:37:13 -0700 (PDT) (envelope-from sgk) Date: Wed, 5 Jun 2013 09:37:13 -0700 From: Steve Kargl To: David Schultz Subject: Re: isinf() on NaN args Message-ID: <20130605163712.GA74110@troutmask.apl.washington.edu> References: <20130605183354.V1244@besplex.bde.org> <20130605160558.GA44189@zim.MIT.EDU> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20130605160558.GA44189@zim.MIT.EDU> User-Agent: Mutt/1.5.21 (2010-09-15) Cc: numerics@freebsd.org, Bruce Evans X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 05 Jun 2013 16:38:38 -0000 On Wed, Jun 05, 2013 at 09:05:58AM -0700, David Schultz wrote: > On Wed, Jun 05, 2013, Bruce Evans wrote: > > Is isinf(x) (or even fpclassify(x)) permitted to raise FE_INVALID if > > x is a (quiet) NaN? It shouldn't, but I can't find where C99 requires > > this. I can only find where C99 requires FE_INVALID to not be raised > > for comparison macros. > > > > It should raise if x is a signaling NaN, but C99 doesn't specify > > signaling NaNs and the hardware behaviour is more varied and the > > software bugs are larger for signaling NaNs, so getting this wrong in > > many cases doesn't matter much. > > The equivalent operation in IEEE-754R never raises an exception, > even for signaling NaNs. If we ignore signaling NaNs for a > moment, I think the general idea is to raise an invalid exception > in two situations: > > 1. at the moment when a NaN is produced from non-NaN inputs > (this is similar to the rule that overflow is raised only > when the overflow first occurs) > > 2. for operations that don't return a floating-point result, > whenever there's no sensible answer to return > > The second case covers things like conversion to integer, ilogb, > and comparison operators. One could argue that classification > macros like isinf() always have a sensible true or false answer, > though, because you're asking a question about what kind of value > you have. IIRC, C99 defers to IEC 60559 (aka IEEE 754) when a behavior is not specified within C99. The copy of IEEE 754 that I have states Some functions, such as the copy operation y := x without change of format, may at the implementor's option be treated as nonarithmetic operations which do not signal the invalid operation exception for signaling NaNs; the functions in question are (1), (2), (6), and (7). ... (7) Isnan(x), or equivalently x != x, returns the value TRUE if x is a NaN, and returns FALSE otherwise. This means that isnan is not required to raised the invalid operation. -- Steve From owner-freebsd-numerics@FreeBSD.ORG Wed Jun 5 19:06:43 2013 Return-Path: Delivered-To: numerics@FreeBSD.org Received: from mx1.freebsd.org (mx1.FreeBSD.org [8.8.178.115]) by hub.freebsd.org (Postfix) with ESMTP id BF051D63 for ; Wed, 5 Jun 2013 19:06:43 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail35.syd.optusnet.com.au (mail35.syd.optusnet.com.au [211.29.133.51]) by mx1.freebsd.org (Postfix) with ESMTP id 609BB1B7F for ; Wed, 5 Jun 2013 19:06:39 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail35.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id r55J5JDU021804 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Thu, 6 Jun 2013 05:05:21 +1000 Date: Thu, 6 Jun 2013 05:05:19 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl Subject: Re: isinf() on NaN args In-Reply-To: <20130605163712.GA74110@troutmask.apl.washington.edu> Message-ID: <20130606031050.Q17887@besplex.bde.org> References: <20130605183354.V1244@besplex.bde.org> <20130605160558.GA44189@zim.MIT.EDU> <20130605163712.GA74110@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.0 cv=e/de0tV/ c=1 sm=1 a=kj9zAlcOel0A:10 a=PO7r1zJSAAAA:8 a=JzwRw_2MAAAA:8 a=CgK77c7O4qAA:10 a=bMBPDBsYS8LiBcrmlpIA:9 a=CjuIK1q_8ugA:10 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 Cc: David Schultz , numerics@FreeBSD.org X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 05 Jun 2013 19:06:43 -0000 On Wed, 5 Jun 2013, Steve Kargl wrote: > On Wed, Jun 05, 2013 at 09:05:58AM -0700, David Schultz wrote: >> On Wed, Jun 05, 2013, Bruce Evans wrote: >>> Is isinf(x) (or even fpclassify(x)) permitted to raise FE_INVALID if >>> x is a (quiet) NaN? It shouldn't, but I can't find where C99 requires >>> this. I can only find where C99 requires FE_INVALID to not be raised >>> for comparison macros. >>> >>> It should raise if x is a signaling NaN, but C99 doesn't specify >>> signaling NaNs and the hardware behaviour is more varied and the >>> software bugs are larger for signaling NaNs, so getting this wrong in >>> many cases doesn't matter much. >> >> The equivalent operation in IEEE-754R never raises an exception, >> even for signaling NaNs. If we ignore signaling NaNs for a >> moment, I think the general idea is to raise an invalid exception >> in two situations: Is this new in 754R? 754 doesn't require it according to Steve's quote. I now remember this part of 754. > IIRC, C99 defers to IEC 60559 (aka IEEE 754) when a behavior is not > specified within C99. The copy of IEEE 754 that I have states > > Some functions, such as the copy operation y := x without change > of format, may at the implementor's option be treated as nonarithmetic > operations which do not signal the invalid operation exception for > signaling NaNs; the functions in question are (1), (2), (6), and (7). > ... > (7) Isnan(x), or equivalently x != x, returns the value TRUE if x is > a NaN, and returns FALSE otherwise. > > This means that isnan is not required to raised the invalid operation. Except for signaling NaNs. These are required to raise invalid for every operation in section 5. These include comparisons but don't include conversions without change of format. If isnan(x) really is required to be equivalent to x != x, IEEE754 isnan() is required to raise invalid for signaling NaNs. For quiet NaNs, isnan() raising invalid is implementation-defineded. Since C99 is older than 754R, it can't require much more than 754, so this must be implementation-defined in C99 too. isnan() on signaling NaNs is even more implementation-defined, since C99 doesn't specify signaling NaNs. All this seems to be implemented correctly by hardware and compilers on x86. Compilers alway use the unordered comparison operations fucom (i387) and ucom (SSE). These always raise invalid for signaling NaNs but never for quiet NaNs (I only checked this for fucom). Most of this is implemented less correctly for all macros at the level of isnan() except the comparison ones. Builtins are used for the latter, so they work. The other macros have various bugs, starting with their slowness. For signaling NaNs, their only (?) bug is that they only raise invalid accidentally (when there is a conversion in the parameter passing. This is very machine-dependent). This conforms to C99 but not to IEEE754. One thing I learned from investigating this recently is that it is not a bug for isnan() and friends to sometimes raise invalid for signaling NaNs. They accidentally conform to IEEE754 when they do this. x86 is easy to fix by using the builtins more. On other arches, it should be safe to use x != x for isnan() and fabs*(x) == INFINITY for isinf(), since even if these raise invalid incorrectly, they just get different cases wrong than now. These expressions can be used on x86 too, so no ifdefs are needed. Perhaps similarly for isfinite() (fabs*(x) < INFINITY?) and isnormal(). fpclassify() and signbit() are not so easy. I investigated this more when emaste@ told me that __isinfl() never worked on arches without real long doubles, except on alpha, since it uses a hard-coded MAX_EXP that is only ifdefed for alpha. Others like __isnanl() never worked on alpha either. Oops. I just remembered that __isinfl() and its long double friends are supposed to be unreacble on arches where they don't apply. The isinf(x) macro translates to isinf(x) for long double x on these arches. It would take either a compiler transalation of isinf(x) to __isinfl(x), or source code hacks that do the same thing, or a direct call to __isinfl() to reach it. There is no weak symbol that translates isinfl() to __isinfl() which would allow smaller code hacks to reach it. Bruce From owner-freebsd-numerics@FreeBSD.ORG Thu Jun 6 06:34:47 2013 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.FreeBSD.org [8.8.178.115]) by hub.freebsd.org (Postfix) with ESMTP id 301ED465 for ; Thu, 6 Jun 2013 06:34:47 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from zim.MIT.EDU (50-196-151-174-static.hfc.comcastbusiness.net [50.196.151.174]) by mx1.freebsd.org (Postfix) with ESMTP id 1A0CA1A61 for ; Thu, 6 Jun 2013 06:34:46 +0000 (UTC) Received: from zim.MIT.EDU (localhost [127.0.0.1]) by zim.MIT.EDU (8.14.7/8.14.2) with ESMTP id r566Yj5x046347 for ; Wed, 5 Jun 2013 23:34:45 -0700 (PDT) (envelope-from das@FreeBSD.ORG) Received: (from das@localhost) by zim.MIT.EDU (8.14.7/8.14.2/Submit) id r566Yj3Z046346 for freebsd-numerics@freebsd.org; Wed, 5 Jun 2013 23:34:45 -0700 (PDT) (envelope-from das@FreeBSD.ORG) Date: Wed, 5 Jun 2013 23:34:45 -0700 From: David Schultz To: freebsd-numerics@freebsd.org Subject: Numerics project webpage Message-ID: <20130606063445.GA46299@zim.MIT.EDU> MIME-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Thu, 06 Jun 2013 06:34:47 -0000 I put together a simple webpage to track the status of the numerics work, and as a repository for any future ideas, links to patches, and so forth: http://wiki.freebsd.org/Numerics As you can see from the table, there are a number of functions that are either trivial or pretty close to done. (Having expl and logl unblocked a lot of low-hanging fruit.) It would be great to get at least half of the remaining rows checked off before the 10.0 release. I put the page on the wiki so that others can edit it more easily. If you're a committer, please bug Simon for an account and help keep it up to date. Also, if anyone else is interested in helping, the page gives a pretty good idea of where things stand and who to contact about various tasks. From owner-freebsd-numerics@FreeBSD.ORG Fri Jun 7 09:16:42 2013 Return-Path: Delivered-To: numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:1900:2254:206a::19:1]) by hub.freebsd.org (Postfix) with ESMTP id B48B15BD for ; Fri, 7 Jun 2013 09:16:42 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail108.syd.optusnet.com.au (mail108.syd.optusnet.com.au [211.29.132.59]) by mx1.freebsd.org (Postfix) with ESMTP id 09AAF129B for ; Fri, 7 Jun 2013 09:16:42 +0000 (UTC) Received: from c122-106-156-23.carlnfd1.nsw.optusnet.com.au (c122-106-156-23.carlnfd1.nsw.optusnet.com.au [122.106.156.23]) by mail108.syd.optusnet.com.au (Postfix) with ESMTPS id 2B78F1A1AAC for ; Fri, 7 Jun 2013 19:16:40 +1000 (EST) Date: Fri, 7 Jun 2013 19:16:38 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: numerics@freebsd.org Subject: acoshl, asinhl, atanhl almost finished Message-ID: <20130607185037.F2756@besplex.bde.org> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Optus-CM-Score: 0 X-Optus-CM-Analysis: v=2.0 cv=Q6eKePKa c=1 sm=1 a=tYOxqLTWYZkA:10 a=kj9zAlcOel0A:10 a=PO7r1zJSAAAA:8 a=JzwRw_2MAAAA:8 a=pVaANWmWtz0A:10 a=i0XaXLn4f-YvG9u7eWMA:9 a=CjuIK1q_8ugA:10 a=ebeQFi2P/qHVC0Yw9JDJ4g==:117 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.14 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 07 Jun 2013 09:16:42 -0000 I converted the fdlibm double precision inverse hyperbolic functions to long double precision and am now testing them. atanhl has passed a lot of test on ld80: % --- e_atanh.c Thu May 30 18:14:16 2013 % +++ e_atanhl.c Fri Jun 7 18:07:32 2013 % @@ -1,2 +1,4 @@ % +/* from: FreeBSD: head/lib/msun/src/e_atanh.c 176451 2008-02-22 02:30:36Z das */ % +/* Conversion to float by Bruce D. Evans. */ % % /* @(#)e_atanh.c 1.3 95/01/18 */ % @@ -14,49 +16,44 @@ % % #include % -__FBSDID("$FreeBSD: src/lib/msun/src/e_atanh.c,v 1.9 2012/11/17 01:50:07 svnexp Exp $"); % +__FBSDID("$FreeBSD$"); % % -/* __ieee754_atanh(x) % - * Method : % - * 1.Reduced x to positive by atanh(-x) = -atanh(x) % - * 2.For x>=0.5 % - * 1 2x x % - * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) % - * 2 1 - x 1 - x % - * % - * For x<0.5 % - * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) % - * % - * Special cases: % - * atanh(x) is NaN if |x| > 1 with signal; % - * atanh(NaN) is that NaN with no signal; % - * atanh(+-1) is +-INF with signal. % - * % - */ % +#include % +#ifdef __i386__ % +#include % +#endif % % #include "math.h" % #include "math_private.h" % % -static const double one = 1.0, huge = 1e300; % +#if LDBL_MAX_EXP != 0x4000 % +/* We also require the usual expsign encoding. */ % +#error "Unsupported long double format" % +#endif % + % +#define BIAS (LDBL_MAX_EXP - 1) % + % +static const double one = 1.0; % +static volatile long double huge = 0x1p10000L; % static const double zero = 0.0; % % -double % -__ieee754_atanh(double x) % +long double % +atanhl(long double x) % { This was supposed to look as much like the double version as the float version does, including style bugs, but almost every line from here on needed changing so I fixed most of the horizontal whitespace. % - double t; % - int32_t hx,ix; % - u_int32_t lx; % - EXTRACT_WORDS(hx,lx,x); % - ix = hx&0x7fffffff; % - if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */ % - return (x-x)/(x-x); Extracting all the bits for this test is too painful, and this test is inefficient anyway, so use an FP comparison for classifying +-1. Only classificatins based on the exponent alone are easy to use and efficient, especially in long double prcision. % - if(ix==0x3ff00000) % - return x/zero; % - if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */ This threshold is badly chosen (and became worse when it was blindly copied for float precision). Change to a higher threshold corresponding to 2**-26 here. x is not a very good approximation to atanh(x) when x is in the changed range, but it is much better than the general approximation used for the rest of the range <= 0.5 (the latter has roundoff errors of 1.0-1.5 ulps over the whole range). % - SET_HIGH_WORD(x,ix); % - if(ix<0x3fe00000) { /* x < 0.5 */ % + long double t; % + uint16_t hx, ix; % + % + ENTERI(); % + GET_LDBL_EXPSIGN(hx, x); % + ix = hx & 0x7fff; % + if (ix >= 0x3fff) /* |x| >= 1 or x is NaN or misnormal */ % + RETURNI(fabsl(x) == 1 ? x / zero : (x - x) / (x - x)); The RETURNI()s give larger diffs. Everything else changed literally on these 2 lines (previously 4 more complicated lines), so change the logic a little too. % + if (ix < BIAS - (LDBL_MANT_DIG + 1) / 2 && (huge + x) > zero) % + RETURNI(x); /* |x| < 2**-(MANT_DIG+1)/2 */ Misnormals are handled in all cases by arranging to do arithmetic on them so that the get converted to a qNaN for return. For example, if x is unnormal here, (huge + x) > zero is false here (x acts like an sNaN), so we don't just return it here, but fall through to cases where it gets converted. % + SET_LDBL_EXPSIGN(x, ix); % + if (ix < 0x3ffe) { /* |x| < 0.5 or x is misnormal */ % t = x+x; % - t = 0.5*log1p(t+t*x/(one-x)); % + t = 0.5*log1pl(t+t*x/(one-x)); % } else Preserving the style bugs includes not adding a comment for this case. % - t = 0.5*log1p((x+x)/(one-x)); % - if(hx>=0) return t; else return -t; % + t = 0.5*log1pl((x+x)/(one-x)); % + RETURNI((hx & 0x8000) == 0 ? t : -t); % } As expected since the algorithm wasn't changed and there few possibilities for numerical accidents, the errors in ulps seem to be about the same as in lower precisions (actually almost half an ulp better because log1pl() is better). Bruce