From owner-freebsd-standards@FreeBSD.ORG Wed Oct 5 10:16:08 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 160AB16A41F for ; Wed, 5 Oct 2005 10:16:08 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailout2.pacific.net.au (mailout2.pacific.net.au [61.8.0.115]) by mx1.FreeBSD.org (Postfix) with ESMTP id 6B7F943D46 for ; Wed, 5 Oct 2005 10:16:07 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy2.pacific.net.au (mailproxy2.pacific.net.au [61.8.0.87]) by mailout2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j95AG5GZ027488; Wed, 5 Oct 2005 20:16:05 +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 j95AG35b005146; Wed, 5 Oct 2005 20:16:03 +1000 Date: Wed, 5 Oct 2005 20:16:03 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20051005032400.GA6736@troutmask.apl.washington.edu> Message-ID: <20051005191936.N51036@delplex.bde.org> 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> <20051005032400.GA6736@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed 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 10:16:08 -0000 On Tue, 4 Oct 2005, Steve Kargl wrote: > On Tue, Oct 04, 2005 at 07:31:43PM +1000, Bruce Evans wrote: >> On Sun, 2 Oct 2005, Steve Kargl wrote: > >> I converted to ccoshf() for testing it and comparing with a simpler >> implementation that uses the formula, then converted back. There are > So which implementation do you prefer: your simpler version or > your fixed version of my first cut at ccosh? If course I prefer the simpler version :-). I think version that handles all the exceptional cases explicitly is useful mainly for bootstrapping. Once things work the exceptional cases are better documented by regression checks for them. I looked at ccosh() in glibc and cephes. The generic version in glibc-2.3.2 is similar to your version except, it uses fpclassify() instead of bit fiddling, feraisexcept (FE_INVALID) instead of invalid operations, and many gccisms. It has the same bugs for ccosh(x+Iy) where x is finite, cosh(x) is infinite and y is 0. glibc doesn't have ccosh() in asm except on m68k's, but it has cexp() in asm on i386's; for that it uses far too much code to handle the exceptional case and has the same bug when cosh(x) is infinite, etc. cephes has all c9x complex functions and some others including cgamma() and zetac(), but has almost no support for exceptional values -- for ccosh() it just uses the formula. >> 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? Depends if isnan() is significantly inefficient. The classification macros are used a lot in your version so I think it is worth avoiding them, but in my version they are only used in rare cases. >> 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. > > :-) The main gccisms in the glibc version are near these assignments. `creal(f) = x' is spelled `__real__ f = x'. This spelling is also used in rvalues. `I' is never used, even in constants. I found one bug, and it was for the above case in both version. The above case is x = 0, y finite; the sign is that of sin(y) but the above assumes that it is that of y (i.e., 1 since we normalized y). One corrected version of the above is: cimag(f) = x * con * copysign(0, sin(y)); but it is simpler to delete the code for this special case and just use the formula. The x finite, y = 0 case is quite different since cosh(x) can overflow so the formula doesn't work, and the sign of sinh(x) is that of x since sinh() doesn't oscillate like sin() so we don't need a copysign() or an sinh() to determine the sign. >> 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? I think it should go in math_private.h only and be an inline function. Here is the handling of the exceptional-exceptional cases (ones which are not handled by the formula) for the simpler version after fixing the (x finite, y finite) case: % if (x == 0 && !isfinite(y)) { % creal(f) = y - y; % cimag(f) = copysign(0, y - y); % return (f); % } % if (y == 0) { % creal(f) = cosh(x); % cimag(f) = isnan(x) ? copysign(0, x - x) : % copysign(0, x) * y; % return (f); % } Bruce