Date: Sun, 12 Aug 2012 22:58:04 -0000 From: Bruce Evans <brde@optusnet.com.au> To: Peter Jeremy <peter@rulingia.com> Cc: Diane Bruce <db@db.net>, Steve Kargl <sgk@troutmask.apl.washington.edu>, John Baldwin <jhb@freebsd.org>, David Chisnall <theraven@freebsd.org>, Stephen Montgomery-Smith <stephen@missouri.edu>, Bruce Evans <bde@freebsd.org>, Bruce Evans <brde@optusnet.com.au>, David Schultz <das@freebsd.org>, Warner Losh <imp@bsdimp.com> Subject: Re: Use of C99 extra long double math functions after r236148 Message-ID: <20120718123627.D1575@besplex.bde.org> Resent-Message-ID: <20120812225756.GK20453@server.rulingia.com> In-Reply-To: <20120718001337.GA87817@server.rulingia.com> 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>
next in thread | previous in thread | raw e-mail | index | archive | help
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 <sgk@troutmask.apl.washington.=
edu> 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--
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120718123627.D1575>