From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:02:54 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id 9F3E9106564A for ; Sun, 12 Aug 2012 23:02:54 +0000 (UTC) (envelope-from peter@rulingia.com) Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by mx1.freebsd.org (Postfix) with ESMTP id 24D778FC0A for ; Sun, 12 Aug 2012 23:02:53 +0000 (UTC) Received: from server.rulingia.com (c220-239-249-137.belrs5.nsw.optusnet.com.au [220.239.249.137]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CN2rbb075616 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:02:54 +1000 (EST) (envelope-from peter@rulingia.com) X-Bogosity: Ham, spamicity=0.000000 Received: from server.rulingia.com (localhost.rulingia.com [127.0.0.1]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CN2lhO021143 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:02:47 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CN2lsU021142 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:02:47 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:02:47 +1000 Resent-Message-ID: <20120812230247.GD20453@server.rulingia.com> Resent-To: freebsd-numerics@freebsd.org Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q6O3I8IZ065193 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Tue, 24 Jul 2012 13:18:08 +1000 (EST) (envelope-from brde@optusnet.com.au) Received: from mail04.syd.optusnet.com.au (mail04.syd.optusnet.com.au [211.29.132.185]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q6O3I8K2017903 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Tue, 24 Jul 2012 13:18:08 +1000 (EST) (envelope-from brde@optusnet.com.au) Received: from c122-106-171-246.carlnfd1.nsw.optusnet.com.au (c122-106-171-246.carlnfd1.nsw.optusnet.com.au [122.106.171.246]) by mail04.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q6O3HqKj021486 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Tue, 24 Jul 2012 13:17:53 +1000 From: Bruce Evans Mail-Followup-To: freebsd-numerics@freebsd.org X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <500DAD41.5030104@missouri.edu> Message-ID: <20120724113214.G934@besplex.bde.org> References: <20120717084457.U3890@besplex.bde.org> <5004A5C7.1040405@missouri.edu> <5004DEA9.1050001@missouri.edu> <20120717040118.GA86840@troutmask.apl.washington.edu> <20120717042125.GF66913@server.rulingia.com> <20120717043848.GB87001@troutmask.apl.washington.edu> <20120717225328.GA86902@server.rulingia.com> <20120717232740.GA95026@troutmask.apl.washington.edu> <20120718001337.GA87817@server.rulingia.com> <20120718123627.D1575@besplex.bde.org> <20120722121219.GC73662@server.rulingia.com> <500DAD41.5030104@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Mailman-Approved-At: Sun, 12 Aug 2012 23:55:59 +0000 Cc: Diane Bruce , Bruce Evans , John Baldwin , David Chisnall , Bruce Evans , Steve Kargl , David Schultz , Peter Jeremy , Warner Losh Subject: Re: Use of C99 extra long double math functions after r236148 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Date: Sun, 12 Aug 2012 23:02:54 -0000 X-Original-Date: Tue, 24 Jul 2012 13:17:52 +1000 (EST) X-List-Received-Date: Sun, 12 Aug 2012 23:02:54 -0000 On Mon, 23 Jul 2012, Stephen Montgomery-Smith wrote: > I just realized that catan(z) = reverse( catanh(reverse (z))), just like > casin relates to casinh (remember reverse(x+I*y) = y+I*x). This is a > consequence of catan and catanh being odd functions, as well as the standard > relation catan(z) = -I*catanh(I*z). > > So I would modify Peter's code by taking out the minus signs. > > Maybe it would make a difference if the answers involved -0. Quite likely. Note that keeping the real and imaginary parts separate is giving perfect reflection across some axes. Reflecting 0 is for other axes and it doesn't fall out so accidentally. > I feel that I am done with these functions for now. I tried to change my > comments to conform to the style given to me by Bruce. However spacing > inside mathematical expressions is something where I am inconsistent. Run it through indent to find all the inconsistencies. > The functions still need a lot of work to handle -0, infs and NaNs correctly. > I will leave that to you guys, because you seem so much better at it than me. > I still don't understand why the proper test is "if (x!=x) return(x+x)" > rather than "if (isnan(x)) return(NAN)". First, signaling NaNs are required to signal, in a strong way: (IEEE doesn't specify complex functions of course, but it requires this for all the functions and most of the operations that it specifies. The main exceptions (in IEEE-854-1987) are: - copying a signaling NaN without changing its precision. Whether this triggers is implementation-defined. So i387 conforms by triggering for changing the precision on loads iff the load does a non-null change to long double precision - the sign of a NaN is not specified. Triggering for signaling NaNs is still required "for every operation listed in Section 5". These operations don't seem to include absolute values or negation. Thus the i387 fabs and fchs conform because they are extensions, but they should never be used!?, since any use doesn't conform to the spirit of the standard.) Every "operation" involving a signaling NaN or invalid operation shall, if no "trap" occurs and a floating point result is to be delivered, deliver "a" quiet NaN as its result. "if (isnan(x)) return(NAN)" does extra work to fail to conform to this. First it uses the slow isnan() classification to be sure that no "trap" occurs. (Whether (x != x) triggers signaling NaNs, and whether the triggering involves a trap or just setting FE_INVALID is delicate. This was broken in the original i8087. Actually, on checking the details, I found that it was for quiet NaNs that was broken. (x != x) should and does set FE_INVALID for signaling NaNs on all x87. This traps iff it is unmasked. Comparison doesn't deliver an FP result, so quietening of x is irrelevant. So (x != x) works right for signaling NaNs on all x87. The i8087 bug was that its comparison operators all set FE_INVALID, etc for quiet NaNs too. Later x87 have comparison operators that work right on all types of NaNs, and compilers use them. Thus you can write isnan() as (x != x) on at least x87 without getting a spurious FE_INVALID for quiet NaNs. I think this is unportable, so isnan() should be used if it is necessary to avoid such FE_INVALID. But where we want FE_INVALID, and we don't care if we get it from both the comparison and (x+x). Second, returning x wouldn't conform, since the returned x might not be quiet if the original x isn't. Returning the `NAN' does conform, provided NaN is quiet, since the standard allows "a" (that is, any quiet NaN to be returned). BTW, C has very bad names for NaNs and Infs. `NaN' should be spelled with its 'a' in lower case, as in the standard. Here the NaN returned should be spelled qNaN or as my dNaN (default NaN) to emphasize that it is a quiet default NaN. The standard then recommends which NaN "should" be returned, and it isn't the C default NaN: Every operation [...] shall [...] if a floating point result is to be delivered, [when this result is specified to be "a" quiet NaN,] then this NaN _should_ be _one_ of its input NaNs. Note that precision conversions might be unable to deliver the same NaN. This is only a "should", but returning the default C NaN fails to conform to at least its spirit, by making no attempt to return the original NaN. Returning the an original one is especially simple when there is only one arg. This arg just needs to be quietened and undergo any necessary conversions. (x+x) does this (modulo excessive conversions due to C's promotion rules). When the arg is a signaling NaN, it is impossible to return exactly it, so returning an unrelated default would be more reasonable, but doing this would take more code. It is best to depend on the hardware doing the right thing for (x+x) on signaling NaN x. The hardware is specified to give a quiet result and would need much the same complications as us return a default instead of x slightly modified for signaling NaNs only. This also specifies NaN mixing: when there are _two_ args, _one_ of them _should_ be returned. Unfortunately, which one is not specified. But hardware operations should pick one, and we can do about as well as the hardware by mixing (= choosing) using (x+y). Bruce