From owner-freebsd-standards@FreeBSD.ORG Tue Jan 3 08:13:29 2006 Return-Path: X-Original-To: freebsd-standards@freebsd.org Delivered-To: freebsd-standards@freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id CD17916A41F for ; Tue, 3 Jan 2006 08:13:29 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailout2.pacific.net.au (mailout2.pacific.net.au [61.8.0.115]) by mx1.FreeBSD.org (Postfix) with ESMTP id 082F243D49 for ; Tue, 3 Jan 2006 08:13:28 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy2.pacific.net.au (mailproxy2.pacific.net.au [61.8.0.87]) by mailout2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id k038DRuR024768; Tue, 3 Jan 2006 19:13:27 +1100 Received: from epsplex.bde.org (katana.zip.com.au [61.8.7.246]) by mailproxy2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id k038DOUW015026; Tue, 3 Jan 2006 19:13:25 +1100 Date: Tue, 3 Jan 2006 19:13:24 +1100 (EST) From: Bruce Evans X-X-Sender: bde@epsplex.bde.org To: Steve Kargl In-Reply-To: <20051031184647.Q38664@delplex.bde.org> Message-ID: <20060103183940.B4246@epsplex.bde.org> References: <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> <20051008052850.S59139@delplex.bde.org> <20051010185153.GA55589@troutmask.apl.washington.edu> <20051012160109.I73531@delplex.bde.org> <20051016184129.GA24651@troutmask.apl.washington.edu> <20051021194624.H598@epsplex.bde.org> <20051030182702.GA18998@troutmask.apl.washington.edu> <20051031184647.Q38664@delplex.bde.org> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@freebsd.org Subject: Re: complex.h math functions X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Tue, 03 Jan 2006 08:13:29 -0000 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 % % - 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 % % 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