Use of C99 extra long double math functions after r236148

[Peter Jeremy]

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 = 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)) == 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*π/2     0+I*π/2     0+I*π/2     0+I*π/2      0+I*π/2
 * NaN    0+I*NaN     0+I*NaN   NaN+I*NaN   NaN+I*NaN    NaN+I*NaN
 *  0     0+I*π/2     0+I*0     Inf+I*0      atanh(x)    NaN+I*NaN
 *  y     0+I*π/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 = creal(z);
	y = cimag(z);

	EXTRACT_WORDS(hx, lx, x);
	EXTRACT_WORDS(hy, ly, y);

	ix = 0x7fffffff & hx;
	iy = 0x7fffffff & hy;

	/* Handle pure real & imaginary cases */
	if ((iy | ly) == 0) {	/* Imaginary part 0 */
		/* z is real - return atanh(x) */
		return (cpack(__ieee754_atanh(x), y));
	}

	if ((ix | lx) == 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) == 0) { /* y is Inf */
		/*
		 * catanh(NaN + I*Inf)    = 0 + I*π/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)    = 0 + I*π/2.
		 * catanh(finite + I*Inf) = 0 + I*π/2.
		 */
		return (cpack(copysign(0.0, x), copysign(M_PI_2, y)));
	}

	if (((ix - 0x7ff00000) | lx) == 0) { /* x is Inf */
		/* catanh(Inf + I*finite) = 0 + I*π/2 */
		if (iy < 0x7ff00000)	/* finite */
			return (cpack(copysign(0.0, x), copysign(M_PI_2, y)));

		/* catanh(Inf + I*NaN) = +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) = d(NaN) + I*dNaN
	 * catanh(NaN + I*NaN)    = d(NaN) + I*d(NaN)
	 * catanh(finite + I*NaN) = dNaN + I*d(NaN)
	 *
	 * Raising "invalid" exception is optional.  Choice = don't
	 * raise, except for signalling NaNs.
	 */
	return(cpack(x + y, y + x));
}

/*
 * Arc-tangent of a complex argument z = x + I*y.
 *
 * catan(z) = -I * catanh(I * z)
 *
 */
double complex
catan(double complex z)
{
	double complex r;

	/* Manually multiply by I to avoid compiler deficiencies. */
	r = catanh(cpack(-cimag(z), creal(z)));
	return (cpack(cimag(r), -creal(r)));
}

-- 
Peter Jeremy
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 196 bytes
Desc: not available
Url : http://lists.freebsd.org/pipermail/freebsd-numerics/attachments/20120812/d02b9c7f/attachment.pgp