From owner-freebsd-standards@FreeBSD.ORG Wed Oct 12 07:23:35 2005 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 24EB716A41F for ; Wed, 12 Oct 2005 07:23:35 +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 5B60443D45 for ; Wed, 12 Oct 2005 07:23:34 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy1.pacific.net.au (mailproxy1.pacific.net.au [61.8.0.86]) by mailout2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j9C7NWZp028238; Wed, 12 Oct 2005 17:23:32 +1000 Received: from katana.zip.com.au (katana.zip.com.au [61.8.7.246]) by mailproxy1.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j9C7NRbd031939; Wed, 12 Oct 2005 17:23:29 +1000 Date: Wed, 12 Oct 2005 17:23:28 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20051010185153.GA55589@troutmask.apl.washington.edu> Message-ID: <20051012160109.I73531@delplex.bde.org> References: <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> <20051007231439.F58005@delplex.bde.org> <20051007190239.GA78674@troutmask.apl.washington.edu> <20051008052850.S59139@delplex.bde.org> <20051010185153.GA55589@troutmask.apl.washington.edu> 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: Wed, 12 Oct 2005 07:23:35 -0000 On Mon, 10 Oct 2005, Steve Kargl wrote: > (Attributions are a mess :-) Now accumulating them again :-). > I fixed this in my version of ccosh and used cpack(,) for > constructing the complex values. I tried to use cpack(,) > in your simplified version, but screwed something up in > that the return value is the complex conjugate from that > version. > > I've attached the latest version. Hopefully, I caught the > rest of the style(9) problems. The cosh.3 man page has been No :-). > updated. I did not hook s_ccosh.c up to the Makefile > because I wasn't sure were to put it in msun/src. If you have > no further comments, can you commit this version (or your > simplified version)? I'll look at the other trigometric and > hyperbolic functions when we converge on s_ccosh.c. I'm beginning to prefer the unsimplified version, but it needs some more editing after the following changes. The simplified version needed another case for ccosh(+-Inf + I NaN) and that already gives it about half as many cases and the other version. It saves code mainly by falling through to the code that works for finite args, but eventually the case of finite args will need to be more special to avoid overflows and then the non-finite cases won't be able to just use it. Listing all the special cases also serves as documentation. I think the current order of special cases is not quite the best, however. Patch relative to your last version: % --- s_ccosh.c~ Wed Oct 12 05:11:13 2005 % +++ s_ccosh.c Wed Oct 12 09:05:44 2005 % @@ -43,5 +43,5 @@ % ccosh(double complex z) % { % - double con, x, y; % + double x, y; % int32_t hx, hy, ix, iy, lx, ly; % % @@ -55,25 +55,10 @@ % iy = 0x7fffffff & hy; /* If iy >= 0x7ff00000, then inf or NaN. */ % % - /* Even function: ccosh(-z) = ccosh(z). */ % - if (hx & 0x80000000 && !isnan(x)) { % - x = -x; % - hx = ix; % - if (!isnan(y)) { % - y = -y; % - hy ^= 0x80000000; % - } % - } % - % - /* Required identity: ccosh(conj(z)) = conj(ccosh(z)). */ % - con = 1; % - if (hy & 0x80000000 && !isnan(y)) { % - con = -1; % - y = -y; % - hy = iy; % - } % - Delete the explicit sign flipping since it mostly just complicates things. % /* Handle the nearly-non-exceptional cases where x and y are finite. */ % - if (ix < 0x7ff00000 && iy < 0x7ff00000) % - return (cpack(cosh(x) * cos(y), con * sinh(x) * sin(y))); % + if (ix < 0x7ff00000 && iy < 0x7ff00000) { % + if ((iy | ly) == 0) % + return (cpack(cosh(x), x * y)); This special case was lost. Without this, cosh(x) * +-0 gives NaN if cosh(x) is +-Inf. x * y for the imaginary part is an efficient way to get +-0 with the right sign. % + return (cpack(cosh(x) * cos(y), sinh(x) * sin(y))); % + } % % /* % @@ -81,16 +66,7 @@ % * invalid exception. % * cosh(+0 + I NaN) = NaN +- I 0, sign of 0 is unspecified. % - * % - * Use the sign of (sinh(original_x) * sin(original_y)) for % - * the optional sign. This expression reduces to (y - y) for % - * all cases checked (mainly i386's). The fdlibm sin(y) is % - * certainly (y - y), but in the +-Inf cases we may have % - * flipped the sign of original_y to get y, and we want to % - * multiply by the sinh() term. These complications have no % - * effect with normal NaN semantics (-Inf - -Inf) = same NaN % - * as (Inf - Inf), and (finite * NaN) = same NaN). % */ Delete this since it became out of date (x and y haven't changed) and getting the sign right has fine details not just here. I don't want comments in the code about all of them. % if ((ix | lx) == 0 && iy >= 0x7ff00000) % - return (cpack(y - y, copysign(0, y - y))); % + return (cpack(y - y, copysign(0, x * (y - y)))); % % /* This change is a no-op on all machines that I know of. It is to get as close as possible to what the hardware would do for sinh(x) * sinh(y). x is +-0 so sinh(x) is the same as x. sinh(y) is (y - y) (a NaN). We can't just muliply these since C99 specifies a result of +-0, so we get as close as possible by taking the sign of the multiple. Note that there is no additional invalid exception for this -- there is already one for the real part. (Actually, there probably is a second one since evaluating (y - y) only once would be an invalid optimization, but these exceptions coalesce.) Multiplying by x has no effect on all machines that I know of. Normally, if y is Inf then (y - y) is some NaN and multiplying by x preserves this, and if y is NaN then (y - y) is the same NaN and that is preseved. % @@ -99,5 +75,5 @@ % */ % if (ix < 0x7ff00000 && iy >= 0x7ff00000) % - return (cpack(y - y, y - y)); % + return (cpack(y - y, x * (y - y))); % /* % * cosh(+Inf + I 0) = +Inf + I 0 Similarly, except here we want a NaN imaginary part so we don't need a copysign(). Cases like this can currently all use the formula. % @@ -108,8 +84,8 @@ % if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { % if ((iy | ly) == 0) % - return (cpack(x * x, con * y)); % - else if (iy >= 0x7ff00000) % - return (cpack(x * x, y - y)); % - return (cpack((x * x) * cos(y), (x * x) * con * sin(y))); % + return (cpack(x * x, copysign(0, x) * y)); This is the only place that deleting the sign flipping made more complicated later. The sign that was in "con" now needs to be determined. % + if (iy >= 0x7ff00000) Minor style fix (don't use "return; else"). % + return (cpack(x * x, x * (y - y))); Let the sign of x possibly affect the NaN (y - y) as usual. % + return (cpack((x * x) * cos(y), x * sin(y))); "con" just goes away (it is in y), provided we use the right sign for the x-related term here. (x * x) was logically wrong -- x is +-Inf and the term is sinh(x) = +-Inf = x. (x * x) only worked because we normalized x so it was +Inf. % } % /* % @@ -119,5 +95,5 @@ % */ % if ((iy | ly) == 0) % - return (cpack(x - x, copysign(0, x - x))); % - return (cpack((x - x) + (y - y), (x - x) + (y - y))); % + return (cpack(x * x, copysign(0, (x + x) * y))); Now x is NaN. (x - x) was logically wrong for cosh(x) -- cosh(x) is (x * x) for both Infs and NaNs -- it's a square so that it works for +-Inf, and this gives a choice of sign and value for NaNs too, depending on what the hardware does for NaN * NaN. Just use the same value here. Similarly for the imaginary part: - sinh(x) is (x + x), not (x - x), for Infs and NaNs, so as to preserve the sign for Infs. Just use the same value. - let the sign of y possibly affect the result like the sign of x did in other cases. % + return (cpack((x * x) * (y - y), (x + x) * (y - y))); % } Similarly. We want to let the sign of sin(y) possibly affect the real part. sin(y) is (y - y), so just use that. For the imaginary part, multiply instead of add to do the right thing with signs; just use sinh(x) = (x + x) not (x - x); sin(y) = (y - y) was already correct. % diff -ruN /usr/src/lib/msun/man/cosh.3 freebsd/src/lib/msun/man/cosh.3 % --- /usr/src/lib/msun/man/cosh.3 Fri Jan 14 15:28:28 2005 % +++ freebsd/src/lib/msun/man/cosh.3 Sat Oct 8 11:46:29 2005 % ... OK; could be more detailed. % diff -ruN /usr/src/lib/msun/src/complex/s_ccosh.c freebsd/src/lib/msun/src/complex/s_ccosh.c % --- /usr/src/lib/msun/src/complex/s_ccosh.c Wed Dec 31 16:00:00 1969 % +++ freebsd/src/lib/msun/src/complex/s_ccosh.c Mon Oct 10 11:10:14 2005 Patch was relative to this. % + % +/* % + * Hyperbolic cosine of a complex argument z = x + i y such that i = sqrt(-1). Other changes: don't define "i" here... % + ... % + ix = 0x7fffffff & hx; /* If ix >= 0x7ff00000, then inf or NaN. */ % + iy = 0x7fffffff & hy; /* If iy >= 0x7ff00000, then inf or NaN. */ % + % + /* Even function: ccosh(-z) = ccosh(z). */ % + if (hx & 0x80000000 && !isnan(x)) { % + x = -x; Still has lots of strange indents -- corrupt tab above, 4-chars only for most second indents. % + /* % + * cosh(x + I Inf) = NaN + I NaN, finite x != 0, invalid exception. % + * cosh(x + I NaN) = NaN + I NaN, finite x != 0, opt. inval. except. % + */ Not quite enough space to describe it all on 1 line; the inputs and outputs get mixed cryptically; might use more formal abbreviations. % + * cosh(+Inf + I y) = +Inf [cos(y) + I sin(y)], y != 0 and finite. Still have the strange brackets here. % diff -ruN /usr/src/lib/msun/src/complex/s_ccoshf.c freebsd/src/lib/msun/src/complex/s_ccoshf.c % --- /usr/src/lib/msun/src/complex/s_ccoshf.c Wed Dec 31 16:00:00 1969 % +++ freebsd/src/lib/msun/src/complex/s_ccoshf.c Fri Oct 7 17:30:36 2005 % ... % +float complex % +ccoshf(float complex z) % +{ % + return ((float complex) ccosh((double complex) z)); % +} I prefer to use implicit conversions. % diff -ruN /usr/src/lib/msun/src/math_private.h freebsd/src/lib/msun/src/math_private.h % --- /usr/src/lib/msun/src/math_private.h Fri Feb 4 12:05:39 2005 % +++ freebsd/src/lib/msun/src/math_private.h Fri Oct 7 17:30:40 2005 % @@ -155,6 +155,36 @@ % } while (0) % % /* % + * Inline functions that can be used to construct complex values. % + */ I moved my comment about the gcc bug -- the reason for existence of these functions -- to here. % +static __inline float complex % +cpackf(float __x, float __y) I needed a namespace hack to make this compile. Clients that don't use this shouldn't need to include , and this file shouldn't include it. I used "#ifdef I". I think the underscores shouldn't be used here (not even for __inline). This is not a public interface so we don't need to be very careful with the namespace. Otherwise good. These interfaces seem to work well. Here is my current "simpler" version. I plan to merge it into the other version when the other version's indentationis fixed. %%% /* * Hyperbolic cosine of a double complex argument. * * Most exceptional values are automatically correctly handled by the * standard formula. See n1124.pdf for details. */ #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. */ 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))); } Here is my program for comparing complex functions. It has been cleaned up a bit (it now uses fixed-width types but not GET_FLOAT_WORD() etc. to access them properly). Compile with -DTESTD for simplistic double precision comparison. Oops, it needs to use cpack*(). %%% #include #include #include #include #include #ifndef FUNC #define FUNC ccosh #endif double complex ccosh(double complex z); double complex ccosh1(double complex z); float complex ccoshf(float complex z); float complex ccosh1f(float complex z); int main(void) { float complex vf0, vf1; float complex z; volatile float vf0i, vf0r, vf1i, vf1r; float x, y; uintmax_t na, tot_er; uint32_t eri, err, max_er, stride, vf0iu, vf0ru, vf1iu, vf1ru, xu, yu; max_er = 0; na = 0; stride = 0x80000; tot_er = 0; for (xu = 0;; xu += stride) { for (yu = 0;; yu += stride) { x = *(float *)&xu; y = *(float *)&yu; /* `z = x + I * y;', unbroken: */ crealf(z) = x; cimagf(z) = y; #ifdef TESTD /* * Hack for double precision functions. Testing them only * on float args is OK, but we should report errors in * double precision for them. */ vf0 = __CONCAT(FUNC, )(z); vf1 = __CONCAT(FUNC, 1)(z); #else vf0 = __CONCAT(FUNC, f)(z); vf1 = __CONCAT(FUNC, 1f)(z); #endif vf0r = crealf(vf0); vf0ru = *(volatile uint32_t *)&vf0r; vf0i = cimagf(vf0); vf0iu = *(volatile uint32_t *)&vf0i; vf1r = crealf(vf1); vf1ru = *(volatile uint32_t *)&vf1r; vf1i = cimagf(vf1); vf1iu = *(volatile uint32_t *)&vf1i; na++; err = abs(vf1ru - vf0ru); eri = abs(vf1iu - vf0iu); tot_er += err + eri; if (max_er < err) max_er = err; if (max_er < eri) max_er = eri; if (err > 4 || eri > 4) { printf("z = %#10x+I*%-#10x %.15g+I*%-.15g\n", xu, yu, x, y); printf( __XSTRING(FUNC) "f(z) = %#10x+I*%-#10x %.15g+I*%-.15g\n", vf0ru, vf0iu, vf0r, vf0i); printf( __XSTRING(FUNC) "1f(z) = %#10x+I*%-#10x %.15g+I*%-.15g\n", vf1ru, vf1iu, vf1r, vf1i); printf("err = %s%#x+I*%s.%#x\n", vf0ru <= vf1ru ? "+" : "-", err, vf0iu <= vf1iu ? "+" : "-", eri); fflush(stdout); } if (na % 10000000 == 1) { printf( "%ju: %#x+I*%#x: %g+I*%g: max_error = %#x, avg_error = %.3f\n", na, xu, yu, x, y, max_er, (double)tot_er / na); fflush(stdout); } if (yu == -stride) break; } if (xu == -stride) break; } printf("%ju cases: max_error = %#x, avg_error = %.3f\n", na, max_er, (double)tot_er / na); } #include "s_ccosh.c" #include "s_ccosh1.c" #include "s_ccoshf.c" #include "s_ccosh1f.c" /* For ccosh(), must avoid using the broken hardware cos() and sin(): */ #undef rcsid #define rcsid rcsid6 #include "../s_cos.c" #undef rcsid #define rcsid rcsid7 #include "../s_sin.c" #ifdef DEBUG #undef rcsid #define rcsid rcsid0 #include "../e_coshf.c" #undef one #define one one0 #undef rcsid #define rcsid rcsid1 #include "../e_sinhf.c" #undef one #define one one1 #undef rcsid #define rcsid rcsid2 #include "../s_tanhf.c" #endif %%% Bruce