Date: Sun, 12 Aug 2012 23:00:41 -0000 From: Peter Jeremy <peter@rulingia.com> To: Bruce Evans <brde@optusnet.com.au> Cc: Diane Bruce <db@db.net>, John Baldwin <jhb@freebsd.org>, David Chisnall <theraven@freebsd.org>, Stephen Montgomery-Smith <stephen@missouri.edu>, Bruce Evans <bde@freebsd.org>, Steve Kargl <sgk@troutmask.apl.washington.edu>, David Schultz <das@freebsd.org>, Warner Losh <imp@bsdimp.com> Subject: Re: Use of C99 extra long double math functions after r236148 Message-ID: <20120722121219.GC73662@server.rulingia.com> Resent-Message-ID: <20120812230033.GT20453@server.rulingia.com> In-Reply-To: <20120718123627.D1575@besplex.bde.org> 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>
next in thread | previous in thread | raw e-mail | index | archive | help
--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 <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <complex.h>
#include <math.h>
#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--
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120722121219.GC73662>