From owner-freebsd-numerics@FreeBSD.ORG Sun Aug 12 23:00:40 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 DBC6E1065677 for ; Sun, 12 Aug 2012 23:00:40 +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 45EAD8FC12 for ; Sun, 12 Aug 2012 23:00:40 +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 q7CN0evV075571 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 13 Aug 2012 09:00:40 +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 q7CN0XdK021040 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 13 Aug 2012 09:00:33 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q7CN0XNd021039 for freebsd-numerics@freebsd.org; Mon, 13 Aug 2012 09:00:33 +1000 (EST) (envelope-from peter) Resent-From: Peter Jeremy Resent-Date: Mon, 13 Aug 2012 09:00:33 +1000 Resent-Message-ID: <20120812230033.GT20453@server.rulingia.com> Resent-To: freebsd-numerics@freebsd.org From: Peter Jeremy Mail-Followup-To: freebsd-numerics@freebsd.org To: Bruce Evans Message-ID: <20120722121219.GC73662@server.rulingia.com> References: <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> <20120718123627.D1575@besplex.bde.org> MIME-Version: 1.0 Content-Type: multipart/signed; micalg=pgp-sha1; protocol="application/pgp-signature"; boundary="E13BgyNx05feLLmH" Content-Disposition: inline In-Reply-To: <20120718123627.D1575@besplex.bde.org> X-PGP-Key: http://www.rulingia.com/keys/peter.pgp User-Agent: Mutt/1.5.21 (2010-09-15) Cc: Diane Bruce , John Baldwin , David Chisnall , Stephen Montgomery-Smith , Bruce Evans , Steve Kargl , 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 23:00:41 -0000 X-Original-Date: Sun, 22 Jul 2012 22:12:19 +1000 X-List-Received-Date: Sun, 12 Aug 2012 23:00:41 -0000 --E13BgyNx05feLLmH Content-Type: text/plain; charset=utf-8 Content-Disposition: inline Content-Transfer-Encoding: quoted-printable OK, here's my next try at the exception handling for catanh(). As before, the real code needs expansion to prevent precision loss. A have simplified the default (NaN + I*NaN) return from catanh() to the minimun to ensure that both real & imaginary parts return as NaN. I've been doing some experiments on mixing NaNs using x87, SSE, SPARC64 and ARM (last on Linux) and have come to the conclusion that there is no standard behaviour: Given x & y as NaNs, (x+y) can return either x or y, possibly with the sign bit from the other operand. depending on the FPU. Inline, as rquested by Steve: -------------------- /*Copyright...*/ #include __FBSDID("$FreeBSD$"); #include #include #include "math_private.h" /* * Hyperbolic arc-tangent of a complex argument z =3D x + I*y. * * Exceptional values are noted in the comments within the source code. * These values and the associated return value were taken from WG14/N1256. * It (and the code comments) generally only refers to behaviour in the * quadrant where both input signs are positive. This is extended to the * remaining quadrants by noting: * a) catanh(conj(z)) =3D=3D conj(catanh(z)) * b) catanh() is odd * therefore the result in each quadrant is the same with the signs of * each part copied from the input to the output. * * Special cases for catanh() based on WG14/N1256: * * y\x Inf 0 1 x NaN * Inf 0+I*=CF=80/2 0+I*=CF=80/2 0+I*=CF=80/2 0+I*=CF=80/2 = 0+I*=CF=80/2 * NaN 0+I*NaN 0+I*NaN NaN+I*NaN NaN+I*NaN NaN+I*NaN * 0 0+I*=CF=80/2 0+I*0 Inf+I*0 atanh(x) NaN+I*NaN * y 0+I*=CF=80/2 I*atan(y) [1] [1] NaN+I*NaN * * [1] clog((1+z)/(1-z))/2 or equivalent. */ double complex catanh(double complex z) { double x, y; /* Real & imaginary parts of argument */ int32_t hx, hy; /* MSW of binary real & imaginary parts */ int32_t ix, iy; /* hx & hy without sign bits */ int32_t lx, ly; /* LSW of binary real & imaginary parts */ x =3D creal(z); y =3D cimag(z); EXTRACT_WORDS(hx, lx, x); EXTRACT_WORDS(hy, ly, y); ix =3D 0x7fffffff & hx; iy =3D 0x7fffffff & hy; /* Handle pure real & imaginary cases */ if ((iy | ly) =3D=3D 0) { /* Imaginary part 0 */ /* z is real - return atanh(x) */ return (cpack(__ieee754_atanh(x), y)); } if ((ix | lx) =3D=3D 0) { /* Real part 0 */ /* z is imaginary - return I*atan(y) */ return (cpack(x, atan(y))); } /* Handle the mostly-non-exceptional cases where x and y are finite. */ if (ix < 0x7ff00000 && iy < 0x7ff00000) { /* Following is possible algorithm, not final implementation */ return ((clog(1.0 + z) - clog(1.0 - z)) / 2.0); } /* x and/or y are not finite */ if (((iy - 0x7ff00000) | ly) =3D=3D 0) { /* y is Inf */ /* * catanh(NaN + I*Inf) =3D 0 + I*=CF=80/2. * The sign of the real part of the result is not * specified by the standard so always return same as x. * catanh(Inf + I*Inf) =3D 0 + I*=CF=80/2. * catanh(finite + I*Inf) =3D 0 + I*=CF=80/2. */ return (cpack(copysign(0.0, x), copysign(M_PI_2, y))); } if (((ix - 0x7ff00000) | lx) =3D=3D 0) { /* x is Inf */ /* catanh(Inf + I*finite) =3D 0 + I*=CF=80/2 */ if (iy < 0x7ff00000) /* finite */ return (cpack(copysign(0.0, x), copysign(M_PI_2, y))); /* catanh(Inf + I*NaN) =3D +0 + I*NaN */ return (cpack(copysign(0.0, x), y+y)); } /* * At this point x and/or y are NaN and these all return NaN + I*NaN * * catanh(NaN + I*finite) =3D d(NaN) + I*dNaN * catanh(NaN + I*NaN) =3D d(NaN) + I*d(NaN) * catanh(finite + I*NaN) =3D dNaN + I*d(NaN) * * Raising "invalid" exception is optional. Choice =3D don't * raise, except for signalling NaNs. */ return(cpack(x + y, y + x)); } /* * Arc-tangent of a complex argument z =3D x + I*y. * * catan(z) =3D -I * catanh(I * z) * */ double complex catan(double complex z) { double complex r; /* Manually multiply by I to avoid compiler deficiencies. */ r =3D catanh(cpack(-cimag(z), creal(z))); return (cpack(cimag(r), -creal(r))); } --=20 Peter Jeremy --E13BgyNx05feLLmH Content-Type: application/pgp-signature -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.19 (FreeBSD) iEYEARECAAYFAlAL7iMACgkQ/opHv/APuIeydwCgr6tuTbwOgBsWgBvK+3bEa/AZ ovsAnjqXmugJGD+ByyBsIHPtRG1zL3pG =kKBH -----END PGP SIGNATURE----- --E13BgyNx05feLLmH--