[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: complex.h math functions
On Wed, 26 Oct 2005, David Schultz wrote:
On Thu, Sep 29, 2005, Steve Kargl wrote:
Is it permissible to implement the complex.h math
functions in terms of functions in math.h. For
example, if z = x + I * y, then
cos(z) = cos(x) * cosh(y) - I sin(x) * sinh(y)
This can be (naively?) implemented as
...
return (cosh(y) * (cos(x) - I * sin(x) * tanh(y)));
So I unfortunately don't have time to look into this right now,
but I have a few brief comments. Yes, all of the complex.h
functions can be implemented in terms of ordinary math.h functions
that operate on the real and imaginary parts. However, with the
"all" == only the C99 ones. The bessel, erf and gamma families
don't have any useful decomposition into real and imaginary parts
AFAIK.
exception of carg() and the functions already implemented in
FreeBSD, implementing them this way could result in significant
errors in corner cases. (For instance, underflow or overflow
might occur at points near the real or imaginary axes.) On the
other hand, doing better than this is no easy task.
We decided not to handle these cases better than ordinary complex
multiplication. Even z*z is hard to implement perfectly, and gcc
doesn't try. However, complex functions in the exp family are easier
to handle than complex multiplication, since they only require
real multiplication and no subtraction. In the above, the
real and imaginary parts can be made accurate independently just
by calculating the real functions accurately, and overflow is
fairly easy to avoid using an frexp() format for exp():
cosh(y) * cos(x) ~ (1/2 * exp(y)) * cos(x) for large y
= (1/2 * 2**N) * ((exp(y) / 2**N) * cos(x))
exp() already has the frexp()-like format internally and its last
step is to add N to the exponent. It just needs a higher overflow
threshold and its internals exported to work in the above (and extra
precision for it and sin() and the last "*" for the final relative
error to be < 1 ulp). The frexp() format is already needed to get
an error of < 1 ulp for real cosh():
cosh(x) ~ 1/2 * exp(y) for large y
= (1/2 * 2**N) * ((exp(y) / 2**N)
but the actual implementation is:
cosh(x) ~ 1/2 * exp(y) for large y
= 1/2 * (exp(y/2)*exp(y/2))
which may have an error of up to 2.5 ulps, just like naive ccosh() for
the non-overflow case. (Its actual max error for coshf() is more like
1.5 ulps).
[1] describes for complex arcsine and arccosine what the
troublesome cases are and how to do better. (You don't need to
Thanks; I'll look at it.
The GNU
C library just uses trivial implementations of the complex math
functions in terms of real math functions. In at least one case,
the implementation is completely bogus, which gives you an idea of
how much people are using this stuff at this point...
Cephes is better (much more complete and less bogus), but is still
very inaccurate. But obviously not many people use even cos(), else
the would notice that the error for i387 cos(x) near pi/2 is about
2**40 ulps :-).
In any case, I'm not up enough on this stuff to even remember what
a complex arctangent is supposed to look like, nor do I have time
to look into it right now. It's fine with me if you want to
Complex inverse trig functions are simpler in some ways than real
ones. They are all given by simple formulas involving other complex
functions. E.g.:
cacos(z) = -I*clog(z + csqrt(z*z - 1))
where the branches of log() and sqrt() must be chosen carefully. This
might be easier to use than a power series, except near z == 1.
Bruce
