complex.h math functions

[Steve Kargl]

On Sat, Oct 01, 2005 at 12:03:01AM +1000, Bruce Evans wrote:
> On Thu, 29 Sep 2005, Steve Kargl wrote:
> 
> >Is it permissible to implement the complex.h math
> >functions in terms of functions in math.h.  For
> 
> There are no namespace problems AFAIK, and the C standard doesn't
> require any particular quality of implementation.

Good.  I found a copy of n869 via google.  Hopefully, the
actual standard and the draft do not greatly disagree.

As you may of guessed, I've put logl() on the back burner
for the moment.  I was hoping the complex.h functions 
would be easier to do.

> >example, if z = x + I * y, then
> >
> >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.  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 can restore the sinh(y).

> (*) The best treatment of exceptional args is unclear, but C99 specifies
> a treatment in detail in at least n869.txt.  This spec probably rules out
> simply using formulae as above -- exceptional args need to be handled
> specially unless you can be sure that applying the formulae results in
> the specified NaNs, etc.  From n869.txt:
> 
> %                ccos(z)     = ccosh(iz)
> % ...
> 
> Some functions, including ccos() are specified in terms of others like this.

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.

> 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.

> %          -- 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.

> 
> % 
> %          -- ccosh(++i)  returns  ++iNaN  and  raises  the   invalid
> %             exception.
> % 
> %          -- ccosh(x+i)  returns  NaN+iNaN  and  raises  the invalid
> %             exception, for finite nonzero x.
> 
> More bugs in n869.txt.
> 
> % 
> %          -- ccosh(++iy) returns +cis(y), for finite nonzero y.
> 
> cis(y) is (cos(y) + I * sin(y)).  This says that if the arg is pure
> imaginary then the implicit real part of the arg (= 0) must not affect
> the result, and it says something about signs.

Again, the pdf file has

ccosh(+Inf+iy) = +Inf cis(y) (finite nonzero y).

which I infer means the sign on Inf is determined from cis(y).

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.

-- 
Steve