From owner-freebsd-standards@FreeBSD.ORG Sat Oct 1 11:38:53 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 1F87A16A41F for ; Sat, 1 Oct 2005 11:38:53 +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 732A043D46 for ; Sat, 1 Oct 2005 11:38:52 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy2.pacific.net.au (mailproxy2.pacific.net.au [61.8.0.87]) by mailout1.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j91BclJ6025302; Sat, 1 Oct 2005 21:38:47 +1000 Received: from katana.zip.com.au (katana.zip.com.au [61.8.7.246]) by mailproxy2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j91BciJ2020242; Sat, 1 Oct 2005 21:38:45 +1000 Date: Sat, 1 Oct 2005 21:38:44 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20050930152734.GA20655@troutmask.apl.washington.edu> Message-ID: <20051001204005.N37704@delplex.bde.org> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <20050930152734.GA20655@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: MULTIPART/MIXED; BOUNDARY="0-1263416119-1128166724=:37704" 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: Sat, 01 Oct 2005 11:38:53 -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-1263416119-1128166724=:37704 Content-Type: TEXT/PLAIN; charset=X-UNKNOWN; format=flowed Content-Transfer-Encoding: QUOTED-PRINTABLE On Fri, 30 Sep 2005, Steve Kargl wrote: > On Sat, Oct 01, 2005 at 12:03:01AM +1000, Bruce Evans wrote: >> On Thu, 29 Sep 2005, Steve Kargl wrote: >>> 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 bug= s >> 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 multultiplicat= ion >> that doesn't take much longer (relatively) but loses precision (relative= ly >> more). > > tanh(y) can't be +-Inf. |tanh(y)| <=3D 1 for all y (except NaN), > so cos(), sin(), and tanh() are bounded and only cosh() diverges > to +Inf for large argument. Oops. > A slightly better implementation than > the above would be > > if (y =3D=3D 0.) return (cos(x)); > if (x =3D=3D 0.) return (cosh(y)); > return (cosh(y) * (cos(x) - I * sin(x) * tanh(y))); I'm not sure if this is better. The special cases may be so rare that checking for them slows down the average case. I think special handling is unnecessary since cos(0) is exactly 1.0, etc, so the checks are just an optimization of the special args. OTOH, the type-generic version should have special checks using something like "if (__builtin_const_p(y) && y =3D=3D 0.)" to classify the pure real case, etc. Inline function args are never constant, so the type-generic version must be implemented as a macro. > ... > 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)) =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 preserve= d >> _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. My main point is that even with round-to-nearest, it is unclear that the rounding goes the right way for standard formulae to work right. >> 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. Except signed zeros (and infinities) don't really work for the complex case. Signed zeros give 8 directions (multiples of Pi/4) for approaching the origin, but there are an infinite number of possible directions (all angles). >> % -- 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. >> >> % >> % -- 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) =3D +inf+i0 > > because the pdf has no ++i0 case. Yes, the infinity sign apparently got lost in translation to ASCII. The previous one, ccosh(+0+i), is actually ccosh(+0+i*Inf), but doesn't look as weird when mistranslated since the infinity is lost at the end. You could also look at the POSIX spec. It is more up to date and has an index. However, it seems to be missing the detail spec for the complex functions in at least the text version of the 2001 draft (this specifies handling of NaNs in detail only for the real functions, and also uses the \xb1 character which I don't understand yet). > 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. Float complex can use double precision with tiny or negative loss of performance given our current low performance float functions. Also, long double precision is almost free for the double exp() family (exp/log/trig) on some machines (i386 and amd64) that can do these functions in hardware given our current not very high performance (hardware or not) double functions. I would use long double (real) versions for the initial implementation of the more exotic complex functions if possible. Unfortunately, it doesn't seem to be possible to implement the more exotic complex functions in terms of real functions. Bruce --0-1263416119-1128166724=:37704--