From owner-freebsd-standards@FreeBSD.ORG Wed Oct 5 03:24:10 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 C152116A41F for ; Wed, 5 Oct 2005 03:24:10 +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 6827343D45 for ; Wed, 5 Oct 2005 03:24:10 +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 j953O64t006902; Tue, 4 Oct 2005 20:24:06 -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 j953O0Lp006901; Tue, 4 Oct 2005 20:24:00 -0700 (PDT) (envelope-from sgk) Date: Tue, 4 Oct 2005 20:24:00 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20051005032400.GA6736@troutmask.apl.washington.edu> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <20050930152734.GA20655@troutmask.apl.washington.edu> <20051001204005.N37704@delplex.bde.org> <20051002191841.GA40568@troutmask.apl.washington.edu> <20051004164618.U47262@delplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20051004164618.U47262@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: Wed, 05 Oct 2005 03:24:10 -0000 On Tue, Oct 04, 2005 at 07:31:43PM +1000, Bruce Evans wrote: > On Sun, 2 Oct 2005, Steve Kargl wrote: > > >One question: Can we use a #define for complex float or do we > >need an trivial C program? That is, in math.h can we do > > > >#define ccoshf(z) ccosh((float complex)(z)) > > A program is needed, since (ccoshf)(z) and &(ccoshf) must work even if > ccoshf() is a macro. OK. > I converted to ccoshf() for testing it and comparing with a simpler > implementation that uses the formula, then converted back. There are > lots of special cases for -0, Infs and NaNs, and there were lots of > bugs, including compiler bugs (gcc is very far from supporting the IEC > 60559 extension). It turns out that the version using formula handles > most of the special cases automatically and is much easier to fix. I > fixed both versions and worked around the compiler bugs, and compared > differences to find the best version of the variant using the formula. So which implementation do you prefer: your simpler version or your fixed version of my first cut at ccosh? Apologies for overlooking -0. > % --- s_ccosh.c~ Mon Oct 3 07:10:59 2005 > % +++ s_ccosh.c Tue Oct 4 16:24:06 2005 > % @@ -40,11 +40,11 @@ > % #include > % > % -#include "math_private.h" > % +#include "../math_private.h" > > Temporary hack. I use a hack in my Makefile. CFLAGS+=-I/usr/src/lib/msun/src. > % double complex > % ccosh(double complex z) > % { > % - int con = 1; > % - double x, y; > % + double complex f; > % + double con, x, y; > % int32_t hx, hy, ix, iy, lx, ly; > % > > Optimize `con' a little, and fix some style bugs (initialization in > declaration, and disordered declarations; most other style bugs are not > fixed or noted). `f' is used to work around gcc bugs. You might say that the style violations are the MUTT style. KNF == kernel normal format (FreeBSD style) KNF2 == kargl normal format (My style) GNU == GCC's FSF style fdlibm == msun style I'll eventually make it conform to KNF. > > I didn't bother optimizing isnan(x), etc. using hx, etc. > For a committable version do you want/prefer the isnan or the bit twiddling? > Handle x == 0, y finite. There are 2 unobvious details: > - getting the sign right when y == 0... > > % + creal(f) = cos(y); > % + cimag(f) = x * con; > % + return f; Wow. I'm used to programming in Fortran and would never have written "creal(f) = ...". This looks like your assigning a value to a function. (For those in the peanut gallery that snicker at the mention of Fortran, please see the Fortran 2003 standard.) > - working around gcc bugs. :-) > Fortunately, we can construct u + I*v by assigning to the real and complex > parts. There should be a macro or inline function to do this. If we go the macro route, do you want it (them?) in math_private.h or complex.h? If creal(f) appears on the LHS, is it a generic reference so that type is forced to match the RHS? Is #define CMPLX((z),(x),(y)) {creal((z)) = (x), cimag((z)) = (y)} acceptable? > The (cosh(x) * tanh(x)) term should be written as sinh(x), especially now > that cosh(x) is evaluated twice. The old way of writing the return value: > > cosh(x) * (cos(y) + I * con * tanh(x) * sin(y)) > > has the interesting property of avoiding some of the compiler bugs. > The second term is always finite, so gcc doesn't mess up Infs and Nans > in it, and multiplying it by an infinite real coshx() does less damage > than multiplying I by an infinite real (sinh(x) * con * siny()).. While I did not know about this property per se. I realized that cos(), sin(), and tanh() are all bounded functions for finite arguments, and cosh() is the only function that could diverge (for finite floating point arguments). > y is Inf or Nan, so we can use the simpler formula y - y to generate a > NaN for it. ((y - y) / (y - y)) is only needed to generate a NaN from > a possibly-finite y. OK. > % } > % + abort(); > > This should be unreachble, so the previous if statement must always be > redundant. Yes, the previous if statement is redundant. I had filtered the exceptional cases first in my first cut implementation, and the finite, nonzero, x and y case were done last. This is how fdlibm appears to things for the float and double functions. With the large number of exceptional cases, I thought that we should handled those last (as an optimization). (2nd version deleted). > Third, a simple test program that compares implementations of ccosh() > on a reasonably large set of float complex args (including +0, +-Inf > and many NaNs): This is quite valuable, and will be stuffed away in my working directory for future references. Once again, thanks for nudging me in the right direction. One of these days, I'll have some of the long double and complex.h functions beaten into submission. -- Steve