From owner-freebsd-numerics@FreeBSD.ORG Mon Jul 30 00:21:06 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id EAC84106564A for ; Mon, 30 Jul 2012 00:21:06 +0000 (UTC) (envelope-from stephen@missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by mx1.freebsd.org (Postfix) with ESMTP id BC5908FC0C for ; Mon, 30 Jul 2012 00:21:06 +0000 (UTC) Received: from [127.0.0.1] (wilberforce.math.missouri.edu [128.206.184.213]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q6U0L5cD051097 for ; Sun, 29 Jul 2012 19:21:06 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <5015D372.2060502@missouri.edu> Date: Sun, 29 Jul 2012 19:21:06 -0500 From: Stephen Montgomery-Smith User-Agent: Mozilla/5.0 (X11; Linux i686; rv:14.0) Gecko/20120714 Thunderbird/14.0 MIME-Version: 1.0 To: freebsd-numerics@freebsd.org References: <500DAD41.5030104@missouri.edu> <20120724113214.G934@besplex.bde.org> <501204AD.30605@missouri.edu> <20120727032611.GB25690@server.rulingia.com> <20120728125824.GA26553@server.rulingia.com> <501460BB.30806@missouri.edu> <20120728231300.GA20741@server.rulingia.com> <50148F02.4020104@missouri.edu> <20120729222706.GA29048@server.rulingia.com> <5015BB9F.90807@missouri.edu> <20120729235302.GA77192@server.rulingia.com> In-Reply-To: <20120729235302.GA77192@server.rulingia.com> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit 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: , X-List-Received-Date: Mon, 30 Jul 2012 00:21:07 -0000 On 07/29/2012 06:53 PM, Peter Jeremy wrote: > [Pruning CC list to keep mailman happy]] I had mailman complain to me in my last post. But it only took the moderator a few minutes to let it through. > > On 2012-Jul-29 17:39:27 -0500, Stephen Montgomery-Smith wrote: >> On 07/29/2012 05:27 PM, Peter Jeremy wrote: >>> WG14/N1256 G.6. I hadn't considered extending that to verifying that >>> purely real or imaginary inputs give purely real or imaginary outputs, >>> with the appropriately signed zero. This might be reasonable but it's >>> not completely trivial to implement in general since the domains of >>> the real part can be different. >> >> Maybe this should be a different program, since its logical structure >> would be quite different. In particular, you wouldn't be checking the >> value of the non-zero parts. > > Adding code to skip checks on the real or imaginary part of the > result is quite easy. > >> Also I forgot that the real part of casinh(0+I*x) isn't always 0. If >> |x|>1, it is something non-zero. And so you need to check that >> creal(casinh(0+I*x)) and creal(casinh(-0+I*x)) have opposite signs in >> this case. > > This is related to my "domains can be different" comment. Adding code > to restrict the domain of the argument to be compatible with the real > function isn't too hard (off the top of my head, I think the domains > are all one of [-1,1], [0,Inf] or (0,Inf]). Handling behaviour > outside that domain requires more special-casing because the behaviour > is less consistent. I have a rather good handle on what that behavior will be for the complex arc-functions (since I have been working hard on them recently). casinh(x+I*y) and casin(x+I*y) have positive real and imaginary parts if x and y are positive. (Positive in this context includes 0, but not -0.) In particular: the sign of creal(casinh(z)) is the same as the sign of x; the sign of cimag(casinh(z)) is the same as the sign of y; the sign of creal(casin(z)) is the same as the sign of x; the sign of cimag(casin(z)) is the same as the sign of y; cacosh(z) and cacos(z) always have positive real part. The imaginary part of cacos(x+I*y) has the opposite sign to x (since it is PI/2 - casin(x+I*y) The imaginary part of cacosh(x+I*y) has the same sign as x. The signs for catanh and catan are exactly the same as for casinh and casin. For clog: The imaginary part of clog(x+I*y) has the same sign as y. The real part of clog will never be -0, and this doesn't have to be checked. I am fairly sure I got these correct. If your program starts spitting out huge numbers of errors, then I am wrong. But it won't take me long to figure out which ones I got wrong. > >>> I'm less sure of the next logical >>> step, which is to check things like >>> casinh(x + I*0) = asinh(x) + I*0 >> >> Does C99 mandate this? > > Nope. They are just mathematical equivalences (at least within the > domains supported by the real function). POLA implies that they > should be true but unless they are special-cased, the complex variant > probably has less accuracy as a result of the additional calculations > to support the imaginary component. > >> My programs probably won't satisfy this, because >> I realized that the computation works in these cases anyway. Of course, >> it would be easy to make it happen. > > It's probably up to the implementation - special casing pure real or > imaginary arguments should give those cases a shorter and simpler (and > therefore faster and more accurate) calculation but it's a matter of > whether this case in common enough to justify the additional test(s) > in all cases. > > It's also just occurred to me that doing so may result in unexpected > output discontinuities between cfoo(x-I*tiny), cfoo(x-I*0), cfoo(x+I*0) > and cfoo(x+I*tiny). Since the functions are correct within a few ULP, you shouldn't see any discontinuities like this. The algorithms I coded are meant to be rather good at handling these situations.