From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:06:37 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 AF41D106566B for ; Sun, 12 Aug 2012 23:06:37 +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 1ABD58FC08 for ; Sun, 12 Aug 2012 23:06:36 +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 q7CN6aat075716 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:06:37 +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 q7CN6U4H021369 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:06:30 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CN6Ugc021368 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:06:30 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:06:30 +1000 Resent-Message-ID: <20120812230630.GZ20453@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 q6IF8Ebf098382 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Thu, 19 Jul 2012 01:08:14 +1000 (EST) (envelope-from stephen@missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q6IF8BGY073152 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Thu, 19 Jul 2012 01:08:13 +1000 (EST) (envelope-from stephen@missouri.edu) 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 q6IF7fou016669; Wed, 18 Jul 2012 10:07:42 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <5006D13D.2080702@missouri.edu> From: Stephen Montgomery-Smith Mail-Followup-To: freebsd-numerics@freebsd.org User-Agent: Mozilla/5.0 (X11; Linux i686; rv:13.0) Gecko/20120615 Thunderbird/13.0.1 MIME-Version: 1.0 To: Bruce Evans References: <20120529045612.GB4445@server.rulingia.com> <20120711223247.GA9964@troutmask.apl.washington.edu> <20120713114100.GB83006@server.rulingia.com> <201207130818.38535.jhb@freebsd.org> <9EB2DA4F-19D7-4BA5-8811-D9451CB1D907@theravensnest.org> <20120713155805.GC81965@zim.MIT.EDU> <20120714120432.GA70706@server.rulingia.com> <20120717084457.U3890@besplex.bde.org> <5004A5C7.1040405@missouri.edu> <5004DEA9.1050001@missouri.edu> <20120717200931.U6624@besplex.bde.org> In-Reply-To: <20120717200931.U6624@besplex.bde.org> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit Cc: Diane Bruce , 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:06:37 -0000 X-Original-Date: Wed, 18 Jul 2012 10:07:41 -0500 X-List-Received-Date: Sun, 12 Aug 2012 23:06:37 -0000 I went on a long road trip yesterday, so I didn't get any code written, but I did have a lot of thoughts about clog and casinh. First, the naive formula (here z=x+I*y) clog(z) = cpack(log(hypot(x,y)),atan2(x,y)) is going to work in a lot more edge cases then one might expect. This is because hypot and atan2, especially atan2, already do a rather good job getting the edge cases right. I am thinking in particular of when x or y are 0 or -0, or one of them is infinity or -infinity. So writing this code will be quite a bit easier than I was expecting. It looks like I will just have to worry about when both x and y are infinite, or when one or both of them are nan. Next, concerning casinh: On 07/17/2012 06:13 AM, Bruce Evans wrote: > I translated this to pari. There was a sign error for the log() term > (assuming that pari asinh() has the same branch cuts as C99): I couldn't spot which sign error Bruce had changed. However I expect it has something to do with what happens when x=0 and fabs(y)>1. This is the reasonable choice of the branch cut. What I think the value should be is casinh(z) = cpack( signum(x)*sqrt(fabs(y)+sqrt(y^2-1)), signum(y)*PI) where the value of signum(x) depends on whether x is 0 or -0. (I might add that I checked against the Mathematica ArcSinh function, and this does NOT follow the above rule. But the document Steve pointed me to says that casinh(conj(z)) = conj(casinh(z)) which means that we cannot follow the Mathematica conventions.) > The most obvious immediate difficulty in translating the above into C is > that y*y and z*z may overflow when the result shouldn't. This will be a lot easier than I originally expected. When we are in conditions when overflow might occur, we can simply make the approximations sqrt(y*y-1) = y csqrt(z*z+1) = signum(x)*z because in floating point arithmetic, these will not be approximations, but true exactly. And I am thinking that the test I will use for when to use these approximations will be (y==y+1) and (z==z+1) respectively. (I would use (z*z==z*z+1) but that test has the overflow problem.) Finally, I want to tell you guys that the reason I used the code: if (x>0) return clog(z+csqrt(z*z+1)); else return -clog(-z+csqrt(z*z+1)); is this. Both formulas are mathematically exactly the same. This is true even if one takes into account the branch cuts for csqrt and clog. The difference between the two formulas is numerical errors. For example, if x<0 and z has very large magnitude, then csqrt(z*z+1) will be very close to -z. In fact in floating point arithmetic, if the magnitude of z is sufficiently large, they will be the same. However, as I am typing this, I realize that the code should really be if (w!=z+1) { w = z*z+1; if (signum(creal(w))==1) return clog(z+csqrt(w)); else return -clog(-z+csqrt(w)); } else /* if (z==z+1) */ { if (x>0) return clog(2*z); else return -clog(-2*z); } where the signum function is defined so that signum(0)==1 and signum(-0)=-1. Next: cacosh and cacos. I had presumed the formula cacosh(z) = I*cacos(z) which can be true depending on how one defines the branch cuts. But this formula won't satisfy the C99 standard which mandates cacosh(conj(z)) = conj(cacosh(z)) That is why in earlier posts I thought there was a mistake in the online documentation, and the range of outputs of cacosh should satisfy imaginary part in [0,pi] rather than [-pi,pi]. Anyway, I am just posting this update so that you see I am thinking about it. I might get the project done in days, but it might also be months. It depends upon how much other stuff I have going on.