Date: Sat, 8 Oct 2005 00:51:22 +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: <20051007231439.F58005@delplex.bde.org> In-Reply-To: <20051006212602.GA40609@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>
next in thread | previous in thread | raw e-mail | index | archive | help
On Thu, 6 Oct 2005, Steve Kargl wrote:
> On Wed, Oct 05, 2005 at 08:16:03PM +1000, Bruce Evans wrote:
>> If course I prefer the simpler version :-). I think version that handles
>> all the exceptional cases explicitly is useful mainly for bootstrapping.
>
> I've implemented some inline functions (see below) for the assignments of
> the real and imaginary parts. After replacing
>
> creal(f) = cosh(x) * cos(y);
> cimag(f) = sinh(x) * sin(y);
> return (f);
>
> in your simpler function by
>
> return (cmplx(cosh(x) * cos(y), sinh(x) * sin(y));
>
> 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.
(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.
(3) s_ccosh1.c is bug for bug compatible with s_ccosh.c here. Passing
the regression test forced it to be incompatible.
(4) using cmplx(u, v) instead of u + I*v fixed s_ccosh.c, and the regression
test shows the difference.
Unless this is fixed in C99[+TC], s_ccosh1.c needs to be broken to be
compatible.
The C99 behaviour has implications for cproj(). cproj(inf+Inan) is
inf+Icopysign(0,nan), but cproj(nan+Inan) is nan+Inan. I think it is
a larger bug to completely lose the nan-ness here since the real part is
not really inf.
BTW, many spurious infs and maybe some spurious nans occur when the
terms overflow but the result doesn't. E.g., cosh(x) might be inf so
cosh(x) * cos(x) is +-inf for all nonzero x, but if cos(x) is say 0.5
then overflow shouldn't occur until slightly larger x. I recently
noticed the great lengths to which fdlibm goes to avoid similar overflows
for the real functions. E.g., for cosh(), cosh(x) ~= exp(x)/2 and the
RHS of this would overflow for x slightly less large than for the LHS.
So instead, for large x fdlibm evaluates the RHS as exp(x/2)/2.*exp(x/2).
This is imperfect due to extra roundoff error from the second
multiplication. In practice, it gives errors of up to ~1.5 ulps in
cosh() and then errors of up ~3 ulps in ccosh1(). So an extra-precision
exp() is apparently needed to get 1-ulp accuracy even for the real
cosh(). ccosh() could creep up on oveflow similarly. It would need
about 53+[1 or 2] bits of extra precision to handle zero points of the
trig terms where cosh() only nneds 1 or 2 bits extra to handle the 1./2
term. I don't think this will be implemented soon.
>>>> 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;
>>>> % + 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.
> --- /mnt1/src/lib/msun/src/math_private.h Fri Feb 4 12:05:39 2005
> +++ ../math_private.h Thu Oct 6 14:06:37 2005
> @@ -155,6 +155,43 @@
> } while (0)
>
> /*
> + * Inline functions that can be used to construct complex values.
> + */
> +
> +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.
Please think of better names. I think "complex" should always be
abbreviated to "c" in libm. Maybe cpackf()?
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20051007231439.F58005>