Daemon News Ezine | BSD News | BSD Mall | BSD Support Forum | BSD Advocacy | BSD Updates |

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

*To*: David Schultz <das@xxxxxxxxxxx>*Subject*: Re: complex.h math functions*From*: Bruce Evans <bde@xxxxxxxxxxx>*Date*: Thu, 27 Oct 2005 15:04:26 +1000 (EST)*Cc*: freebsd-standards@xxxxxxxxxxx, Steve Kargl <sgk@xxxxxxxxxxxxxxxxxxxxxxxxxxxx>*Delivered-to*: freebsd-standards@xxxxxxxxxxx*In-reply-to*: <20051026190015.GA52635@xxxxxxxxxxxx>*List-archive*: <http://lists.freebsd.org/pipermail/freebsd-standards>*List-help*: <mailto:freebsd-standards-request@freebsd.org?subject=help>*List-id*: Standards compliance <freebsd-standards.freebsd.org>*List-post*: <mailto:freebsd-standards@freebsd.org>*List-subscribe*: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, <mailto:freebsd-standards-request@freebsd.org?subject=subscribe>*List-unsubscribe*: <http://lists.freebsd.org/mailman/listinfo/freebsd-standards>, <mailto:freebsd-standards-request@freebsd.org?subject=unsubscribe>*References*: <20050929195552.GA14982@xxxxxxxxxxxxxxxxxxxxxxxxxxxx> <20051026190015.GA52635@xxxxxxxxxxxx>

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

**References**:**complex.h math functions***From:*Steve Kargl

**Re: complex.h math functions***From:*David Schultz

- Prev by Date:
**Re: complex.h math functions** - Next by Date:
**Re: complex.h math functions** - Previous by thread:
**Re: complex.h math functions** - Next by thread:
**Re: kern/64875: [libc] [patch] [feature request] add a system call: fdatasync()** - Index(es):