[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: complex.h math functions
On Mon, 31 Oct 2005, Bruce Evans wrote:
On Sun, 30 Oct 2005, Steve Kargl wrote:
I've re-arrange the code to do x = 0 and then y = 0 cases
after the nearly-non-exceptional case. You may want to
flip x = 0 and y = 0. I did not do this because I did not
want to mess up your explanation of choice of signs.
I'll merge with it and see if there is more to clean up.
I cleaned it up not long after writing the above, but didn't get
around to replying until now. The folowing patches are supposed
to be relative to the latest versions in this thread.
% diff -c2 s_ccosh.c~ s_ccosh.c
% *** s_ccosh.c~ Mon Oct 31 21:07:23 2005
% --- s_ccosh.c Mon Oct 31 22:15:16 2005
% ***************
% *** 79,86 ****
% return (cpack(y - y, copysign(0, x * (y - y))));
%
% ! /*
% * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
% *
% ! * cosh(NaN + I 0) = NaN +- I 0.
% * The sign of 0 in the result is unspecified.
% */
% --- 79,86 ----
% return (cpack(y - y, copysign(0, x * (y - y))));
%
% ! /*
% * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
% *
% ! * cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0.
% * The sign of 0 in the result is unspecified.
% */
Whitespace, and fix a lost sign. d(NaN) is the default unary result
with a NaN arg; on i386's it is the same as the NaN arg for a quiet
NaN and the NaN with the same bits except for the quiet bit for a
signaling NaN: d(NaN) = NaN | QUIETBIT in bits. d(NaN, +-0) is
the corresponding default binary result. The quiet bit always gets
set when one arg is a NaN. IIRC, the result is the largest (in binary)
NaN if 2 NaNs are involved, and the sign bit only affects the order.
I forget if negative NaNs are larger.
% ***************
% *** 97,101 ****
% * cosh(x + I NaN) = d(NaN) + I d(NaN).
% * Optionally raises the invalid floating-point exception for finite
% ! * nonzero x. Choice = raise.
% */
% if (ix < 0x7ff00000 && iy >= 0x7ff00000)
% --- 97,101 ----
% * cosh(x + I NaN) = d(NaN) + I d(NaN).
% * Optionally raises the invalid floating-point exception for finite
% ! * nonzero x. Choice = don't raise (except for signaling NaNs).
% */
% if (ix < 0x7ff00000 && iy >= 0x7ff00000)
The comment was wrong. All these comments need to be checked.
% ***************
% *** 117,126 ****
% }
%
% ! /*
% ! * cosh(NaN + I NaN) = NaN + I NaN.
% *
% ! * cosh(NaN + I y) = NaN + I NaN.
% ! * Raise the invalid floating-point exception for finite nonzero y.
% *
% */
% return (cpack((x * x) * (y - y), (x + x) * (y - y)));
% --- 117,130 ----
% }
%
% ! /*
% ! * cosh(NaN + I NaN) = d(NaN) + I d(NaN).
% *
% ! * cosh(NaN +- I Inf) = d(NaN) + I d(NaN).
% ! * Optionally raises the invalid floating-point exception.
% ! * Choice = raise.
% *
% + * cosh(NaN + I y) = d(NaN) + I d(NaN).
% + * Optionally raises the invalid floating-point exception for finite
% + * nonzero y. Choice = don't raise (except for signaling NaNs).
% */
% return (cpack((x * x) * (y - y), (x + x) * (y - y)));
Whitespace, and +- and d() and backwards comment as above, and restore (?)
comment for the +- I Inf case.
% diff -c2 s_ccoshf.c~ s_ccoshf.c
% *** ss_ccoshf.c~ Wed Oct 12 05:10:40 2005
% --- ss_ccoshf.c Mon Oct 17 12:26:13 2005
% ***************
% *** 26,39 ****
%
% /*
% ! * Hyperbolic cosine of a complex argument z = x + i y such that i = sqrt(-1).
% */
%
% #include <complex.h>
%
% - double complex ccosh(double complex);
% -
% float complex
% ccoshf(float complex z)
% {
% ! return ((float complex) ccosh((double complex) z));
% }
% --- 26,38 ----
%
% /*
% ! * Hyperbolic cosine of a float complex argument.
% */
%
% #include <complex.h>
%
% float complex
% ccoshf(float complex z)
% {
% !
% ! return (ccosh(z));
% }
Style fixes.
The main thing missing is correct multiplication of cosh(x)*cos(y),
etc., when cosh(x) is Inf but the result should be finite. I have
made some progress fixing rounding errors in the corresponding
almost-correct multiplication for cosh(x) ~= exp(x)*0.5 when exp(x)
is Inf but the result is finite. fdlibm uses cosh(x) ~= exp(x/2)*exp(x/2)*0.5
when exp(x) is Inf but the result should be finite, but this gives
errors of > 1 ulp (up to 1 for each exp term and another 0.5 for the
multiplication, so up to 2.5 altogether in theory, and in practice
2.5029 for both coshf() and sinhf()). My fix for sinhf() uses sinhf(x)
~= P1(x - 88.875) when |x - 88.875| <= 0.25 and sinhf(x) ~= P2(x - 89.125)
when |x - 89.125| < 0.25, where P1 and P2 are similar to the
Taylor expansions about 88.875 and 89.125, respectively. This method
doesn't work for cosh(x)*cos(y) but was easy to implement for sinhf()
once I had similar approximations for |x| < 4 in many subintervals.
Bruce
