complex.h math functions

[Steve Kargl]

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