From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:13:46 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 ECC5A106566B for ; Sun, 12 Aug 2012 23:13:45 +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 720888FC17 for ; Sun, 12 Aug 2012 23:13:45 +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 q7CNDjlQ075977 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:13:45 +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 q7CNDcnQ021937 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:13:38 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CNDcJV021936 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:13:38 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:13:38 +1000 Resent-Message-ID: <20120812231338.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 q6J6cJUx011155 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Thu, 19 Jul 2012 16:38:19 +1000 (EST) (envelope-from brde@optusnet.com.au) Received: from mail06.syd.optusnet.com.au (mail06.syd.optusnet.com.au [211.29.132.187]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q6J6cJmp079623 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Thu, 19 Jul 2012 16:38:19 +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 mail06.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q6J6bqRG014261 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Thu, 19 Jul 2012 16:37: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: <5006D13D.2080702@missouri.edu> Message-ID: <20120719144432.N1596@besplex.bde.org> 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> <5006D13D.2080702@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Mailman-Approved-At: Sun, 12 Aug 2012 23:56:00 +0000 Cc: Diane Bruce , Steve Kargl , John Baldwin , David Chisnall , Bruce Evans , Bruce Evans , 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:13:46 -0000 X-Original-Date: Thu, 19 Jul 2012 16:37:52 +1000 (EST) X-List-Received-Date: Sun, 12 Aug 2012 23:13:46 -0000 On Wed, 18 Jul 2012, Stephen Montgomery-Smith wrote: > 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. Right, clog is deceptively simple. This is because it decomposes perfectly into 2 real functions of 2 real variables and both of these functions are standard and already implemented almost as well as possible. ISTR das saying that it had a complicated case, but I don't see even one. atan2() is supposed to handle all combinations of +-0. Now I remember a potential problem. Complex functions should have only poles and zeros, with projective infinity and "projective zero" (= inverse of projective infinity). Real functions can and do have affine infinities and zeros (+-Inf and +-0), with more detailed special cases. It's just impossible to have useful, detailed special cases for all the ways of approaching complex (projective) infinity and 0. I think Kahan wanted projective infinity in IEEE7xx in ~1980. Intel 8087 had both projective infinity and affine infinities, but projective infinity didn't make it into the first IEEEE7xx, and hardly anyone understood it and it was eventually dropped from Intel FPUs (I think it was in 80287; then in i486 it was reduced to a bit in the control word that can never be cleared (the bit is to set affine infinities); then in SSE the bit went away too). However, C99 tries too hard to make complex functions reduce to real functions when everything is purely real or purely complex. So most of the special cases for +-0 and +-Inf affect complex functions (for other directions of approaching 0 and infinity, not much is specified but you should try to be as continuous as possible, where continuity has delicate unclear meanings since it is related to discontinuous sign functions). Hopefully, the specification of imag(clog()) is that it has the same sign behaviour as atan2(), so you can just use atan2(). The sign conventions for both are arbitrary, but they shouldn't be gratuitously different. You still have to check that they aren't non-gratuitously different, because different conventions became established. cexp() was not quite as simple as clog(). In polar coordinates for the result it is even simpler: cexp(x, y) = (exp(x), y) (the real functions are independent), but in affine coordinates you have to multiply real functions and avoid underflow and overflow in the multiplications. clog() is simpler than that since it an addition instead of multiplications, and the addition doesn't even mix the real functions, but keeps them separate as real and immaginary parts. The real functions are more complicated since they take 2 args, but they are already implemented. > 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.) I had forgotten that pari doesn't support -0 at all (AFAIK). I certainly had to change a sign to get match the pari result for 0+y*I, but it was the sign of y. Your original code seems to have y where it should have x: > if (x==0) { > if (fabs(y)<=1) return I*asin(y); > else return signum(y)* ( > log(fabs(y)+sqrt(y*y-1)) > + I*PI/2); > } Here x is +-0 and there is no sign test for it. Pari probably doesn't distinguish +-0, and produces a result for +0. You choose a sign that only depends on and on conventions for asin() when fabs(y) <= 1. Then the choice of signs given by signum(y) is necessary for making the imaginary part agree with the choice when fabs(y) <= 1. The sign error that I got was for the real part: pari asinh(2*I) = -1.316... + Pi/2*I above asinh(2*I) = 1.316... + Pi/2*I C99 doesn't specify casinh(0+y*I) directly, However, it specifies that casinh(x+Inf*I) = +Inf + Pi/2*I for finite positive x, so it needs your choice of sign to be continuous for y in [1, Inf]. C99 specifies lots of weird behaviours for various affine infinities, but all of them have nonnegative real parts, so the pari choice of sign is not wanted. >> 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.) The problem is that y*y overflows when casinh() doesn't. A simpler case is that |y| never overflows for finite y, but sqrt(y*y) overflows at about sqrt(DBL_MAX). So sqrt(y*y) is quite different from |y| in floating point, though it is the same in real arithmeric. Adding +-1 to y*y moves the points of spurious and real overflow slightly. > 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. Probably even more sophistication is needed. I don't have much idea of the actual errors involved. At some point, someone should do an error analysis, and or test the worst cases to verify that the error is not unreasonably large (say < 100 ulps). BTW, Intel and AMD docs have always claimed errors of < 1 ulp for trig functions, but the actual errors are multi-giga-ulps near multiples of Pi/2 (my tests routinely find ones of 17 Gulps in float precision and 4503 Gulps in double precision). Such errors are normal near zeros of analytic functions unless the zeros are known to a precision of hyndreds or thousands of bits and extra-precision code with this many bits is used near them. Intel and AMD only use ~68 bits for Pi/2 in i387 functions. FreeBSD (fdlibm) uses the necessary thousands of bits for Pi/2. FreeBSD (fdlibm) doesn't do this for zeros of Bessel or gamma functions, so it has large errors for them too (I see the same errors for Bessel functions as for trig functions, but for lgammaf I only see errors 10 mega-ulps routinely; after switching from the i387 trig functions, I only see errors of 2 mega-ulps for Bessel functions). exp() has no zeros, so it doesn't have the probles of trig functions (or rather, it has essentially the same problems, but you agree not to notice them by looking at it in polar coordinates: exp(y*I) = cos(y) + sin(y)*I Near a zero of cos(y), sin(y) is +-1 and the relative error of exp(y*I) is small relative to |exp(y*I)| = 1, though it may be 4503 Gulps for the relative error of the real part relative to the real part. For cexp(), we easily implement accurate real and immaginary parts relative to themselves, limited mainly by the accuracy of the real functions, since the real functions are mixed nicely. Similarly for clog(). This is too hard for a general analytic function. I don't know what happens with zeros of complex inverse trig functions. I think they don't have many (like log()), but their real and imaginary parts do, and they are too general for accurate behaviour of the real and imaginary parts relative to themselves to fall out. The arbitrary mediocre error bound of 100 ulps is modestly less than 4503 Gulps. Bruce