From owner-freebsd-standards@FreeBSD.ORG Fri Oct 7 14:51:27 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 56BB716A424 for ; Fri, 7 Oct 2005 14:51:27 +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 5CB3843D48 for ; Fri, 7 Oct 2005 14:51:26 +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 j97EpOmi031653; Sat, 8 Oct 2005 00:51:24 +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 j97EpMWX028252; Sat, 8 Oct 2005 00:51:22 +1000 Date: Sat, 8 Oct 2005 00:51:22 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20051006212602.GA40609@troutmask.apl.washington.edu> Message-ID: <20051007231439.F58005@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> <20051005191936.N51036@delplex.bde.org> <20051006212602.GA40609@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: Fri, 07 Oct 2005 14:51:27 -0000 On Thu, 6 Oct 2005, Steve Kargl wrote: > On Wed, Oct 05, 2005 at 08:16:03PM +1000, Bruce Evans wrote: >> If course I prefer the simpler version :-). I think version that handles >> all the exceptional cases explicitly is useful mainly for bootstrapping. > > I've implemented some inline functions (see below) for the assignments of > the real and imaginary parts. After replacing > > creal(f) = cosh(x) * cos(y); > cimag(f) = sinh(x) * sin(y); > return (f); > > in your simpler function by > > return (cmplx(cosh(x) * cos(y), sinh(x) * sin(y)); > > your test program is generating a large number of > > z = 0x7f800000+I0x7fa00000 inf+Inan > ccosh(z) = 0x7f800000+I0x7fe00000 inf+Inan > ccosh1(z) = 0x7fe00000+I0x7fe00000 nan+Inan > err = +0x600000+I+.0 > > where ccosh1 is your simpler function and ccosh is your fixed > version of my ccosh (both have been converted to use the inline > functions). n869.pdf says > > ccosh(inf+Inan) = inf+Inan This is because: (1) n869.pdf is broken here. For +inf*nan, if we knew that the nan meant an indeterminate but positive number (possibly including inf), then +inf*nan should be +inf, but in both of the expressions cosh(+inf) * cos(nan) and sinh(+inf) * sin(nan) we don't know the signs of the trig terms so the result should be nan for both the real and the imaginary parts. (2) s_ccosh.c still uses essentially your version for the return value in this case. It is "x * x + I * (y - y)". x * x is supposed to give inf, but the gcc bug propagates the nan y * y into the real part. (3) s_ccosh1.c is bug for bug compatible with s_ccosh.c here. Passing the regression test forced it to be incompatible. (4) using cmplx(u, v) instead of u + I*v fixed s_ccosh.c, and the regression test shows the difference. Unless this is fixed in C99[+TC], s_ccosh1.c needs to be broken to be compatible. The C99 behaviour has implications for cproj(). cproj(inf+Inan) is inf+Icopysign(0,nan), but cproj(nan+Inan) is nan+Inan. I think it is a larger bug to completely lose the nan-ness here since the real part is not really inf. BTW, many spurious infs and maybe some spurious nans occur when the terms overflow but the result doesn't. E.g., cosh(x) might be inf so cosh(x) * cos(x) is +-inf for all nonzero x, but if cos(x) is say 0.5 then overflow shouldn't occur until slightly larger x. I recently noticed the great lengths to which fdlibm goes to avoid similar overflows for the real functions. E.g., for cosh(), cosh(x) ~= exp(x)/2 and the RHS of this would overflow for x slightly less large than for the LHS. So instead, for large x fdlibm evaluates the RHS as exp(x/2)/2.*exp(x/2). This is imperfect due to extra roundoff error from the second multiplication. In practice, it gives errors of up to ~1.5 ulps in cosh() and then errors of up ~3 ulps in ccosh1(). So an extra-precision exp() is apparently needed to get 1-ulp accuracy even for the real cosh(). ccosh() could creep up on oveflow similarly. It would need about 53+[1 or 2] bits of extra precision to handle zero points of the trig terms where cosh() only nneds 1 or 2 bits extra to handle the 1./2 term. I don't think this will be implemented soon. >>>> 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.) > > Apparently, my version of gcc does not like the above > code. How did you compile your s_ccosh1.c? I used gcc-3.3.3. The problem is that creal(f) never was an lvalue since it is declared as a function in , and gcc has finally fixed its its default of violating the C standard's requirements for lvalues. The uglier syntax used in glibc (e.g., __real__ f = cos(y)) still works. > --- /mnt1/src/lib/msun/src/math_private.h Fri Feb 4 12:05:39 2005 > +++ ../math_private.h Thu Oct 6 14:06:37 2005 > @@ -155,6 +155,43 @@ > } while (0) > > /* > + * Inline functions that can be used to construct complex values. > + */ > + > +static __inline float complex > +cmplxf(float __x, float __y) > +{ > + float *__ptr; > + float complex __z; > + __ptr = (float *) &__z; > + *__ptr++ = __x; > + *__ptr = __y; > + return (__z); > +} Use __real__/__imag or something less ugly if possible. These and creal()/cimag() are mentioned in gcc.info but neither is really documented there. Neither is mentioned to work as an lvalue. But the layout of "float complex" needed for the pointer hacks in the above to work is documented to not always work there -- it is documented that the the floats may noncontiguous, even with 1 in memory and the other in a register. Please think of better names. I think "complex" should always be abbreviated to "c" in libm. Maybe cpackf()? Bruce