complex.h math functions

[Bruce Evans]

On Fri, 30 Sep 2005, Steve Kargl wrote:

> On Sat, Oct 01, 2005 at 12:03:01AM +1000, Bruce Evans wrote:
>> On Thu, 29 Sep 2005, Steve Kargl wrote:
>>> cos(z) = cos(x) * cosh(y) - I sin(x) * sinh(y)
>>>
>>> This can be (naively?) implemented as
>>>
>>> double complex
>>> ccos(double complex z)
>>> {
>>>  double x, y;
>>>  x = creal(z);
>>>  y = cimag(y);
>>>  return (cosh(y) * (cos(x) - I * sin(x) * tanh(y)));
>>> }
>>
>> The original formula gives a much better implementation.  It is missing
>> serious bugs at points where tanh(y) == +-Inf; it may be missing bugs
>> at points where y is a NaN (*); tanh(y) takes more work to compute than
>> sinh(y), and then multiplying it by cosh(y) gives 1 extra multultiplication
>> that doesn't take much longer (relatively) but loses precision (relatively
>> more).
>
> tanh(y) can't be +-Inf.  |tanh(y)| <= 1 for all y (except NaN),
> so cos(), sin(), and tanh() are bounded and only cosh() diverges
> to +Inf for large argument.

Oops.

> A slightly better implementation than
> the above would be
>
>        if (y == 0.) return (cos(x));
>        if (x == 0.) return (cosh(y));
>        return (cosh(y) * (cos(x) - I * sin(x) * tanh(y)));

I'm not sure if this is better.  The special cases may be so rare that
checking for them slows down the average case.  I think special handling
is unnecessary since cos(0) is exactly 1.0, etc, so the checks are just
an optimization of the special args.  OTOH, the type-generic version
should have special checks using something like
"if (__builtin_const_p(y) && y == 0.)" to classify the pure real case,
etc.  Inline function args are never constant, so the type-generic
version must be implemented as a macro.

> ...
> OK, I'll use the implementation details of ccosh to infer the
> handling of exceptional values.
>
>> %        G.5.2.4  The ccosh functions
>> %
>> %          -- ccosh(conj(z)) = conj(ccosh(z)) and ccosh is even.
>>
>> A strict reading says that we must ensure that if a formula is used, the
>> rounding must always go in such a way that these symmetries are preserved
>> _exactly_.  Rounding to nearest may give this, but this is not clear.
>> Rounding towards Inf obviously wouldn't give this.
>
> I haven't looked closely enough at libm, but I suspect that it
> it at best handles round-to-nearest.  The default rounding in
> FreeBSD is round-to-nearest.  To place greater restrictions on
> the complex.h with respect to proper handling of all rounding
> modes than we have on libm would seem to be an undue requirement.
> If a programmer sets another rounding mode, then its her responsibility
> o test that libm has the expected behavior.

My main point is that even with round-to-nearest, it is unclear that
the rounding goes the right way for standard formulae to work right.

>> I just realized that many of the complications here are for type-generic
>> functions.  0.0 is quite different from (0.0 + I * 0.0) since it has
>> type double but the latter has type Complex.  When combined with +0
>> being slightly different from -0, this gives zillions of cases.
>
> Signed zero are going to be important especially for those complex.h
> functions with branch cuts (e.g., csqrt).  We do not want to return
> a value on the wrong Riemann sheet.  Of course, this is an implementation
> detail that I'll worry about when I consider a fucntion with a branch
> cut.

Except signed zeros (and infinities) don't really work for the complex
case.  Signed zeros give 8 directions (multiples of Pi/4) for approaching
the origin, but there are an infinite number of possible directions (all
angles).

>> %          -- ccosh(+0+i) returns  NaN±i0  (where  the  sign  of  the
>> %             imaginary part of the result is unspecified) and raises
>> %             the invalid exception.
>>
>> This is apparently a bug in n869.txt.  ccosh() on finite values is never
>> a NaN.
>>
>> %
>> %          -- ccosh(++i0) returns ++i0.
>>
>> Here ++ is not the C operator, but something like +-.  It may be another
>> bug that this says ++ and not +-.  Shouldn't the result reverse the sign?
>
> I found n869.pdf.  I wonder if the above is suppose to be
>
> ccosh(+Inf+i0) = +inf+i0
>
> because the pdf has no ++i0 case.

Yes, the infinity sign apparently got lost in translation to ASCII.
The previous one, ccosh(+0+i), is actually ccosh(+0+i*Inf), but doesn't
look as weird when mistranslated since the infinity is lost at the end.

You could also look at the POSIX spec.  It is more up to date and has
an index.  However, it seems to be missing the detail spec for the
complex functions in at least the text version of the 2001 draft (this
specifies handling of NaNs in detail only for the real functions, and
also uses the \xb1 character which I don't understand yet).

> Thanks for visiting the accuracy issue.  I came to essentially
> the same conclusion that you have.  If we had reasonable
> implementations of the long double complex function, we could
> use these to compute float complex and double complex to achieve
> better accuracy.  Of course, performance would suffer.

Float complex can use double precision with tiny or negative loss of
performance given our current low performance float functions.  Also,
long double precision is almost free for the double exp() family
(exp/log/trig) on some machines (i386 and amd64) that can do these
functions in hardware given our current not very high performance
(hardware or not) double functions.  I would use long double (real)
versions for the initial implementation of the more exotic complex
functions if possible.  Unfortunately, it doesn't seem to be possible
to implement the more exotic complex functions in terms of real
functions.

Bruce