Skip site navigation (1)Skip section navigation (2)
Date:      Sun, 12 Aug 2012 23:13:46 -0000
From:      Bruce Evans <brde@optusnet.com.au>
To:        Stephen Montgomery-Smith <stephen@missouri.edu>
Cc:        Diane Bruce <db@db.net>, Steve Kargl <sgk@troutmask.apl.washington.edu>, John Baldwin <jhb@freebsd.org>, David Chisnall <theraven@freebsd.org>, Bruce Evans <bde@freebsd.org>, Bruce Evans <brde@optusnet.com.au>, 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:  <20120719144432.N1596@besplex.bde.org>
Resent-Message-ID: <20120812231338.GD20453@server.rulingia.com>
In-Reply-To: <5006D13D.2080702@missouri.edu>
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> <5006D13D.2080702@missouri.edu>

next in thread | previous in thread | raw e-mail | index | archive | help
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



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120719144432.N1596>