Date: Sun, 12 Aug 2012 23:06:37 -0000 From: Stephen Montgomery-Smith <stephen@missouri.edu> To: Bruce Evans <brde@optusnet.com.au> Cc: Diane Bruce <db@db.net>, John Baldwin <jhb@freebsd.org>, David Chisnall <theraven@freebsd.org>, Bruce Evans <bde@freebsd.org>, Steve Kargl <sgk@troutmask.apl.washington.edu>, David Schultz <das@freebsd.org>, Peter Jeremy <peter@rulingia.com>, Warner Losh <imp@bsdimp.com> Subject: Re: Use of C99 extra long double math functions after r236148 Message-ID: <5006D13D.2080702@missouri.edu> Resent-Message-ID: <20120812230630.GZ20453@server.rulingia.com> In-Reply-To: <20120717200931.U6624@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> <C527B388-3537-406F-BA6D-2FA45B9EAA3B@FreeBSD.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>
next in thread | previous in thread | raw e-mail | index | archive | help
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.
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?5006D13D.2080702>