Date: Fri, 7 Oct 2005 12:02:39 -0700 From: Steve Kargl <sgk@troutmask.apl.washington.edu> To: Bruce Evans <bde@zeta.org.au> Cc: freebsd-standards@FreeBSD.org Subject: Re: complex.h math functions Message-ID: <20051007190239.GA78674@troutmask.apl.washington.edu> In-Reply-To: <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> <20051007231439.F58005@delplex.bde.org>
next in thread | previous in thread | raw e-mail | index | archive | help
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 <complex.h>, 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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051007190239.GA78674>