Date: Fri, 21 Aug 2009 19:45:12 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: Steve Kargl <sgk@troutmask.apl.washington.edu> Cc: freebsd-arch@freebsd.org Subject: Re: C99 Complex Math Functions Message-ID: <20090821190006.E36897@delplex.bde.org> In-Reply-To: <20090820230400.GA61066@troutmask.apl.washington.edu> References: <20090815015947.GA4682@server.vk2pj.dyndns.org> <20090817201319.GA22437@zim.MIT.EDU> <20090820230400.GA61066@troutmask.apl.washington.edu>
next in thread | previous in thread | raw e-mail | index | archive | help
On Thu, 20 Aug 2009, Steve Kargl wrote:
> On Mon, Aug 17, 2009 at 04:13:19PM -0400, David Schultz wrote:
>
> In a private response to Peter, I directed him to NetBSD. NetBSD
> contains some/most(?)/all(?) of the simple formula implementations.
> I sent Peter the guts of ccosh as implemented by NetBSD (the simple
> formula) and the current version of ccosh that Bruce and I had written.
> The netbsd version is 15 lines (excluding the license). The bde+me
> version is 91 lines (excluding license and initial comment). To be
> fair, a large portion of those 91 lines are comments of the form
>
> /*
> * cosh(+-Inf + I NaN) = +Inf + I d(NaN).
> *
> * cosh(+-Inf +- I Inf) = +Inf + I dNaN.
> * The sign of Inf in the result is unspecified. Choice = always +.
> * Raise the invalid floating-point exception.
> *
> * cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y)
> */
> ...
We also have a ~15 line version for ccosh[f]() only (after removing
comments that are less detailed than the above). Howver, the error handling
has only been checked for the float version, and the overflow handling is
mostly nonexistent for all versions. Here is a simple version:
% /*
% * Hyperbolic cosine of a double complex argument.
% *
% * Most exceptional values are automatically correctly handled by the
% * standard formula. See n1124.pdf for details.
^^^^^^^^^
This (draft) C99+ standard gives details like the above, except of course
our choice of sign and our details of NaN handling).
% */
%
% #include <complex.h>
% #include <math.h>
%
% #include "../math_private.h"
%
% double complex
% ccosh1(double complex z)
% {
% double x, y;
%
% x = creal(z);
% y = cimag(z);
%
% /*
% * This is subtler than it looks. We handle the cases x == 0 and
% * y == 0 directly not for efficiency, but to avoid multiplications
% * that don't work like we need. In these cases, the result must
% * be almost pure real or a pure imaginary, except it has type
% * float complex and its impure part may be -0. Multiplication of
% * +-0 by an infinity or a NaN would give a NaN for the impure part,
% * and would raise an unwanted exception.
% *
% * We depend on cos(y) and sin(y) never being precisely +-0 except
% * when y is +-0 to prevent getting NaNs from other cases of
% * +-Inf*+-0. This is true in infinite precision (since pi is
% * irrational), and checking shows that it is also true after
% * rounding to float precision.
This is impossible to check in a reasonable time for double and higher
precisions, but our implementation of accuracy for trig functions already
depends on more (on certain bounding of the results away from 0).
% */
% if (x == 0 && !isfinite(y))
% return (cpack(y - y, copysign(0, x * (y - y))));
% if (y == 0)
% return (cpack(cosh(x), isnan(x) ? copysign(0, (x + x) * y) :
% copysign(0, x) * y));
% if (isinf(x) && !isfinite(y))
% return (cpack(x * x, x * (y - y)));
% return (cpack(cosh(x) * cos(y), sinh(x) * sin(y)));
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Standard formulas like this almost work for this and for all elementary
transcendental functions. There are obvious accuracy and overflow
problems in such formulas. Here cosh(x) may overflow although all of
the infinitely precise results are less than DBL_MAX. All these problems
can be reduced by evaluating expressions in extra precision, but that
would be slow, and no extra precision is availble for the long double
versions.
% }
Bruce
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20090821190006.E36897>