From owner-freebsd-standards@FreeBSD.ORG Fri Sep 30 15:27:40 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 82F6116A41F for ; Fri, 30 Sep 2005 15:27:40 +0000 (GMT) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.FreeBSD.org (Postfix) with ESMTP id 3196F43D48 for ; Fri, 30 Sep 2005 15:27:40 +0000 (GMT) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost [127.0.0.1]) by troutmask.apl.washington.edu (8.13.4/8.13.4) with ESMTP id j8UFRdnQ020840; Fri, 30 Sep 2005 08:27:39 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.13.4/8.13.1/Submit) id j8UFRYPG020839; Fri, 30 Sep 2005 08:27:34 -0700 (PDT) (envelope-from sgk) Date: Fri, 30 Sep 2005 08:27:34 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20050930152734.GA20655@troutmask.apl.washington.edu> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=unknown-8bit Content-Disposition: inline Content-Transfer-Encoding: 8bit In-Reply-To: <20050930221846.T34595@delplex.bde.org> User-Agent: Mutt/1.4.2.1i 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 15:27:40 -0000 On Sat, Oct 01, 2005 at 12:03:01AM +1000, Bruce Evans wrote: > 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. Good. I found a copy of n869 via google. Hopefully, the actual standard and the draft do not greatly disagree. As you may of guessed, I've put logl() on the back burner for the moment. I was hoping the complex.h functions would be easier to do. > >example, if z = x + I * y, then > > > >cos(z) = 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 = creal(z); > > y = 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) == +-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). tanh(y) can't be +-Inf. |tanh(y)| <= 1 for all y (except NaN), so cos(), sin(), and tanh() are bounded and only cosh() diverges to +Inf for large argument. A slightly better implementation than the above would be if (y == 0.) return (cos(x)); if (x == 0.) return (cosh(y)); return (cosh(y) * (cos(x) - I * sin(x) * tanh(y))); I can restore the sinh(y). > (*) 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) = ccosh(iz) > % ... > > Some functions, including ccos() are specified in terms of others like this. 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)) = 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. 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. > 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. > % -- ccosh(+0+i) returns NaNħi0 (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) = +inf+i0 because the pdf has no ++i0 case. > > % > % -- ccosh(++i) returns ++iNaN and raises the invalid > % exception. > % > % -- ccosh(x+i) returns NaN+iNaN and raises the invalid > % exception, for finite nonzero x. > > More bugs in n869.txt. > > % > % -- 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 (= 0) must not affect > the result, and it says something about signs. Again, the pdf file has ccosh(+Inf+iy) = +Inf cis(y) (finite nonzero y). which I infer means the sign on Inf is determined from cis(y). 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. -- Steve