From owner-freebsd-standards@FreeBSD.ORG Fri Oct 7 19:02: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 5DA2E16A41F for ; Fri, 7 Oct 2005 19:02: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 D3B1843D48 for ; Fri, 7 Oct 2005 19:02:39 +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 j97J2dJi056751; Fri, 7 Oct 2005 12:02: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 j97J2d0V056750; Fri, 7 Oct 2005 12:02:39 -0700 (PDT) (envelope-from sgk) Date: Fri, 7 Oct 2005 12:02:39 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20051007190239.GA78674@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> <20051005191936.N51036@delplex.bde.org> <20051006212602.GA40609@troutmask.apl.washington.edu> <20051007231439.F58005@delplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20051007231439.F58005@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, 07 Oct 2005 19:02:40 -0000 On Sat, Oct 08, 2005 at 12:51:22AM +1000, Bruce Evans wrote: > On Thu, 6 Oct 2005, Steve Kargl wrote: > > >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. It's still broken in n1124.pdf. > (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. I replaced the "x * x + I * (y - y)" with cmplx(x * x, y - y). > (4) using cmplx(u, v) instead of u + I*v fixed s_ccosh.c, and the regression > test shows the difference. Ah, this is the key. (ULP discussion of cosh(x) deleted). > >>>>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; I think the above is wrong. We need "x * con * sin(y)" because sin(y) is between -1 and 1. For example we may have x * con * sin(y) = (-0) * (1) * (-0.5) = +0 but I need to check n1124.pdf to see what the behavior of signed zero arithmetic does. > >>>>% + 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. > > >+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. Ugh. If I read Section 6.2.5(13), 6.2.5(20), and 6.2.6.1(2) correctly, then z in the above needs to occupy contiguous memory. I'll use the __real__/__imag__ gccism. > Please think of better names. I think "complex" should always be > abbreviated to "c" in libm. Maybe cpackf()? You don't do Fortran, do you? :-) program a complex z real :: x = 1., y = 2. z = cmplx(x,y) end program a OK, I'll use cpackf, cpack, and cpackl. -- Steve