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>