From owner-freebsd-standards@FreeBSD.ORG Tue Oct 4 09:31:49 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 A6DA916A41F for ; Tue, 4 Oct 2005 09:31:49 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailout1.pacific.net.au (mailout1.pacific.net.au [61.8.0.84]) by mx1.FreeBSD.org (Postfix) with ESMTP id E7E0743D48 for ; Tue, 4 Oct 2005 09:31:48 +0000 (GMT) (envelope-from bde@zeta.org.au) Received: from mailproxy2.pacific.net.au (mailproxy2.pacific.net.au [61.8.0.87]) by mailout1.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j949Vldb020422; Tue, 4 Oct 2005 19:31:47 +1000 Received: from katana.zip.com.au (katana.zip.com.au [61.8.7.246]) by mailproxy2.pacific.net.au (8.13.4/8.13.4/Debian-3) with ESMTP id j949Vhl0012775; Tue, 4 Oct 2005 19:31:43 +1000 Date: Tue, 4 Oct 2005 19:31:43 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20051002191841.GA40568@troutmask.apl.washington.edu> Message-ID: <20051004164618.U47262@delplex.bde.org> References: <20050929195552.GA14982@troutmask.apl.washington.edu> <20050930221846.T34595@delplex.bde.org> <20050930152734.GA20655@troutmask.apl.washington.edu> <20051001204005.N37704@delplex.bde.org> <20051002191841.GA40568@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: Tue, 04 Oct 2005 09:31:49 -0000 On Sun, 2 Oct 2005, Steve Kargl wrote: >>> ... >>> OK, I'll use the implementation details of ccosh to infer the >>> handling of exceptional values. >>> > > At the risk of embarassing myself once again on freebsd-standards, > here's an implementation of ccosh. If this is acceptable, I'll > proceed with implementations for csinh, ctanh, ccos, csin, and > ctan. > > One question: Can we use a #define for complex float or do we > need an trivial C program? That is, in math.h can we do > > #define ccoshf(z) ccosh((float complex)(z)) A program is needed, since (ccoshf)(z) and &(ccoshf) must work even if ccoshf() is a macro. I converted to ccoshf() for testing it and comparing with a simpler implementation that uses the formula, then converted back. There are lots of special cases for -0, Infs and NaNs, and there were lots of bugs, including compiler bugs (gcc is very far from supporting the IEC 60559 extension). It turns out that the version using formula handles most of the special cases automatically and is much easier to fix. I fixed both versions and worked around the compiler bugs, and compared differences to find the best version of the variant using the formula. First, fixes for your version: % --- s_ccosh.c~ Mon Oct 3 07:10:59 2005 % +++ s_ccosh.c Tue Oct 4 16:24:06 2005 % @@ -40,11 +40,11 @@ % #include % % -#include "math_private.h" % +#include "../math_private.h" Temporary hack. % % double complex % ccosh(double complex z) % { % - int con = 1; % - double x, y; % + double complex f; % + double con, x, y; % int32_t hx, hy, ix, iy, lx, ly; % Optimize `con' a little, and fix some style bugs (initialization in declaration, and disordered declarations; most other style bugs are not fixed or noted). `f' is used to work around gcc bugs. % @@ -52,15 +52,4 @@ % y = cimag(z); % % - if (x < 0.) { % - /* Even function: ccosh(-z) = ccosh(z). */ % - x = -x; % - y = -y; % - /* Required identity: ccosh(conj(z)) = conj(ccosh(z)). */ % - if (y < 0.) { % - con = -1; % - y = -y; % - } % - } % - Move this to after extracting words so that we can look at the sign bit in ix and iy to classify -0. Note that n869.pdf only describes the behaviour explicitly for +0 and +Inf so that it doesn't have to describe many more cases in detail, but ut is clear that the evenness and conjugacy rules to the sign bit in +-0 although this isn't stated explicitly at least near the spec for ccosh(). Evenness is stated explicitly for cos(+-0) and oddness is stated explicitly for sin(+-0) (e.g., sin(-0) is specified to be -0, and this is implemented by fdlibm and i387 hardware sin). Although changing signs so that all the sign bits are 1 simplifies n869.pdf, it mostly complicates things for us since we can depend on multiplication to do the right things with signs. We need complications to avoid flipping the sign bit for NaNs so as to set the sign bits in the correct MD way for the cases where the standard doesn't specify the sign. Most operations on NaNs don't change the sign bit on most (?) machines, but -NaN changes it on i386's at least (0-NaN is different from -NaN on i386's!). % EXTRACT_WORDS(hx, lx, x); % EXTRACT_WORDS(hy, ly, y); % @@ -69,10 +58,37 @@ % iy = 0x7fffffff&hy; /* if iy >= 0x7ff00000, then inf or NaN */ % % - /* Filter out the non-exceptional cases. x and y are finite. */ This was changed to match the code, but didn't really move. % + /* 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; % + } The symmetry for conjugacy is independent of the one for evenness, so flipping the sign for it must be moved out of the if for flipping the sign(s) for evenness. Avoid flipping the sign bit for NaNs. I didn't bother optimizing isnan(x), etc. using hx, etc. % + % + /* Handle the nearly-non-exceptional cases where x and y are finite. */ Change comment to match code. % if (ix < 0x7ff00000 && iy < 0x7ff00000) { % - /* cosh(+0 + I 0) = 1. */ % - if ((ix|lx) == 0 && (iy|ly) == 0) % - return 1.; Adjust this to handle all cases where x is 0 (and y is finite). This fixes handling of cosh(+-0 + I*+-0) and handles nonzero y at negative extra cost. The result of cosh(+-0 + I*+-0) must be (1. + I*(+-0)*(+-0)). We have forgotten the original x and y, else we could return (1. + I*x*y) after working around the compiler bugs. (x * con) would give a correctly-signed 0, but it is simpler to use this in new code that handles all finite y (and zero x). % - return (cosh(x) * (cos(y) + I * con * tanh(x) * sin(y))); This changes but doesn't really move. % + if ((ix | lx) == 0) { Handle x == 0, y finite. There are 2 unobvious details: - getting the sign right when y == 0... % + creal(f) = cos(y); % + cimag(f) = x * con; % + return f; ... - working around gcc bugs. gcc's support for IEC 60559 is so noexistent that u + I*v is unusable in general to construct a double complex from 2 doubles. Here the main bug is that I*-0 gets corrupted to I*+0. gcc generates code that does extra work to implement this bug. gcc's handling of full complex operations is more reasonable -- it generates code that doesn't do the extra, very expensive work to avoid related bugs. See n869.pdf section G.4.1 for how difficult it is to implement complex multiplication perfectly correctly. The example there handles Infs and and NaNs and some overflows, but not underflows or signed 0's. The complications are related to the ones here and to the ones that I meantioned for x*x. E.g., the real part of (a+ib)(c+id) is ac-bd on real numbers, but for finite floating point numbers the multiplications can overflow without the their difference overflowing, and their difference can underflow; Infs, NaNs and signed 0's give additional complications. Fortunately, we can construct u + I*v by assigning to the real and complex parts. There should be a macro or inline function to do this. % + } % + if ((iy | ly) == 0) { % + creal(f) = cosh(x); % + cimag(f) = con * y; % + return f; % + } This special case for x finite and y zero is much more critical. cosh(x) can be infinite for finite x, and then we must not perform or use the multiplication (cosh(x) * tanh(x) * con * sin(y)), since that would be (Inf * +-0) so it would raise the invalid exception and give result NaN, but we want a result of +-0 for the imaginary part (just like for x infinite and y 0). We work around the compiler bug as usual. % + creal(f) = cosh(x) * cos(y); % + cimag(f) = cosh(x) * tanh(x) * con * sin(y); % + return f; % } % /* The only change for finite nonzero x and y is to avoid the compiler bug. When sinh(x) is Inf and y is nonzero, it is correct for the imaginary part to be NaN. The (cosh(x) * tanh(x)) term should be written as sinh(x), especially now that cosh(x) is evaluated twice. The old way of writing the return value: cosh(x) * (cos(y) + I * con * tanh(x) * sin(y)) has the interesting property of avoiding some of the compiler bugs. The second term is always finite, so gcc doesn't mess up Infs and Nans in it, and multiplying it by an infinite real coshx() does less damage than multiplying I by an infinite real (sinh(x) * con * siny()).. Using cosh(x) * tanh(x) instead of sinh(x) is also wrong because in theory cosh(x) might be Inf when sinh(x) is finite. I think this doesn't happen in practice. Similarly, sinh(x) might be infinite while the infinite precision sinh(x) * sin(y) is finite -- this is a subcase of the problem of multiplying 2 double complexes perfectly correctly. Regression tests showed that ccoshf() implemented using float precision differed significantly from ccosh() implemented using double precision only because the double precision version mostly overflow here (except when the result overflows a float). % @@ -81,6 +97,19 @@ % * cosh(+0 + I NaN) = NaN +- I 0, sign of 0 is unspecified. % */ % - if ((ix|lx) == 0 && iy >= 0x7ff00000) % - return (y - y) / (y - y); % + if ((ix|lx) == 0 && iy >= 0x7ff00000) { % + creal(f) = y - y; Work around compiler bugs as usual. y is Inf or Nan, so we can use the simpler formula y - y to generate a NaN for it. ((y - y) / (y - y)) is only needed to generate a NaN from a possibly-finite y. % + /* % + * 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). % + */ % + cimag(f) = copysign(0, y - y); % + return f; % + } % /* % * cosh(x + I Inf) = NaN + I NaN, finite x != 0, invalid exception. Intricacies to set the sign to the least surprising/most reproducible value. % @@ -88,5 +117,5 @@ % */ % if (ix < 0x7ff00000 && iy >= 0x7ff00000) % - return (x - x) / (x - x) + I * (x - x) / (x - x); % + return y - y + I * (y - y); % /* % * cosh(+Inf + I 0) = +Inf + I 0 Use the right arg in the return value so that NaNs get propagated to the result. Here x is finite (nonzero) and y is +-Inf or NaN. For y = +-Inf, y - y gives a standard fixed NaN "real indefinite" on i386's at least, the same one as the expression involving x, but for y = NaN we want to return the same NaN and not change to "real indefinite". Use a simpler expression for the resulting NaN, as above. Testing showed that working around the compiler bugs is not needed here or in some other places. Here this is because everything ends up as NaN, so the spurious NaNs introduced by gcc have low enough precedence and/or differences that they don't affect the result. % @@ -96,8 +125,10 @@ % */ % if (ix >= 0x7ff00000 && ((hx&0xfffff)|lx) == 0) { % - if ((iy|ly) == 0) % - return (x * x); % - else if (iy >= 0x7ff00000) % - return (x * x) + I * (y - y) / (y - y); % + if ((iy|ly) == 0) { % + creal(f) = x * x; % + cimag(f) = con * y; % + return f; The usual workaround, plus fix the sign for the imaginary part. % + } else if (iy >= 0x7ff00000) % + return x * x + I * (y - y); Use fewer parentheses for the real part and a simpler formula for the imaginary part. % return (x * x) * (cos(y) + I * con * sin(y)); % } % @@ -108,7 +139,11 @@ % */ % if (ix >= 0x7ff00000 && ((hx&0xfffff)|lx) != 0) { % - if ((iy|ly) == 0) % - return x * x; % - return x * x + I * x * x; % + if ((iy|ly) == 0) { % + creal(f) = x - x; % + cimag(f) = copysign(0, x - x); % + return f; The usual workaround, plus fix the sign for the imaginary part. x is NaN, so `con' doesn't apply. % + } % + return (x - x) + (y - y) + I * ((x - x) + (y - y)); Here x and y are both NaNs. Mix both of them into the output. The parentheses are essential, perhaps due to extra precision in intermediate values. On i386's, this results in the largest (with the NaNs considered to be (signed?) 64-bit integers) being propagated. Propagation of NaNs doesn't always preserve them. On i386's the only changes is to convert signaling NaNs to quiet NaNs by setting the quietness bit. We use the formula (y - y) extensively to convert infinities to NaNs and get whatever MD conversions the hardware wants to do to propagate NaNs. % } % + abort(); This should be unreachble, so the previous if statement must always be redundant. % } --- Second, a simple version: % /* % * Hyperbolic cosine of a complex argument. % * % * Most exceptional values are automatically correctly handled by the % * standard formula. See n689.pdf for details. % */ % % #include % #include % % double complex % ccosh1(double complex z) % { % double complex f; % 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. The result must be almost pure % * real or a pure imaginary, except it has type double 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. % * % * In theory, any of sinh(x), cos(y) and sin(y) may be +-0 for a % * nonzero arg, because although the infinitely precise value cannot % * even be rational (Gelfond-Schneider theorem?), the rounded result % * may be 0. For such values, we would have to avoid multiplying % * +-0 by an infinity or a NaN, and cos() and sin() would have to % * be careful to get the sign right. Exhaustive checking shows that % * this never happens for args that are representable as floats. % * % * Note: the expression u + I*v is unusable since gcc introduces % * many overflow, underflow, sign and efficiency bugs by rewriting % * I*v as (0+I)*(v+0*I) and laboriously computing the full complex % * product. In particular, I*Inf is corrupted to NaN+I*Inf, and % * I*-0 is corrupted to 0+I*0. % */ % if (x == 0) { % creal(f) = cos(y); % cimag(f) = isfinite(y) ? x * copysign(1, y) : % copysign(0, y - y); % return (f); % } % if (y == 0) { % creal(f) = cosh(x); % cimag(f) = !isnan(x) ? copysign(1, x) * y : % copysign(0, x - x); % return (f); % } % % creal(f) = cosh(x) * cos(y); % cimag(f) = sinh(x) * sin(y); % return (f); % } I hope this is described in too much detail in its main comment. Third, a simple test program that compares implementations of ccosh() on a reasonably large set of float complex args (including +0, +-Inf and many NaNs): % #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; % vf0 = FUNC(z); % vf1 = __CONCAT(FUNC, 1)(z); % 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 > 2 || eri > 2) { % printf("z = %#8x+I%-#8x %.15g+I%-.15g\n", % xu, yu, x, y); % printf( % __XSTRING(FUNC) "(z) = %#8x+I%-#8x %.15g+I%-.15g\n", % vf0ru, vf0iu, vf0r, vf0i); % printf( % __XSTRING(FUNC) "1(z) = %#8x+I%-#8x %.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 == 0) { % 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" #ifdef nothere % #include "s_ccoshf.c" % #include "s_ccosh1f.c" #endif % % /* 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 This must be built in a subdir of lib/msun/src. I used this mainly with ccosh[1]f() and hacked it to work with ccosh[1](). float complex arg space is too large to check exhaustively, and double complex arg space is much larger and harder to step through and display in the error messages, so I kept using float complex arg space. The coverage with stride 0x80000 should be more than enough for NaNs. Testing with stride 0x80000 showed a maximum error of 2 ulps (for the real and imaginary parts of the result) from using coshf(x) * tanhf(x) instead of sinhf(x) if the non-complex float functions are compiled with -O, but the error is 3-4 ulps if they are compiled with -O0. Bruce