From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 22:58:04 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id 9C384106566B for ; Sun, 12 Aug 2012 22:58:04 +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 EE1AA8FC18 for ; Sun, 12 Aug 2012 22:58:03 +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 q7CMw3Hw075516 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 08:58:03 +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 q7CMvuk2020909 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 08:57:56 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CMvuQ9020908 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 08:57:56 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 08:57:56 +1000 Resent-Message-ID: <20120812225756.GK20453@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 q6I41xco090771 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Wed, 18 Jul 2012 14:01:59 +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 q6I41xkX071239 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Wed, 18 Jul 2012 14:01:59 +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 q6I41gtH021603 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Wed, 18 Jul 2012 14:01:43 +1000 From: Bruce Evans Mail-Followup-To: freebsd-numerics@freebsd.org X-X-Sender: bde@besplex.bde.org To: Peter Jeremy In-Reply-To: <20120718001337.GA87817@server.rulingia.com> Message-ID: <20120718123627.D1575@besplex.bde.org> References: <20120713155805.GC81965@zim.MIT.EDU> <20120714120432.GA70706@server.rulingia.com> <20120717084457.U3890@besplex.bde.org> <5004A5C7.1040405@missouri.edu> <5004DEA9.1050001@missouri.edu> <20120717040118.GA86840@troutmask.apl.washington.edu> <20120717042125.GF66913@server.rulingia.com> <20120717043848.GB87001@troutmask.apl.washington.edu> <20120717225328.GA86902@server.rulingia.com> <20120717232740.GA95026@troutmask.apl.washington.edu> <20120718001337.GA87817@server.rulingia.com> MIME-Version: 1.0 Content-Type: MULTIPART/MIXED; BOUNDARY="0-1667853581-1342584102=:1575" X-Mailman-Approved-At: Sun, 12 Aug 2012 23:55:59 +0000 Cc: Diane Bruce , Steve Kargl , John Baldwin , David Chisnall , Stephen Montgomery-Smith , Bruce Evans , Bruce Evans , David Schultz , 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 22:58:04 -0000 X-Original-Date: Wed, 18 Jul 2012 14:01:42 +1000 (EST) X-List-Received-Date: Sun, 12 Aug 2012 22:58:04 -0000 This message is in MIME format. The first part should be readable text, while the remaining parts are likely unreadable without MIME-aware tools. --0-1667853581-1342584102=:1575 Content-Type: TEXT/PLAIN; charset=X-UNKNOWN; format=flowed Content-Transfer-Encoding: QUOTED-PRINTABLE On Wed, 18 Jul 2012, Peter Jeremy wrote: > On 2012-Jul-17 16:27:40 -0700, Steve Kargl wrote: >> I won't have time to go over the code in detail until >> this weekend, but a quick peek showed some issues. The >> first is style. Although fdlibm has a rather interest >> coding style, new code should use KNF. > > I hope that was only the function declaration lines. I think > the rest is KNF. From=20style(9): % /* % * Multi-line comments look like this. Make them real sentences. F= ill % * them so they look like real paragraphs. % */ (beware than man(1) has mangled the indentation from 0 to 5 spaces. I added my usual quoting for code of "% ". From=20the new file: % /* % ** Calculate complex arc tangent using the identity: % ** catan(z) =3D -i catanh(iz) % */ Multi-line comments don't look like this. Another style point visible in this comment is how to write 'i' and multiplication. Multiplication by juxtaposition (iz) doesn't work near C code with long identifiers which might be named iz. It probably requires all variable names to be 1 letter in a special font for this use. I tried to use "I z" consistently in comments in c_ccosh*.c and to get everyone to follow this convention, but there are already some inconsistencies, and I now wonder if "z I" is better. The pari presentation uses "*" and puts "I" last, and uses spaces for "+" but not for "*" (e.g., "1 + 2*I"). % double complex % catanh(double complex z) % { % =09double zr, zi; % =09int cr, ci; %=20 % =09zr =3D creal(z); % =09cr =3D fpclassify(zr); % =09zi =3D cimag(z); % =09ci =3D fpclassify(zi); The standard classification macros are good for developing things, but they are very slow. All (?) committed complex functions use hard-coded bit test. These are almost as easy to write as the classification macros. Copy them from c_ccosh.c. Copy variable names from c_ccosh.c too. %=20 % =09/* % =09 * catanh(+0 + i0) returns +0 + i0. % =09 * catanh(+0 + iNaN) returns +0 + iNaN. % =09 */ This looks like the description in C99. ccosh.c uses something like: % =09/* % =09 * ctanh(+0 + I 0) =3D +0 + I 0. % =09 * atanh(+0 + I NaN) =3D +0 + I d(NaN). % =09 */ d(NaN) documents that the NaN returned is not necessarily the same as the NaN in the arg. It should be the original NaN with default conversions. The default conversion should be to only quieten signaling NaNs. Most arches have a quiet bit and the conversion sets this. Other conversions are very MD and related to bugs like giving different results depending on the precision used to pass args. d(NaN) is quite different from dNaN. The latter is a fixed default NaN that normally results from any invalid operation on non-NaNs. This is arcane and I probably got it wrong in many cases. My hope was that someday all of these comments could be turned into meta-info that is used to generate test vectors and assertions and maybe man pages. They don't belong in the code. But to generate test vectors and assertions, they need to be very formal and correct. For man pages, I think I prefer to hard-code the documentation but test that it agrees with the meta-info. % =09if (cr =3D=3D FP_NAN) { % =09=09/* % =09=09 * catanh(NaN + iInf) returns =C2=B10 + i=CF=80/2 % =09=09 * the sign of the real part of the result is not % =09=09 * specified by the standard so return +0. % =09=09 */ The UTF is similar to in C99 where it is used for the "+-" amd "infinity" symbols. It messes up n869.txt too (C and POSIX working group translations to text are poor. IIRC, "+-" gets mangled to "+", and "infinity gets mangled to "0"). Why Inf for the arg and not for the result? Here are my current fixes for committed versions of complex functions. They are mostly to fix comments and to remove redundant classifications. This is far from complete in getting all the functions written in a consistent style. They have some fixes for missing or wrong NaN conversions. Most of the bugs for NaN conversions were found by diffing source files and noticing style differences that were actually code bugs. % Index: s_ccoshf.c % =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D % RCS file: /home/ncvs/src/lib/msun/src/s_ccoshf.c,v % retrieving revision 1.2 % diff -u -2 -r1.2 s_ccoshf.c % --- s_ccoshf.c=0921 Oct 2011 06:29:32 -0000=091.2 % +++ s_ccoshf.c=0925 Oct 2011 14:04:11 -0000 % @@ -26,5 +26,5 @@ %=20 % /* % - * Hyperbolic cosine of a complex argument. See s_ccosh.c for details. % + * Hyperbolic cosine of a float complex argument. See s_ccosh.c for det= ails. % */ %=20 % @@ -63,5 +63,5 @@ % =09=09if (ix < 0x42b17218) { % =09=09=09/* x < 88.7: expf(|x|) won't overflow */ % -=09=09=09h =3D expf(fabsf(x)) * 0.5f; % +=09=09=09h =3D expf(fabsf(x)) * 0.5F; % =09=09=09return (cpackf(h * cosf(y), copysignf(h, x) * sinf(y))); % =09=09} else if (ix < 0x4340b1e7) { % @@ -76,20 +76,28 @@ % =09} %=20 % -=09if (ix =3D=3D 0 && iy >=3D 0x7f800000) % +=09if (ix =3D=3D 0)=09=09=09/* && iy >=3D 0x7f800000 */ % =09=09return (cpackf(y - y, copysignf(0, x * (y - y)))); %=20 % -=09if (iy =3D=3D 0 && ix >=3D 0x7f800000) { % -=09=09if ((hx & 0x7fffff) =3D=3D 0) % +=09if (iy =3D=3D 0) {=09=09=09/* && ix >=3D 0x7f800000 */ % +=09=09if (ix =3D=3D 0x7f800000) % =09=09=09return (cpackf(x * x, copysignf(0, x) * y)); % +=09=09/* % +=09=09 * This does a lot of work to get the same sign as the NaN % +=09=09 * sinhf(x =3D NaN) * (y =3D +-0). % +=09=09 * % +=09=09 * In s_ccosh.c, we document that the sign is unspecified, % +=09=09 * but don't document our choice, and d(NaN) is possibly % +=09=09 * not parenthesized correctly. % +=09=09 */ Float versions are supposed to look as much like double versions as possible, except routine and/or comments should be left out of the float version. IIRC, I wrote the large comment here because I originally wanted to limit this set of changes to mainly this file. % =09=09return (cpackf(x * x, copysignf(0, (x + x) * y))); % =09} %=20 % -=09if (ix < 0x7f800000 && iy >=3D 0x7f800000) % +=09if (ix < 0x7f800000)=09=09/* && iy >=3D 0x7f800000 */ % =09=09return (cpackf(y - y, x * (y - y))); %=20 % -=09if (ix >=3D 0x7f800000 && (hx & 0x7fffff) =3D=3D 0) { % +=09if (ix =3D=3D 0x7f800000) { % =09=09if (iy >=3D 0x7f800000) % -=09=09=09return (cpackf(x * x, x * (y - y))); % -=09=09return (cpackf((x * x) * cosf(y), x * sinf(y))); % +=09=09=09return (cpackf(INFINITY, x * (y - y))); % +=09=09return (cpackf(INFINITY * cosf(y), x * sinf(y))); The handling of infinities was wrong and/or unecessarily complicated. In general, we use the expression (x * x) to turn both +Inf and -Inf into +Inf, and also to generate d(NaN) or dNaN if x is invalid, but if we know that the result is always +Inf then we can just use INFINITY. Similarly, we use the expression (x + x) if we want to preserve the sign of +-Inf but generate NaNs. IIRC, there is a fix for a sign error somewhere in these patches, with the bug caused by using the wrong expressi= on. % =09} %=20 % Index: s_csinh.c % =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D % RCS file: /home/ncvs/src/lib/msun/src/s_csinh.c,v % retrieving revision 1.2 % diff -u -2 -r1.2 s_csinh.c % --- s_csinh.c=0921 Oct 2011 06:29:32 -0000=091.2 % +++ s_csinh.c=093 Jan 2012 06:20:18 -0000 % @@ -26,4 +26,11 @@ %=20 % /* % + * XXX TODO: % + * Change x +-I y to x + I (+-y) or vice versa? We currently use the % + * former former for args and the latter for results. % + * s/the invalid floating-point exception/FE_INVALID/g % + */ Oops, it was this and not putting "I" last that I want to change. "I" in the middle goes better with omitting "*". It serves as a delimiter for the real and complex parts. But when the "+-" os attached to "I", it is not arrached to the complex part. % + % +/* % * Hyperbolic sine of a complex argument z =3D x + i y. % * % @@ -64,5 +71,5 @@ % =09=09if ((iy | ly) =3D=3D 0) % =09=09=09return (cpack(sinh(x), y)); % -=09=09if (ix < 0x40360000)=09/* small x: normal case */ % +=09=09if (ix < 0x40360000)=09/* |x| < 22: normal case */ % =09=09=09return (cpack(sinh(x) * cos(y), cosh(x) * sin(y))); %=20 % @@ -84,14 +91,14 @@ %=20 % =09/* % -=09 * sinh(+-0 +- I Inf) =3D sign(d(+-0, dNaN))0 + I dNaN. % -=09 * The sign of 0 in the result is unspecified. Choice =3D normally % -=09 * the same as dNaN. Raise the invalid floating-point exception. % -=09 * % -=09 * sinh(+-0 +- I NaN) =3D sign(d(+-0, NaN))0 + I d(NaN). % -=09 * The sign of 0 in the result is unspecified. Choice =3D normally % -=09 * the same as d(NaN). % +=09 * sinh(+-0 +- I Inf) =3D +-0 + I dNaN. % +=09 * The sign of 0 in the result is unspecified. Choice =3D same sign % +=09 * as the argument. Raise the invalid floating-point exception. % +=09 * % +=09 * sinh(+-0 +- I NaN) =3D +-0 + I d(NaN). % +=09 * The sign of 0 in the result is unspecified. Choice =3D same sign % +=09 * as the argument. % =09 */ % -=09if ((ix | lx) =3D=3D 0 && iy >=3D 0x7ff00000) % -=09=09return (cpack(copysign(0, x * (y - y)), y - y)); % +=09if ((ix | lx) =3D=3D 0)=09=09/* && iy >=3D 0x7ff00000 */ % +=09=09return (cpack(x, y - y)); %=20 % =09/* Since the sign of the 0 is unspecified, don't use a complicated rule to get a particular one. % @@ -100,9 +107,6 @@ % =09 * sinh(NaN +- I 0) =3D d(NaN) + I +-0. % =09 */ % -=09if ((iy | ly) =3D=3D 0 && ix >=3D 0x7ff00000) { % -=09=09if (((hx & 0xfffff) | lx) =3D=3D 0) % -=09=09=09return (cpack(x, y)); % -=09=09return (cpack(x, copysign(0, y))); % -=09} % +=09if ((iy | ly) =3D=3D 0)=09=09/* && ix >=3D 0x7ff00000 */ % +=09=09return (cpack(x + x, y)); %=20 % =09/* This was missing an operation on x, so signaling NaNs were not properly quieted. I probably missed this bug in early testing because I mostly tested on i386, and on i386 signaling NaNs in < long double precision are quieted just by loading them. This was unnecessarily complicated. % @@ -114,21 +118,21 @@ % =09 * nonzero x. Choice =3D don't raise (except for signaling NaNs). % =09 */ % -=09if (ix < 0x7ff00000 && iy >=3D 0x7ff00000) % -=09=09return (cpack(y - y, x * (y - y))); % +=09if (ix < 0x7ff00000)=09=09/* && iy >=3D 0x7ff00000 */ % +=09=09return (cpack(y - y, y - y)); %=20 % =09/* % =09 * sinh(+-Inf + I NaN) =3D +-Inf + I d(NaN). % -=09 * The sign of Inf in the result is unspecified. Choice =3D normally % -=09 * the same as d(NaN). % +=09 * The sign of Inf in the result is unspecified. Choice =3D same sig= n % +=09 * as the argument. % =09 * % =09 * sinh(+-Inf +- I Inf) =3D +Inf + I dNaN. % -=09 * The sign of Inf in the result is unspecified. Choice =3D always += =2E % -=09 * Raise the invalid floating-point exception. % +=09 * The sign of Inf in the result is unspecified. Choice =3D same sig= n % +=09 * as the argument. Raise the invalid floating-point exception. % =09 * % =09 * sinh(+-Inf + I y) =3D +-Inf cos(y) + I Inf sin(y) % =09 */ % -=09if (ix >=3D 0x7ff00000 && ((hx & 0xfffff) | lx) =3D=3D 0) { % +=09if (ix =3D=3D 0x7ff00000 && lx =3D=3D 0) { % =09=09if (iy >=3D 0x7ff00000) % -=09=09=09return (cpack(x * x, x * (y - y))); % +=09=09=09return (cpack(x, y - y)); % =09=09return (cpack(x * cos(y), INFINITY * sin(y))); % =09} % @@ -145,5 +149,5 @@ % =09 * nonzero y. Choice =3D don't raise (except for signaling NaNs). % =09 */ % -=09return (cpack((x * x) * (y - y), (x + x) * (y - y))); % +=09return (cpack((x + x) * (y - y), (x * x) * (y - y))); % } %=20 % @@ -152,5 +156,5 @@ % { %=20 % -=09/* csin(z) =3D -I * csinh(I * z) */ % +=09/* csin(z) =3D -I * csinh(I * z). */ % =09z =3D csinh(cpack(-cimag(z), creal(z))); % =09return (cpack(cimag(z), -creal(z))); % Index: s_csinhf.c % =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D % RCS file: /home/ncvs/src/lib/msun/src/s_csinhf.c,v % retrieving revision 1.2 % diff -u -2 -r1.2 s_csinhf.c % --- s_csinhf.c=0921 Oct 2011 06:29:32 -0000=091.2 % +++ s_csinhf.c=093 Jan 2012 06:20:30 -0000 % @@ -26,5 +26,5 @@ %=20 % /* % - * Hyperbolic sine of a complex argument z. See s_csinh.c for details. % + * Float version of csinh(). See s_csinh.c for details. % */ %=20 % @@ -57,5 +57,5 @@ % =09=09if (iy =3D=3D 0) % =09=09=09return (cpackf(sinhf(x), y)); % -=09=09if (ix < 0x41100000)=09/* small x: normal case */ % +=09=09if (ix < 0x41100000)=09/* |x| < 9: normal case */ % =09=09=09return (cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y))); %=20 % @@ -63,5 +63,5 @@ % =09=09if (ix < 0x42b17218) { % =09=09=09/* x < 88.7: expf(|x|) won't overflow */ % -=09=09=09h =3D expf(fabsf(x)) * 0.5f; % +=09=09=09h =3D expf(fabsf(x)) * 0.5F; % =09=09=09return (cpackf(copysignf(h, x) * cosf(y), h * sinf(y))); % =09=09} else if (ix < 0x4340b1e7) { % @@ -76,23 +76,20 @@ % =09} %=20 % -=09if (ix =3D=3D 0 && iy >=3D 0x7f800000) % -=09=09return (cpackf(copysignf(0, x * (y - y)), y - y)); % +=09if (ix =3D=3D 0)=09=09=09/* && iy >=3D 0x7f800000 */ % +=09=09return (cpackf(x, y - y)); %=20 % -=09if (iy =3D=3D 0 && ix >=3D 0x7f800000) { % -=09=09if ((hx & 0x7fffff) =3D=3D 0) % -=09=09=09return (cpackf(x, y)); % -=09=09return (cpackf(x, copysignf(0, y))); % -=09} % +=09if (iy =3D=3D 0)=09=09=09/* && ix >=3D 0x7f800000 */ % +=09=09return (cpackf(x + x , y)); %=20 % -=09if (ix < 0x7f800000 && iy >=3D 0x7f800000) % -=09=09return (cpackf(y - y, x * (y - y))); % +=09if (ix < 0x7f800000)=09=09/* && iy >=3D 0x7f800000 */ % +=09=09return (cpackf(y - y, y - y)); %=20 % -=09if (ix >=3D 0x7f800000 && (hx & 0x7fffff) =3D=3D 0) { % +=09if (ix =3D=3D 0x7f800000) { % =09=09if (iy >=3D 0x7f800000) % -=09=09=09return (cpackf(x * x, x * (y - y))); % +=09=09=09return (cpackf(x, y - y)); % =09=09return (cpackf(x * cosf(y), INFINITY * sinf(y))); % =09} %=20 % -=09return (cpackf((x * x) * (y - y), (x + x) * (y - y))); % +=09return (cpackf((x + x) * (y - y), (x * x) * (y - y))); The last change (swapping "*"s and "+"s might be only to be logically correct. The sign behaviour for sinh() is quited different than for cosh(), and too much code was copied from c_cosh*.c to create c_sinh*.c. % } %=20 % Index: s_csqrt.c % =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D % RCS file: /home/ncvs/src/lib/msun/src/s_csqrt.c,v % retrieving revision 1.4 % diff -u -2 -r1.4 s_csqrt.c % --- s_csqrt.c=098 Aug 2008 00:15:16 -0000=091.4 % +++ s_csqrt.c=0925 Oct 2011 14:49:27 -0000 % @@ -64,5 +64,5 @@ % =09if (isnan(a)) { % =09=09t =3D (b - b) / (b - b);=09/* raise invalid if b is not a NaN */ % -=09=09return (cpack(a, t));=09/* return NaN + NaN i */ % +=09=09return (cpack(a + a, t));=09/* return NaN + NaN i */ % =09} % =09if (isinf(a)) { All the csqrt() functions were missing NaN quieting. Their comment style is still quite different. % Index: s_csqrtf.c % =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D % RCS file: /home/ncvs/src/lib/msun/src/s_csqrtf.c,v % retrieving revision 1.3 % diff -u -2 -r1.3 s_csqrtf.c % --- s_csqrtf.c=098 Aug 2008 00:15:16 -0000=091.3 % +++ s_csqrtf.c=0925 Oct 2011 14:49:51 -0000 % @@ -55,5 +55,5 @@ % =09if (isnan(a)) { % =09=09t =3D (b - b) / (b - b);=09/* raise invalid if b is not a NaN */ % -=09=09return (cpackf(a, t));=09/* return NaN + NaN i */ % +=09=09return (cpackf(a + a, t));=09/* return NaN + NaN i */ % =09} % =09if (isinf(a)) { % Index: s_csqrtl.c % =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D % RCS file: /home/ncvs/src/lib/msun/src/s_csqrtl.c,v % retrieving revision 1.2 % diff -u -2 -r1.2 s_csqrtl.c % --- s_csqrtl.c=098 Aug 2008 00:15:16 -0000=091.2 % +++ s_csqrtl.c=0925 Oct 2011 14:50:01 -0000 % @@ -64,5 +64,5 @@ % =09if (isnan(a)) { % =09=09t =3D (b - b) / (b - b);=09/* raise invalid if b is not a NaN */ % -=09=09return (cpackl(a, t));=09/* return NaN + NaN i */ % +=09=09return (cpackl(a + a, t));=09/* return NaN + NaN i */ % =09} % =09if (isinf(a)) { % Index: s_ctanh.c % =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D % RCS file: /home/ncvs/src/lib/msun/src/s_ctanh.c,v % retrieving revision 1.2 % diff -u -2 -r1.2 s_ctanh.c % --- s_ctanh.c=0921 Oct 2011 06:30:16 -0000=091.2 % +++ s_ctanh.c=0925 Oct 2011 14:30:18 -0000 % @@ -86,4 +86,15 @@ %=20 % =09/* % +=09 * XXX this is missing the dNaN/d(NaN) notation, which tells us the % +=09 * following: % +=09 * dNaN is a default NaN unrelated to any NaN args % +=09 * d(NaN) is a unary conversion (usually quieting) of the arg `NaN' % +=09 * % +=09 * XXX everything is missing: % +=09 * d(NaN1, NaN2) and d(NaN, y) % +=09 * which should be used for binary conversions. % +=09 * % +=09 * XXX this misspells I as i. % +=09 * % =09 * ctanh(NaN + i 0) =3D NaN + i 0 % =09 * % @@ -103,5 +114,5 @@ % =09if (ix >=3D 0x7ff00000) { % =09=09if ((ix & 0xfffff) | lx)=09/* x is NaN */ % -=09=09=09return (cpack(x, (y =3D=3D 0 ? y : x * y))); % +=09=09=09return (cpack(x + x, y =3D=3D 0 ? y : x * y)); % =09=09SET_HIGH_WORD(x, hx - 0x40000000);=09/* x =3D copysign(1, x) */ % =09=09return (cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y)))); % Index: s_ctanhf.c % =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D= =3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D % RCS file: /home/ncvs/src/lib/msun/src/s_ctanhf.c,v % retrieving revision 1.2 % diff -u -2 -r1.2 s_ctanhf.c % --- s_ctanhf.c=0921 Oct 2011 06:30:16 -0000=091.2 % +++ s_ctanhf.c=0925 Oct 2011 14:30:57 -0000 % @@ -52,5 +52,5 @@ % =09if (ix >=3D 0x7f800000) { % =09=09if (ix & 0x7fffff) % -=09=09=09return (cpackf(x, (y =3D=3D 0 ? y : x * y))); % +=09=09=09return (cpackf(x + x, y =3D=3D 0 ? y : x * y)); % =09=09SET_FLOAT_WORD(x, hx - 0x40000000); % =09=09return (cpackf(x, s_ctanh*.c is also missing NaN quieting. I only got my tests in good enough shape to find the NaN quieting bugs for complex functions a little after the above were committed. The tests are still painful to configure for complex and long double functions. Recently I started using ones written in pari for long double functions. They run thousands of times slower than the C tests but work to higher precisions. One for logl looks like this: % \p 100 % \r ../lib/ftoa.gp % \r ../lib/ftoe.gp % \r ../lib/ftof.gp % \r ../lib/roundn.gp %=20 % show(x, y, z) =3D % { % =09local(libm_v, pari_v, relerr, rlibm_v, rpari_v, rx, ry, rz, ulps); %=20 % =09rx =3D roundn(x, prec); % =09ry =3D roundn(y, prec); % =09rz =3D roundn(z, prec); % =09libm_v =3D ry + rz; % =09pari_v =3D log(rx); % =09rlibm_v =3D roundn(libm_v, prec); % =09rpari_v =3D roundn(pari_v, prec); % =09if (rlibm_v !=3D rpari_v || % =09 roundn(libm_v, prec + 7) !=3D roundn(pari_v, prec + 7), % =09=09print("rx =3D " ftoe(rx, 36)); % =09=09print("rx =3D " ftoa(rx, prec)); % =09=09print("ry =3D " ftoa(ry, prec)); % =09=09print("rz =3D " ftoa(rz, prec)); % =09=09print("libm log(x) =3D " ftoa(libm_v, prec + 16)); % =09=09print("rlibm log(x) =3D " ftoa(rlibm_v, prec)); % =09=09print("pari log(x) =3D " ftoa(pari_v, prec + 16)); % =09=09print("rpari log(x) =3D " ftoa(rpari_v, prec)); % =09=09relerr =3D log(abs(libm_v / pari_v - 1)) / log(2); % =09=09print("relerr =3D 2**" ftof(relerr, 5)); % =09=09ulps =3D abs(rlibm_v - pari_v) / % =09=09 2^(ceil(log(abs(pari_v)) / log(2)) - prec); % =09=09print("rnderr =3D " ftof(ulps, 5) " ulps"); % =09); % } %=20 % prec =3D 113; % \r gin This takes an input file named gin that looks like this: % x =3D 0.998031648000081380208333333333333327; y =3D -0.001968431290575624= 91416520767945158578; z =3D -1.86045996676022698958009120338776959e-06; sho= w(x, y, z); % x =3D 0.998031618197758992513020833333333327; y =3D -0.001968461209313334= 43641252213257658578; z =3D -1.86040232903924733902992614348864124e-06; sho= w(x, y, z); % [... thousands or millions of lines. Takes to long to do the billions of cases routinely done in lower precisions by C tests] It produces output like this: % rx =3D 9.98039664824803670247395833333333327e-1 % rx =3D 0x1feff0e1111111111111111111111.0p-113 % ry =3D -0x100f38b08ee99e27e6493479dddef.0p-121 % rz =3D -0x1f7958665a30c56debcb8d67bfc00.0p-132 % libm log(x) =3D -0x101327db9bb4e4409406adeb8ad6e8000.0p-137 % rlibm log(x) =3D -0x101327db9bb4e4409406adeb8ad6e.0p-121 % pari log(x) =3D -0x101327db9bb4e4409406adeb8ad6e802e.0p-137 % rpari log(x) =3D -0x101327db9bb4e4409406adeb8ad6f.0p-121 % relerr =3D 2**-122.49 % rnderr =3D 0.50070 ulps Bruce --0-1667853581-1342584102=:1575--