complex.h math functions

[Bruce Evans]

On Fri, 7 Oct 2005, Steve Kargl wrote:

> On Sat, Oct 08, 2005 at 12:51:22AM +1000, Bruce Evans wrote:
>> On Thu, 6 Oct 2005, Steve Kargl wrote:
>>> ... 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.

I just checked that too.

The rationale (both the 1998 n850.pdf one and the new one near n1124.pdf)
says a bit about this.  Infinities are supposed to be preserved as much as
possible.  In particular, it is almost right for Inf * NaN to be Inf,
not NaN, since NaN is best thought of as an indeterminate number and
Inf * non_NaN is +-Inf except when non_NaN is 0.  There is a problem with
the sign if the resulting Inf being indeterminate, but that vanishes if
we use projective infinity.  However the problem with Inf * 0 is too large,
so both Inf * 0 and Inf * NaN are NaN for IEEE floating point.

The rationale for ccosh(inf+Inan) = inf+Inan is partly that after
projectivizing the value is unambigously the unique projective infinity.
I think this is true for the finite precision case because cos(y) and
sin(y) are never (?) zero for real rational x != 0.  Thus the cosh(x)
and sinh(x) terms eventually grow large enough to give +-Inf when
multiplied by the cos(y) and sin(y) terms (calculating everything in
infinite precision until the final step of rounding to double precision
gives +-Inf).

There is a similar rationale for pow(-2.0, Inf) being +Inf (its sign
is uniquely determined since the domain is limited to even values if
the base is 2).  I like this rationale for pow() but not for ccosh().
Where pow(-2.0, y) increases monotonically to +Inf as y increases to
+Inf, ccosh(z) has an essentially singularity as z -> projective
infinity, so in infinite precision it takes on all values infinitely
many times and in finite precision it oscilates wildly.

>>>>>> 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.

Right.  I found this bug a few days ago and mentioned it in previous
mail.  I have only fixed it for the float versions.

>> 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.

I'm used to C, and forget what Fortran does.  I think of (x, y) -> x+I*y
as a constructor but couldn't think of a good (short) name with
"constructor" in it.

Other ideas that don't quite seem to work:
- only provide an imaginifying constructor from y to I*y.  I think this
   would give pessimizations from extra additions for addition x to it.
- use C++ and overload + and * in x+I*y with non-broken operators so as
   to keep the sources clean and not have to clean them up when x+I*y is
   fixed in gcc.  I would have to learn C++.

Bruce