From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:09:59 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 27BDE1065674 for ; Sun, 12 Aug 2012 23:09:59 +0000 (UTC) (envelope-from peter@rulingia.com) Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by mx1.freebsd.org (Postfix) with ESMTP id A1EB98FC16 for ; Sun, 12 Aug 2012 23:09:58 +0000 (UTC) Received: from server.rulingia.com (c220-239-249-137.belrs5.nsw.optusnet.com.au [220.239.249.137]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CN9waE075800 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:09:58 +1000 (EST) (envelope-from peter@rulingia.com) X-Bogosity: Ham, spamicity=0.000000 Received: from server.rulingia.com (localhost.rulingia.com [127.0.0.1]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q7CN9qYR021565 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:09:52 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CN9qaS021564 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:09:52 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:09:52 +1000 Resent-Message-ID: <20120812230952.GQ20453@server.rulingia.com> Resent-To: freebsd-numerics@freebsd.org Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q6K7xrkM027009 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Fri, 20 Jul 2012 17:59:53 +1000 (EST) (envelope-from brde@optusnet.com.au) Received: from mail05.syd.optusnet.com.au (mail05.syd.optusnet.com.au [211.29.132.186]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q6K7xrK5085896 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Fri, 20 Jul 2012 17:59:53 +1000 (EST) (envelope-from brde@optusnet.com.au) Received: from c122-106-171-246.carlnfd1.nsw.optusnet.com.au (c122-106-171-246.carlnfd1.nsw.optusnet.com.au [122.106.171.246]) by mail05.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id q6K7xboh024498 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Fri, 20 Jul 2012 17:59:38 +1000 From: Bruce Evans Mail-Followup-To: freebsd-numerics@freebsd.org X-X-Sender: bde@besplex.bde.org To: Stephen Montgomery-Smith In-Reply-To: <50085441.4090305@missouri.edu> Message-ID: <20120720162953.N2162@besplex.bde.org> References: <20120714120432.GA70706@server.rulingia.com> <20120717084457.U3890@besplex.bde.org> <5004A5C7.1040405@missouri.edu> <5004DEA9.1050001@missouri.edu> <20120717200931.U6624@besplex.bde.org> <5006D13D.2080702@missouri.edu> <20120718205625.GA409@troutmask.apl.washington.edu> <500725F2.7060603@missouri.edu> <20120719025345.GA1376@troutmask.apl.washington.edu> <50077987.1080307@missouri.edu> <20120719032706.GA1558@troutmask.apl.washington.edu> <5007826D.7060806@missouri.edu> <5007AD41.9070000@missouri.edu> <20120719205347.T2601@besplex.bde.org> <50084322.7020401@missouri.edu> <20120720035001.W4053@besplex.bde.org> <50085441.4090305@missouri.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed X-Mailman-Approved-At: Sun, 12 Aug 2012 23:55:59 +0000 Cc: Diane Bruce , Bruce Evans , John Baldwin , David Chisnall , Bruce Evans , Steve Kargl , David Schultz , Peter Jeremy , Warner Losh Subject: Re: Use of C99 extra long double math functions after r236148 X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Date: Sun, 12 Aug 2012 23:09:59 -0000 X-Original-Date: Fri, 20 Jul 2012 17:59:37 +1000 (EST) X-List-Received-Date: Sun, 12 Aug 2012 23:09:59 -0000 On Thu, 19 Jul 2012, Stephen Montgomery-Smith wrote: First a bit from a previous reply: bde: >> WIll look more closely later. I see that you already use log1p() and much >> more. Apple clog() uses not so much more, mainly by depending on extra >> precision in hardware. The above also avoids overflow and use of hypot() >> for all finite x and and y, but is probably too simple. Now looked more closely. I put the real part of it it in my test framework as loghypot() and tested this against alternative versions like ordinary hypot(). stephen: > The Apple solution has a problem. The two invocations of log might > produce results that are nearly identical, but with opposite signs. > Think about x = y = 1/sqrt(2). Simple log(hypot()) has huge errors, but the Apple version seemed to be generally accurate and was much more accurate near |z| = 1. This pointed to some bugs in your version. These show up before any cancelation bugs in the Apple version. > I think their solution merely avoids the overflow/underflow problem, and was > not meant to address the problem I worked on. However, their solution will > fail if y=1e100 and y=1e-100. > > This caused me to realize that my solution failed to account for underflow, > so here is my next iteration. Tested this version. First I reformatted it all (the comment indentation was furthest from KNF). indent(1) did a good job and I only reformatted a couple of long lines manually: % #include % #include % #include % #include % % #include "math_private.h" % % /* % * gcc doesn't implement complex multiplication or division correctly, so we % * need to handle infinities specially. We turn on this pragma to notify % * conforming c99 compilers that the fast-but-incorrect code that gcc % * generates is acceptable, since the special cases have already been % * handled. % */ % #pragma STDC CX_LIMITED_RANGE ON % % double complex % clog(double complex z) % { % double x, y, h, t1, t2, t3; % double x0, y0, x1, y1; % % x = creal(z); % y = cimag(z); % % #define NANMIX_APPLE_CLOG_COMPAT 1 % /* Handle special cases when x or y is NAN. */ % if (isnan(x)) { % if (isinf(y)) % return (cpack(INFINITY, NAN)); % else { % #if NANMIX_HYPOTF_COMPAT % y = fabs(y + 0); % t1 = (y - y) / (y - y); /* Raise invalid flag if y is % * not NaN */ % t1 = fabs(x + 0) + t1; /* Mix NaN(s). */ % #elif NANMIX_APPLE_CLOG_COMPAT % /* No actual mixing. */ % return (cpack(x, copysign(x, y))); % #else % t1 = (y - y) / (y - y); /* Raise invalid flag if y is % * not NaN */ % #endif % return (cpack(t1, t1)); % } % } else if (isnan(y)) { % if (isinf(x)) % return (cpack(INFINITY, NAN)); % else { % #ifdef NANMIX_HYPOTF_COMPAT % x = fabs(x + 0); % t1 = (x - x) / (x - x); /* Raise invalid flag if x is % * not NaN */ % t1 = t1 + fabs(y + 0); /* Mix NaN(s). */ % #elif NANMIX_APPLE_CLOG_COMPAT % /* No actual mixing. */ % return (cpack(y, y)); % #else % t1 = (x - x) / (x - x); /* Raise invalid flag if x is % * not NaN */ % #endif % return (cpack(t1, t1)); % } Here I added sqillions of ifdefs so that the results come out the same for NaNs as with log(hypot()) and with apple clog(). This is not very important, but my tests report all differences and it is hard to see important differences when there are many for NaNs. In general, the result of an operation with 1 NaN are should be the same NaN quieted, and the result of an operation with 2 NaN args should combine both args in a similar way to a hardware operation. For hypot(x, y), the result should be the same as from the naive formula sqrt(x*x+y*y) executed by the hardware. The real part of clog() should be the result of log() on this (according to the simpler rules for 1 NaN). Unfortunately, fdlibm hypot() historically first replaces x and y by |x| and |y|. fabs() clears the sign bit even for NaNs, and the result is quite different from that given by the naive formula. The above adds some fabs()'s and other tweaks for consistency with hypot(). I don't try hard to make the imaginary part consistent with atan2(). It would be simpler to just call hypot() for NaN cases. Apple clog() mostly doesn't mix pairs of NaNs. It mostly doesn't even try to apply the default conversions to a single NaN. The above adds another set of tweaks for consistency with it. % } else if (isfinite(x) && isfinite(y) && % (fabs(x) > 1e308 || fabs(y) > 1e308)) % /* % * To avoid unnecessary overflow, if x or y are very large, % * divide x and y by M_E, and then add 1 to the logarithm. % * This depends on M_E being larger than sqrt(2). % */ % return (cpack(log(hypot(x / M_E, y / M_E) + 1), atan2(y, x))); % else if (fabs(x) < 1e-50 || fabs(y) < 1e-50 || % fabs(x) > 1e50 || fabs(y) > 1e50) % /* % * Because atan2 and hypot conform to C99, this also covers % * all the edge cases when x or y are 0 or infinite. % */ % return (cpack(log(hypot(x, y)), atan2(y, x))); % else { % /* We don't need to worry about overflow in x*x+y*y. */ % h = x * x + y * y; % if (h < 0.1 || h > 10) % return (cpack(log(h) / 2, atan2(y, x))); % /* Take extra care if h is moderately close to 1 */ % else { % /* % * x0 and y0 are good approximations to x and y, but % * have their bits trimmed so that double precision % * floating point is capable of calculating x0*x0 + % * y0*y0 - 1 exactly. % */ % x0 = ldexp(floor(ldexp(x, 24)), -24); % x1 = x - x0; % y0 = ldexp(floor(ldexp(y, 24)), -24); % y1 = y - y0; This fails for example when x is near 1 and y is tiny. Say x == 1 and y == +-1e-19. The result for the real part should be about y**2/2 (almost exactly for y this tiny I think), but the above makes the following messes, since the ldexp()s don't actually trim to 24 bits: - for tiny y > 0, y times 2**24 is < 1 and floor() gives 0; x0 is 0 and x1 is x. The error has been moved to a worse place. - for tiny y < 0, y times 2**24 is > -1 and floor() gives -1; x0 is -2**-24, which is much larger in magnitude than y, and x1 is about as large to compensate. The error has been increased as well as moved. Since y is tiny, the final result was just 0 in all cases that I looked at. Apple clog() handles this case by using the y**2/2 approximation. The result of 0 instead of y**2/2 is about 9 GULPS of error, even when eveything is reduced to float precision (multiply this by 2**29 for the error in in double precision. I don't remember the SI prefixes that high). I tried simple fixes, but they made the total number of errors larger. ldexp() is slow, and the following are standard ways of trimming low bits much mure efficiently and correctly: - fdlibm mostly uses GET/EXTRACT to extract the bits, then masks low bits, then uses SET/INSERT to put the changed bits in an FP variable. (The naming differences are historical bugs.) - for double precision, most hi+lo decompositions split the double almost in half, with 26 hi bits and 27 lo bits. But often you only need 24 hi bits. Then it is usually faster to assign the double to a float. to get 24 hi bits. This requires that the double be small enough to be almost representable as a float: x0 = (float)x; x1 = x - x0; Here x0 can be float or the same type as x, or even different. On i386, there are compiler bugs which may require using STRICT_ASSIGN() here (at least if x is an expression). The compiler bugs are best avoided by working throughout with variables of type double_t (for double functions). - a more general method is to add and subtract a magic number like 0x1p40 or 0x1.8p40: x0 = x + 0x1p40 - 0x1p40; x1 = x - x0; This is usually faster, but more fragile. The magic 40 trims to 24 (+-1?) bits, provided the evaluation precision is 64 bits and there are no compiler bugs. On i386, the evaluation precision is variable but usually 53, and libm uses code like the above with magic for 53 in a couple of places. In my attempted simple fix, I used the float cast method. This failed because y = +-1e-19 is too tiny for anything good to happen. The actual y was already exactly representable in 24 bits, so y0 was y and y1 was 0... % /* Notice that mathematically, h = t1*(1+t3). */ % t1 = x0 * x0 + y0 * y0; ... t1 = 1 + 1e-38 rounds to 1, and there is no extra precision. I tried logl(hypotl()). The extra precision provided by long doubles was very far from making much difference. hypot()s methods don't seem to apply, since for tiny y hypot(1, y) ~= 1 + y**2/2 ~= 1, so the case of tiny y is especially easy for hypot(). % t2 = 2 * x0 * x1 + x1 * x1 + 2 * y0 * y1 + y1 * y1; % t3 = t2 / t1; % return (cpack((log1p(t1 - 1) + log1p(t3)) / 2, % atan2(y, x))); Isn't log1p(t1 -1) no different from log(t1)? % } % } % } I was going to leave all the fixes to you, but now want to try the approximation for tiny y :-). Bruce