From owner-freebsd-standards@FreeBSD.ORG Fri Sep 30 14:03:07 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 B63D916A41F for ; Fri, 30 Sep 2005 14:03:07 +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 1952943D49 for ; Fri, 30 Sep 2005 14:03:06 +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 j8UE34Jn001971; Sat, 1 Oct 2005 00:03:04 +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 j8UE31MA017499; Sat, 1 Oct 2005 00:03:02 +1000 Date: Sat, 1 Oct 2005 00:03:01 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20050929195552.GA14982@troutmask.apl.washington.edu> Message-ID: <20050930221846.T34595@delplex.bde.org> References: <20050929195552.GA14982@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: MULTIPART/MIXED; BOUNDARY="0-1213393512-1128088981=:34595" Cc: freebsd-standards@freebsd.org 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: Fri, 30 Sep 2005 14:03:07 -0000 This message is in MIME format. The first part should be readable text, while the remaining parts are likely unreadable without MIME-aware tools. --0-1213393512-1128088981=:34595 Content-Type: TEXT/PLAIN; charset=X-UNKNOWN; format=flowed Content-Transfer-Encoding: QUOTED-PRINTABLE 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. > example, if z =3D x + I * y, then > > cos(z) =3D 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 =3D creal(z); > y =3D 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) =3D=3D +-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). (*) 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) =3D ccosh(iz) % ... Some functions, including ccos() are specified in terms of others like this= =2E % G.5.2.4 The ccosh functions %=20 % [#1] %=20 % -- ccosh(conj(z)) =3D 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. In the following, "I" is spelled "i" and "*" (multiplication) is omitted. %=20 % -- ccosh(+0+i0) returns 1+i0. The formula should give this automatically provided the real functions aren't broken. Handly of -0 and -i0 (in various combinations with +0's) doesn't seem to be specified for this function, although it is for others. The sign of the real 0 doesn't matter because the result is 1. I suppose the sign for ccosh(-i0) should track the sign of the arg in the same way as for similar real functions. Since ccosh(1) =3D=3D 1, the sign of the real 0 can't affect the sign of the imaginary part of the result, unlike for some other functions. 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. %=20 % -- ccosh(+0+i) returns NaN=B1i0 (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. %=20 % -- 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? %=20 % -- ccosh(++i) returns ++iNaN and raises the invalid % exception. %=20 % -- ccosh(x+i) returns NaN+iNaN and raises the invalid % exception, for finite nonzero x. More bugs in n869.txt. %=20 % -- 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 (=3D 0) must not affect the result, and it says something about signs. %=20 % -- ccosh(+0+iNaN) returns NaN=B1i0 (where the sign of the % imaginary part of the result is unspecified). This requires not propagating a NaN in the imaginary part to the real part if the real part is precisely 0. %=20 % -- ccosh(++iNaN) returns ++iNaN. The case of a pure imaginary NaN may be slightly different. I don't understand the non-ASCII character in the spec for ccosh(+0+iNaN). It is '\xb1'. %=20 % -- ccosh(x+iNaN) returns NaN+iNaN and optionally raises % the invalid exception, for finite nonzero x. Nonzero x causes the NaN to be propagated to the real part. Formulae tend to do this automatically; it is the previous 2 cases that using the formula for ccosh() breaks. %=20 % -- ccosh(NaN+i0) returns NaN=B1i0 (where the sign of the % imaginary part of the result is unspecified). %=20 % -- ccosh(NaN+iy) returns NaN+iNaN and optionally raises % the invalid exception, for all nonzero numbers y. Similarly if only the real part of the arg is a NaN. %=20 % -- ccosh(NaN+iNaN) returns NaN+iNaN. This requirement is obviously right and is hard to avoid satisfy using formulae. About accuratcy of complex functions: It is unreasonable to require or expect the result to be as accurate as possible for the real and imaginary parts independently, since even ordinary operations don't have this propery. E.g., consider the simple polynomial function f(z) =3D z*z. For z =3D x+Iy, f(z) =3D (x*x-y*y)+2Ix*y= =2E Here 2Ix*y is almost as accurate as x and y, but x*x-y*y only has 1 or 2 bits of accuracy when x is as close as possible to y (but different). It might be useful for cpow(z, 2) to be much more accurate than z*z, but this would be hard to implement. Standard polynomial approximation methods just won't give enough accuracy. OTOH, for cexp() and probably for most related functions, equal (relative) accuracy for the real and imaginary parts occurs automatically using standard formulae (since cexp(z) depends on Re(z) and Im(z) independently): cexp(z) =3D exp(x) * (cos(x) + I * sin(x)) The multiplication expands the roundoff errors for the real and imaginary parts independently, so the final relative error for each is < 2 ulps if the errors for exp(), cos() and sin() are always < 1 ulp. I'm not happy with an error of 2 ulps of course, but don't know any morwe accurate way implementing cexp() short of doing everything in more precision. Standard polynomial approximations with complex args might be more accurate=20 relative to |cexp(z)| =3D exp(x), but couldn't be very accurate for the real and imaginary parts independently. This corresponds to the arg reduction requirements for cos(x) and sin(x) being very different. Remez methods also seem to be completely unsuitable for for complex args since according to my limited understanding of them they depend on delicate cancellations which are only likely to occur when their complex arg is in a discrete set of directions. For ccos(z) implemented as: ccos(z) =3D cos(x) * cosh(y) - I sin(x) * sinh(y) the error is again < 2 ulps for each of the real and imaginary parts provided it is < 1 ulp for all the real functions. Bruce --0-1213393512-1128088981=:34595--