Skip site navigation (1)Skip section navigation (2)
Date:      Sat, 8 Oct 2005 06:35:26 +1000 (EST)
From:      Bruce Evans <bde@zeta.org.au>
To:        Steve Kargl <sgk@troutmask.apl.washington.edu>
Cc:        freebsd-standards@freebsd.org
Subject:   Re: complex.h math functions
Message-ID:  <20051008052850.S59139@delplex.bde.org>
In-Reply-To: <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> <20051007190239.GA78674@troutmask.apl.washington.edu>

next in thread | previous in thread | raw e-mail | index | archive | help
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



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051008052850.S59139>