From owner-freebsd-standards@FreeBSD.ORG Thu Oct 27 05:04:31 2005 Return-Path: X-Original-To: freebsd-standards@freebsd.org Delivered-To: freebsd-standards@freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id 948E016A41F; Thu, 27 Oct 2005 05:04:31 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailout1.pacific.net.au (mailout1.pacific.net.au [61.8.0.84]) by mx1.FreeBSD.org (Postfix) with ESMTP id E2EA943D45; Thu, 27 Oct 2005 05:04:30 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy1.pacific.net.au (mailproxy1.pacific.net.au [61.8.0.86]) by mailout1.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j9R54T2f005424; Thu, 27 Oct 2005 15:04:29 +1000 Received: from katana.zip.com.au (katana.zip.com.au [61.8.7.246]) by mailproxy1.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j9R54QVK022594; Thu, 27 Oct 2005 15:04:27 +1000 Date: Thu, 27 Oct 2005 15:04:26 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: David Schultz In-Reply-To: <20051026190015.GA52635@VARK.MIT.EDU> Message-ID: <20051027145316.Y23976@delplex.bde.org> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20051026190015.GA52635@VARK.MIT.EDU> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@freebsd.org, Steve Kargl Subject: Re: complex.h math functions X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Thu, 27 Oct 2005 05:04:31 -0000 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