Date: Wed, 5 Oct 2005 20:16:03 +1000 (EST) From: Bruce Evans <bde@zeta.org.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: freebsd-standards@freebsd.org Subject: Re: complex.h math functions Message-ID: <20051005191936.N51036@delplex.bde.org> In-Reply-To: <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> <20051005032400.GA6736@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051005191936.N51036>