From owner-freebsd-arch@FreeBSD.ORG Fri Aug 21 09:45:18 2009 Return-Path: Delivered-To: freebsd-arch@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 6E2AA106568B for ; Fri, 21 Aug 2009 09:45:18 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail04.syd.optusnet.com.au (mail04.syd.optusnet.com.au [211.29.132.185]) by mx1.freebsd.org (Postfix) with ESMTP id 08C2C8FC63 for ; Fri, 21 Aug 2009 09:45:17 +0000 (UTC) Received: from c122-106-152-1.carlnfd1.nsw.optusnet.com.au (c122-106-152-1.carlnfd1.nsw.optusnet.com.au [122.106.152.1]) by mail04.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id n7L9jDfQ032411 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Fri, 21 Aug 2009 19:45:14 +1000 Date: Fri, 21 Aug 2009 19:45:12 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20090820230400.GA61066@troutmask.apl.washington.edu> Message-ID: <20090821190006.E36897@delplex.bde.org> References: <20090815015947.GA4682@server.vk2pj.dyndns.org> <20090817201319.GA22437@zim.MIT.EDU> <20090820230400.GA61066@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-arch@freebsd.org Subject: Re: C99 Complex Math Functions X-BeenThere: freebsd-arch@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Discussion related to FreeBSD architecture List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 21 Aug 2009 09:45:18 -0000 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 % #include % % #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