I'm happy to finally commit stephen@'s implementations of cacos,

cacosh, casin, casinh, catan, and catanh. Thanks to stephen@ and bde@
for working on these.

Submitted by:	stephen@
Reviewed by:	bde
This commit is contained in:
David Schultz 2013-05-30 04:49:26 +00:00
parent 6bba248bee
commit e4afa19c33
9 changed files with 1198 additions and 15 deletions

View File

@ -63,9 +63,21 @@ __BEGIN_DECLS
double cabs(double complex);
float cabsf(float complex);
long double cabsl(long double complex);
double complex cacos(double complex);
float complex cacosf(float complex);
double complex cacosh(double complex);
float complex cacoshf(float complex);
double carg(double complex);
float cargf(float complex);
long double cargl(long double complex);
double complex casin(double complex);
float complex casinf(float complex);
double complex casinh(double complex);
float complex casinhf(float complex);
double complex catan(double complex);
float complex catanf(float complex);
double complex catanh(double complex);
float complex catanhf(float complex);
double complex ccos(double complex);
float complex ccosf(float complex);
double complex ccosh(double complex);

View File

@ -105,7 +105,8 @@ COMMON_SRCS+= e_acosl.c e_asinl.c e_atan2l.c e_fmodl.c \
.endif
# C99 complex functions
COMMON_SRCS+= s_ccosh.c s_ccoshf.c s_cexp.c s_cexpf.c \
COMMON_SRCS+= catrig.c catrigf.c \
s_ccosh.c s_ccoshf.c s_cexp.c s_cexpf.c \
s_cimag.c s_cimagf.c s_cimagl.c \
s_conj.c s_conjf.c s_conjl.c \
s_cproj.c s_cprojf.c s_creal.c s_crealf.c s_creall.c \
@ -126,7 +127,7 @@ SRCS= ${COMMON_SRCS} ${ARCH_SRCS}
INCS+= fenv.h math.h
MAN= acos.3 acosh.3 asin.3 asinh.3 atan.3 atan2.3 atanh.3 \
ceil.3 ccos.3 ccosh.3 cexp.3 \
ceil.3 cacos.3 ccos.3 ccosh.3 cexp.3 \
cimag.3 copysign.3 cos.3 cosh.3 csqrt.3 erf.3 exp.3 fabs.3 fdim.3 \
feclearexcept.3 feenableexcept.3 fegetenv.3 \
fegetround.3 fenv.3 floor.3 \
@ -144,6 +145,9 @@ MLINKS+=atan.3 atanf.3 atan.3 atanl.3
MLINKS+=atanh.3 atanhf.3
MLINKS+=atan2.3 atan2f.3 atan2.3 atan2l.3 \
atan2.3 carg.3 atan2.3 cargf.3 atan2.3 cargl.3
MLINKS+=cacos.3 cacosf.3 cacos.3 cacosh.3 cacos.3 cacoshf.3 \
cacos.3 casin.3 cacos.3 casinf.3 cacos.3 casinh.3 cacos.3 casinhf.3 \
cacos.3 catan.3 cacos.3 catanf.3 cacos.3 catanh.3 cacos.3 catanhf.3
MLINKS+=ccos.3 ccosf.3 ccos.3 csin.3 ccos.3 csinf.3 ccos.3 ctan.3 ccos.3 ctanf.3
MLINKS+=ccosh.3 ccoshf.3 ccosh.3 csinh.3 ccosh.3 csinhf.3 \
ccosh.3 ctanh.3 ccosh.3 ctanhf.3

View File

@ -237,6 +237,18 @@ FBSD_1.3 {
fegetround;
fesetround;
fesetenv;
cacos;
cacosf;
cacosh;
cacoshf;
casin;
casinf;
casinh;
casinhf;
catan;
catanf;
catanh;
catanhf;
csin;
csinf;
csinh;

128
lib/msun/man/cacos.3 Normal file
View File

@ -0,0 +1,128 @@
.\" Copyright (c) 2013 David Schultz <das@FreeBSD.org>
.\" All rights reserved.
.\"
.\" Redistribution and use in source and binary forms, with or without
.\" modification, are permitted provided that the following conditions
.\" are met:
.\" 1. Redistributions of source code must retain the above copyright
.\" notice, this list of conditions and the following disclaimer.
.\" 2. Redistributions in binary form must reproduce the above copyright
.\" notice, this list of conditions and the following disclaimer in the
.\" documentation and/or other materials provided with the distribution.
.\"
.\" THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
.\" ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
.\" IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
.\" ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
.\" FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
.\" DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
.\" OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
.\" HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
.\" LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
.\" OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
.\" SUCH DAMAGE.
.\"
.\" $FreeBSD$
.\"
.Dd May 27, 2013
.Dt CACOS 3
.Os
.Sh NAME
.Nm cacos ,
.Nm cacosf ,
.Nm cacosh ,
.Nm cacoshf ,
.Nm casin ,
.Nm casinf
.Nm casinh ,
.Nm casinhf
.Nm catan ,
.Nm catanf
.Nm catanh ,
.Nm catanhf
.Nd complex arc trigonometric and hyperbolic functions
.Sh LIBRARY
.Lb libm
.Sh SYNOPSIS
.In complex.h
.Ft double complex
.Fn cacos "double complex z"
.Ft float complex
.Fn cacosf "float complex z"
.Ft double complex
.Fn cacosh "double complex z"
.Ft float complex
.Fn cacoshf "float complex z"
.Ft double complex
.Fn casin "double complex z"
.Ft float complex
.Fn casinf "float complex z"
.Ft double complex
.Fn casinh "double complex z"
.Ft float complex
.Fn casinhf "float complex z"
.Ft double complex
.Fn catan "double complex z"
.Ft float complex
.Fn catanf "float complex z"
.Ft double complex
.Fn catanh "double complex z"
.Ft float complex
.Fn catanhf "float complex z"
.Sh DESCRIPTION
The
.Fn cacos ,
.Fn casin ,
and
.Fn catan
functions compute the principal value of the inverse cosine, sine,
and tangent of the complex number
.Fa z ,
respectively.
The
.Fn cacosh ,
.Fn casinh ,
and
.Fn catanh
functions compute the principal value of the inverse hyperbolic
cosine, sine, and tangent.
The
.Fn cacosf ,
.Fn casinf ,
.Fn catanf
.Fn cacoshf ,
.Fn casinhf ,
and
.Fn catanhf
functions perform the same operations in
.Fa float
precision.
.Pp
.ie '\*[.T]'utf8'
. ds Un \[cu]
.el
. ds Un U
.
There is no universal convention for defining the principal values of
these functions. The following table gives the branch cuts, and the
corresponding ranges for the return values, adopted by the C language.
.Bl -column ".Sy Function" ".Sy (-\*(If*I, -I) \*(Un (I, \*(If*I)" ".Sy [-\*(Pi/2*I, \*(Pi/2*I]"
.It Sy Function Ta Sy Branch Cut(s) Ta Sy Range
.It cacos Ta (-\*(If, -1) \*(Un (1, \*(If) Ta [0, \*(Pi]
.It casin Ta (-\*(If, -1) \*(Un (1, \*(If) Ta [-\*(Pi/2, \*(Pi/2]
.It catan Ta (-\*(If*I, -i) \*(Un (I, \*(If*I) Ta [-\*(Pi/2, \*(Pi/2]
.It cacosh Ta (-\*(If, 1) Ta [-\*(Pi*I, \*(Pi*I]
.It casinh Ta (-\*(If*I, -i) \*(Un (I, \*(If*I) Ta [-\*(Pi/2*I, \*(Pi/2*I]
.It catanh Ta (-\*(If, -1) \*(Un (1, \*(If) Ta [-\*(Pi/2*I, \*(Pi/2*I]
.El
.Sh SEE ALSO
.Xr ccos 3 ,
.Xr ccosh 3 ,
.Xr complex 3 ,
.Xr cos 3 ,
.Xr math 3 ,
.Xr sin 3 ,
.Xr tan 3
.Sh STANDARDS
These functions conform to
.St -isoC-99 .

View File

@ -69,6 +69,7 @@ functions perform the same operations in
.Fa float
precision.
.Sh SEE ALSO
.Xr cacos 3 ,
.Xr ccosh 3 ,
.Xr complex 3 ,
.Xr cos 3 ,

View File

@ -69,6 +69,7 @@ functions perform the same operations in
.Fa float
precision.
.Sh SEE ALSO
.Xr cacosh 3 ,
.Xr ccos 3 ,
.Xr complex 3 ,
.Xr cosh 3 ,

View File

@ -89,6 +89,12 @@ creal compute the real part
.\" Section 7.3.5-6 of ISO C99 standard
.Ss Trigonometric and Hyperbolic Functions
.Cl
cacos arc cosine
cacosh arc hyperbolic cosine
casin arc sine
casinh arc hyperbolic sine
catan arc tangent
catanh arc hyperbolic tangent
ccos cosine
ccosh hyperbolic cosine
csin sine
@ -111,20 +117,8 @@ The
functions described here conform to
.St -isoC-99 .
.Sh BUGS
The inverse trigonometric and hyperbolic functions
.Fn cacos ,
.Fn cacosh ,
.Fn casin ,
.Fn casinh ,
.Fn catan ,
and
.Fn catanh
are not implemented.
.Pp
The logarithmic functions
.Fn clog
are not implemented.
.Pp
The power functions
and the power functions
.Fn cpow
are not implemented.

643
lib/msun/src/catrig.c Normal file
View File

@ -0,0 +1,643 @@
/*-
* Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <complex.h>
#include <float.h>
#include "math.h"
#include "math_private.h"
#undef isinf
#define isinf(x) (fabs(x) == INFINITY)
#undef isnan
#define isnan(x) ((x) != (x))
#define raise_inexact() do { volatile float junk = 1 + tiny; } while(0)
#undef signbit
#define signbit(x) (__builtin_signbit(x))
/* We need that DBL_EPSILON^2/128 is larger than FOUR_SQRT_MIN. */
static const double
A_crossover = 10, /* Hull et al suggest 1.5, but 10 works better */
B_crossover = 0.6417, /* suggested by Hull et al */
FOUR_SQRT_MIN = 0x1p-509, /* >= 4 * sqrt(DBL_MIN) */
QUARTER_SQRT_MAX = 0x1p509, /* <= sqrt(DBL_MAX) / 4 */
m_e = 2.7182818284590452e0, /* 0x15bf0a8b145769.0p-51 */
m_ln2 = 6.9314718055994531e-1, /* 0x162e42fefa39ef.0p-53 */
pio2_hi = 1.5707963267948966e0, /* 0x1921fb54442d18.0p-52 */
RECIP_EPSILON = 1 / DBL_EPSILON,
SQRT_3_EPSILON = 2.5809568279517849e-8, /* 0x1bb67ae8584caa.0p-78 */
SQRT_6_EPSILON = 3.6500241499888571e-8, /* 0x13988e1409212e.0p-77 */
SQRT_MIN = 0x1p-511; /* >= sqrt(DBL_MIN) */
static const volatile double
pio2_lo = 6.1232339957367659e-17; /* 0x11a62633145c07.0p-106 */
static const volatile float
tiny = 0x1p-100;
static double complex clog_for_large_values(double complex z);
/*
* Testing indicates that all these functions are accurate up to 4 ULP.
* The functions casin(h) and cacos(h) are about 2.5 times slower than asinh.
* The functions catan(h) are a little under 2 times slower than atanh.
*
* The code for casinh, casin, cacos, and cacosh comes first. The code is
* rather complicated, and the four functions are highly interdependent.
*
* The code for catanh and catan comes at the end. It is much simpler than
* the other functions, and the code for these can be disconnected from the
* rest of the code.
*/
/*
* ================================
* | casinh, casin, cacos, cacosh |
* ================================
*/
/*
* The algorithm is very close to that in "Implementing the complex arcsine
* and arccosine functions using exception handling" by T. E. Hull, Thomas F.
* Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
* Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
* http://dl.acm.org/citation.cfm?id=275324.
*
* Throughout we use the convention z = x + I*y.
*
* casinh(z) = sign(x)*log(A+sqrt(A*A-1)) + I*asin(B)
* where
* A = (|z+I| + |z-I|) / 2
* B = (|z+I| - |z-I|) / 2 = y/A
*
* These formulas become numerically unstable:
* (a) for Re(casinh(z)) when z is close to the line segment [-I, I] (that
* is, Re(casinh(z)) is close to 0);
* (b) for Im(casinh(z)) when z is close to either of the intervals
* [I, I*infinity) or (-I*infinity, -I] (that is, |Im(casinh(z))| is
* close to PI/2).
*
* These numerical problems are overcome by defining
* f(a, b) = (hypot(a, b) - b) / 2 = a*a / (hypot(a, b) + b) / 2
* Then if A < A_crossover, we use
* log(A + sqrt(A*A-1)) = log1p((A-1) + sqrt((A-1)*(A+1)))
* A-1 = f(x, 1+y) + f(x, 1-y)
* and if B > B_crossover, we use
* asin(B) = atan2(y, sqrt(A*A - y*y)) = atan2(y, sqrt((A+y)*(A-y)))
* A-y = f(x, y+1) + f(x, y-1)
* where without loss of generality we have assumed that x and y are
* non-negative.
*
* Much of the difficulty comes because the intermediate computations may
* produce overflows or underflows. This is dealt with in the paper by Hull
* et al by using exception handling. We do this by detecting when
* computations risk underflow or overflow. The hardest part is handling the
* underflows when computing f(a, b).
*
* Note that the function f(a, b) does not appear explicitly in the paper by
* Hull et al, but the idea may be found on pages 308 and 309. Introducing the
* function f(a, b) allows us to concentrate many of the clever tricks in this
* paper into one function.
*/
/*
* Function f(a, b, hypot_a_b) = (hypot(a, b) - b) / 2.
* Pass hypot(a, b) as the third argument.
*/
static inline double
f(double a, double b, double hypot_a_b)
{
if (b < 0)
return ((hypot_a_b - b) / 2);
if (b == 0)
return (a / 2);
return (a * a / (hypot_a_b + b) / 2);
}
/*
* All the hard work is contained in this function.
* x and y are assumed positive or zero, and less than RECIP_EPSILON.
* Upon return:
* rx = Re(casinh(z)) = -Im(cacos(y + I*x)).
* B_is_usable is set to 1 if the value of B is usable.
* If B_is_usable is set to 0, sqrt_A2my2 = sqrt(A*A - y*y), and new_y = y.
* If returning sqrt_A2my2 has potential to result in an underflow, it is
* rescaled, and new_y is similarly rescaled.
*/
static inline void
do_hard_work(double x, double y, double *rx, int *B_is_usable, double *B,
double *sqrt_A2my2, double *new_y)
{
double R, S, A; /* A, B, R, and S are as in Hull et al. */
double Am1, Amy; /* A-1, A-y. */
R = hypot(x, y + 1); /* |z+I| */
S = hypot(x, y - 1); /* |z-I| */
/* A = (|z+I| + |z-I|) / 2 */
A = (R + S) / 2;
/*
* Mathematically A >= 1. There is a small chance that this will not
* be so because of rounding errors. So we will make certain it is
* so.
*/
if (A < 1)
A = 1;
if (A < A_crossover) {
/*
* Am1 = fp + fm, where fp = f(x, 1+y), and fm = f(x, 1-y).
* rx = log1p(Am1 + sqrt(Am1*(A+1)))
*/
if (y == 1 && x < DBL_EPSILON*DBL_EPSILON / 128) {
/*
* fp is of order x^2, and fm = x/2.
* A = 1 (inexactly).
*/
*rx = sqrt(x);
} else if (x >= DBL_EPSILON * fabs(y - 1)) {
/*
* Underflow will not occur because
* x >= DBL_EPSILON^2/128 >= FOUR_SQRT_MIN
*/
Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
*rx = log1p(Am1 + sqrt(Am1 * (A + 1)));
} else if (y < 1) {
/*
* fp = x*x/(1+y)/4, fm = x*x/(1-y)/4, and
* A = 1 (inexactly).
*/
*rx = x / sqrt((1 - y) * (1 + y));
} else /* if (y > 1) */ {
/*
* A-1 = y-1 (inexactly).
*/
*rx = log1p((y - 1) + sqrt((y - 1) * (y + 1)));
}
} else {
*rx = log(A + sqrt(A * A - 1));
}
*new_y = y;
if (y < FOUR_SQRT_MIN) {
/*
* Avoid a possible underflow caused by y/A. For casinh this
* would be legitimate, but will be picked up by invoking atan2
* later on. For cacos this would not be legitimate.
*/
*B_is_usable = 0;
*sqrt_A2my2 = A * (2 / DBL_EPSILON);
*new_y = y * (2 / DBL_EPSILON);
return;
}
/* B = (|z+I| - |z-I|) / 2 = y/A */
*B = y / A;
*B_is_usable = 1;
if (*B > B_crossover) {
*B_is_usable = 0;
/*
* Amy = fp + fm, where fp = f(x, y+1), and fm = f(x, y-1).
* sqrt_A2my2 = sqrt(Amy*(A+y))
*/
if (y == 1 && x < DBL_EPSILON / 128) {
/*
* fp is of order x^2, and fm = x/2.
* A = 1 (inexactly).
*/
*sqrt_A2my2 = sqrt(x) * sqrt((A + y) / 2);
} else if (x >= DBL_EPSILON * fabs(y - 1)) {
/*
* Underflow will not occur because
* x >= DBL_EPSILON/128 >= FOUR_SQRT_MIN
* and
* x >= DBL_EPSILON^2 >= FOUR_SQRT_MIN
*/
Amy = f(x, y + 1, R) + f(x, y - 1, S);
*sqrt_A2my2 = sqrt(Amy * (A + y));
} else if (y > 1) {
/*
* fp = x*x/(y+1)/4, fm = x*x/(y-1)/4, and
* A = y (inexactly).
*
* y < RECIP_EPSILON. So the following
* scaling should avoid any underflow problems.
*/
*sqrt_A2my2 = x * (4 / DBL_EPSILON / DBL_EPSILON) * y /
sqrt((y + 1) * (y - 1));
*new_y = y * (4 / DBL_EPSILON / DBL_EPSILON);
} else /* if (y < 1) */ {
/*
* fm = 1-y >= DBL_EPSILON, fp is of order x^2, and
* A = 1 (inexactly).
*/
*sqrt_A2my2 = sqrt((1 - y) * (1 + y));
}
}
}
/*
* casinh(z) = z + O(z^3) as z -> 0
*
* casinh(z) = sign(x)*clog(sign(x)*z) + O(1/z^2) as z -> infinity
* The above formula works for the imaginary part as well, because
* Im(casinh(z)) = sign(x)*atan2(sign(x)*y, fabs(x)) + O(y/z^3)
* as z -> infinity, uniformly in y
*/
double complex
casinh(double complex z)
{
double x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
int B_is_usable;
double complex w;
x = creal(z);
y = cimag(z);
ax = fabs(x);
ay = fabs(y);
if (isnan(x) || isnan(y)) {
/* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */
if (isinf(x))
return (cpack(x, y + y));
/* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
if (isinf(y))
return (cpack(y, x + x));
/* casinh(NaN + I*0) = NaN + I*0 */
if (y == 0)
return (cpack(x + x, y));
/*
* All other cases involving NaN return NaN + I*NaN.
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
/* Bruce Evans tells me this is the way to do this: */
return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
/* clog...() will raise inexact unless x or y is infinite. */
if (signbit(x) == 0)
w = clog_for_large_values(z) + m_ln2;
else
w = clog_for_large_values(-z) + m_ln2;
return (cpack(copysign(creal(w), x), copysign(cimag(w), y)));
}
/* Avoid spuriously raising inexact for z = 0. */
if (x == 0 && y == 0)
return (z);
/* All remaining cases are inexact. */
raise_inexact();
if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
return (z);
do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
if (B_is_usable)
ry = asin(B);
else
ry = atan2(new_y, sqrt_A2my2);
return (cpack(copysign(rx, x), copysign(ry, y)));
}
/*
* casin(z) = reverse(casinh(reverse(z)))
* where reverse(x + I*y) = y + I*x = I*conj(z).
*/
double complex
casin(double complex z)
{
double complex w = casinh(cpack(cimag(z), creal(z)));
return (cpack(cimag(w), creal(w)));
}
/*
* cacos(z) = PI/2 - casin(z)
* but do the computation carefully so cacos(z) is accurate when z is
* close to 1.
*
* cacos(z) = PI/2 - z + O(z^3) as z -> 0
*
* cacos(z) = -sign(y)*I*clog(z) + O(1/z^2) as z -> infinity
* The above formula works for the real part as well, because
* Re(cacos(z)) = atan2(fabs(y), x) + O(y/z^3)
* as z -> infinity, uniformly in y
*/
double complex
cacos(double complex z)
{
double x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
int sx, sy;
int B_is_usable;
double complex w;
x = creal(z);
y = cimag(z);
sx = signbit(x);
sy = signbit(y);
ax = fabs(x);
ay = fabs(y);
if (isnan(x) || isnan(y)) {
/* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
if (isinf(x))
return (cpack(y + y, -INFINITY));
/* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
if (isinf(y))
return (cpack(x + x, -y));
/* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
if (x == 0)
return (cpack(pio2_hi + pio2_lo, y + y));
/*
* All other cases involving NaN return NaN + I*NaN.
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
/* clog...() will raise inexact unless x or y is infinite. */
w = clog_for_large_values(z);
rx = fabs(cimag(w));
ry = creal(w) + m_ln2;
if (sy == 0)
ry = -ry;
return (cpack(rx, ry));
}
/* Avoid spuriously raising inexact for z = 1. */
if (x == 1 && y == 0)
return (cpack(0, -y));
/* All remaining cases are inexact. */
raise_inexact();
if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON/4)
return (cpack(pio2_hi - (x - pio2_lo), -y));
do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
if (B_is_usable) {
if (sx==0)
rx = acos(B);
else
rx = acos(-B);
} else {
if (sx==0)
rx = atan2(sqrt_A2mx2, new_x);
else
rx = atan2(sqrt_A2mx2, -new_x);
}
if (sy == 0)
ry = -ry;
return (cpack(rx, ry));
}
/*
* cacosh(z) = I*cacos(z) or -I*cacos(z)
* where the sign is chosen so Re(cacosh(z)) >= 0.
*/
double complex
cacosh(double complex z)
{
double complex w;
double rx, ry;
w = cacos(z);
rx = creal(w);
ry = cimag(w);
/* cacosh(NaN + I*NaN) = NaN + I*NaN */
if (isnan(rx) && isnan(ry))
return (cpack(ry, rx));
/* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
/* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
if (isnan(rx))
return (cpack(fabs(ry), rx));
/* cacosh(0 + I*NaN) = NaN + I*NaN */
if (isnan(ry))
return (cpack(ry, ry));
return (cpack(fabs(ry), copysign(rx, cimag(z))));
}
/*
* Optimized version of clog() for |z| finite and larger than ~RECIP_EPSILON.
*/
static double complex
clog_for_large_values(double complex z)
{
double x, y;
double ax, ay, t;
x = creal(z);
y = cimag(z);
ax = fabs(x);
ay = fabs(y);
if (ax < ay) {
t = ax;
ax = ay;
ay = t;
}
/*
* Avoid overflow in hypot() when x and y are both very large.
* Divide x and y by E, and then add 1 to the logarithm. This depends
* on E being larger than sqrt(2).
* Dividing by E causes an insignificant loss of accuracy; however
* this method is still poor since it is uneccessarily slow.
*/
if (ax > DBL_MAX / 2)
return (cpack(log(hypot(x / m_e, y / m_e)) + 1, atan2(y, x)));
/*
* Avoid overflow when x or y is large. Avoid underflow when x or
* y is small.
*/
if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
return (cpack(log(hypot(x, y)), atan2(y, x)));
return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x)));
}
/*
*=============================================================================
*/
/*
* =================
* | catanh, catan |
* =================
*/
/*
* sum_squares(x,y) = x*x + y*y (or just x*x if y*y would underflow).
* Assumes x*x and y*y will not overflow.
* Assumes x and y are finite.
* Assumes y is non-negative.
* Assumes fabs(x) >= DBL_EPSILON.
*/
static inline double
sum_squares(double x, double y)
{
/* Avoid underflow when y is small. */
if (y < SQRT_MIN)
return (x * x);
return (x * x + y * y);
}
/*
* real_part_reciprocal(x, y) = Re(1/(x+I*y)) = x/(x*x + y*y).
* Assumes x and y are not NaN, and one of x and y is larger than
* RECIP_EPSILON. We avoid unwarranted underflow. It is important to not use
* the code creal(1/z), because the imaginary part may produce an unwanted
* underflow.
* This is only called in a context where inexact is always raised before
* the call, so no effort is made to avoid or force inexact.
*/
static inline double
real_part_reciprocal(double x, double y)
{
double scale;
uint32_t hx, hy;
int32_t ix, iy;
/*
* This code is inspired by the C99 document n1124.pdf, Section G.5.1,
* example 2.
*/
GET_HIGH_WORD(hx, x);
ix = hx & 0x7ff00000;
GET_HIGH_WORD(hy, y);
iy = hy & 0x7ff00000;
#define BIAS (DBL_MAX_EXP - 1)
/* XXX more guard digits are useful iff there is extra precision. */
#define CUTOFF (DBL_MANT_DIG / 2 + 1) /* just half or 1 guard digit */
if (ix - iy >= CUTOFF << 20 || isinf(x))
return (1 / x); /* +-Inf -> +-0 is special */
if (iy - ix >= CUTOFF << 20)
return (x / y / y); /* should avoid double div, but hard */
if (ix <= (BIAS + DBL_MAX_EXP / 2 - CUTOFF) << 20)
return (x / (x * x + y * y));
scale = 1;
SET_HIGH_WORD(scale, 0x7ff00000 - ix); /* 2**(1-ilogb(x)) */
x *= scale;
y *= scale;
return (x / (x * x + y * y) * scale);
}
/*
* catanh(z) = log((1+z)/(1-z)) / 2
* = log1p(4*x / |z-1|^2) / 4
* + I * atan2(2*y, (1-x)*(1+x)-y*y) / 2
*
* catanh(z) = z + O(z^3) as z -> 0
*
* catanh(z) = 1/z + sign(y)*I*PI/2 + O(1/z^3) as z -> infinity
* The above formula works for the real part as well, because
* Re(catanh(z)) = x/|z|^2 + O(x/z^4)
* as z -> infinity, uniformly in x
*/
double complex
catanh(double complex z)
{
double x, y, ax, ay, rx, ry;
x = creal(z);
y = cimag(z);
ax = fabs(x);
ay = fabs(y);
/* This helps handle many cases. */
if (y == 0 && ax <= 1)
return (cpack(atanh(x), y));
/* To ensure the same accuracy as atan(), and to filter out z = 0. */
if (x == 0)
return (cpack(x, atan(y)));
if (isnan(x) || isnan(y)) {
/* catanh(+-Inf + I*NaN) = +-0 + I*NaN */
if (isinf(x))
return (cpack(copysign(0, x), y + y));
/* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
if (isinf(y)) {
return (cpack(copysign(0, x),
copysign(pio2_hi + pio2_lo, y)));
}
/*
* All other cases involving NaN return NaN + I*NaN.
* C99 leaves it optional whether to raise invalid if one of
* the arguments is not NaN, so we opt not to raise it.
*/
return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
return (cpack(real_part_reciprocal(x, y),
copysign(pio2_hi + pio2_lo, y)));
}
if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
/*
* z = 0 was filtered out above. All other cases must raise
* inexact, but this is the only only that needs to do it
* explicitly.
*/
raise_inexact();
return (z);
}
if (ax == 1 && ay < DBL_EPSILON)
rx = (log(ay) - m_ln2) / -2;
else
rx = log1p(4 * ax / sum_squares(ax - 1, ay)) / 4;
if (ax == 1)
ry = atan2(2, -ay) / 2;
else if (ay < DBL_EPSILON)
ry = atan2(2 * ay, (1 - ax) * (1 + ax)) / 2;
else
ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
return (cpack(copysign(rx, x), copysign(ry, y)));
}
/*
* catan(z) = reverse(catanh(reverse(z)))
* where reverse(x + I*y) = y + I*x = I*conj(z).
*/
double complex
catan(double complex z)
{
double complex w = catanh(cpack(cimag(z), creal(z)));
return (cpack(cimag(w), creal(w)));
}

388
lib/msun/src/catrigf.c Normal file
View File

@ -0,0 +1,388 @@
/*-
* Copyright (c) 2012 Stephen Montgomery-Smith <stephen@FreeBSD.ORG>
* All rights reserved.
*
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions
* are met:
* 1. Redistributions of source code must retain the above copyright
* notice, this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright
* notice, this list of conditions and the following disclaimer in the
* documentation and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
* ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
* FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
* DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
* OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
* HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
* LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*/
/*
* The algorithm is very close to that in "Implementing the complex arcsine
* and arccosine functions using exception handling" by T. E. Hull, Thomas F.
* Fairgrieve, and Ping Tak Peter Tang, published in ACM Transactions on
* Mathematical Software, Volume 23 Issue 3, 1997, Pages 299-335,
* http://dl.acm.org/citation.cfm?id=275324.
*
* The code for catrig.c contains complete comments.
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
#include <complex.h>
#include <float.h>
#include "math.h"
#include "math_private.h"
#undef isinf
#define isinf(x) (fabsf(x) == INFINITY)
#undef isnan
#define isnan(x) ((x) != (x))
#define raise_inexact() do { volatile float junk = 1 + tiny; } while(0)
#undef signbit
#define signbit(x) (__builtin_signbitf(x))
static const float
A_crossover = 10,
B_crossover = 0.6417,
FOUR_SQRT_MIN = 0x1p-61,
QUARTER_SQRT_MAX = 0x1p61,
m_e = 2.7182818285e0, /* 0xadf854.0p-22 */
m_ln2 = 6.9314718056e-1, /* 0xb17218.0p-24 */
pio2_hi = 1.5707962513e0, /* 0xc90fda.0p-23 */
RECIP_EPSILON = 1 / FLT_EPSILON,
SQRT_3_EPSILON = 5.9801995673e-4, /* 0x9cc471.0p-34 */
SQRT_6_EPSILON = 8.4572793338e-4, /* 0xddb3d7.0p-34 */
SQRT_MIN = 0x1p-63;
static const volatile float
pio2_lo = 7.5497899549e-8, /* 0xa22169.0p-47 */
tiny = 0x1p-100;
static float complex clog_for_large_values(float complex z);
static inline float
f(float a, float b, float hypot_a_b)
{
if (b < 0)
return ((hypot_a_b - b) / 2);
if (b == 0)
return (a / 2);
return (a * a / (hypot_a_b + b) / 2);
}
static inline void
do_hard_work(float x, float y, float *rx, int *B_is_usable, float *B,
float *sqrt_A2my2, float *new_y)
{
float R, S, A;
float Am1, Amy;
R = hypotf(x, y + 1);
S = hypotf(x, y - 1);
A = (R + S) / 2;
if (A < 1)
A = 1;
if (A < A_crossover) {
if (y == 1 && x < FLT_EPSILON * FLT_EPSILON / 128) {
*rx = sqrtf(x);
} else if (x >= FLT_EPSILON * fabsf(y - 1)) {
Am1 = f(x, 1 + y, R) + f(x, 1 - y, S);
*rx = log1pf(Am1 + sqrtf(Am1 * (A + 1)));
} else if (y < 1) {
*rx = x / sqrtf((1 - y)*(1 + y));
} else {
*rx = log1pf((y - 1) + sqrtf((y - 1) * (y + 1)));
}
} else {
*rx = logf(A + sqrtf(A * A - 1));
}
*new_y = y;
if (y < FOUR_SQRT_MIN) {
*B_is_usable = 0;
*sqrt_A2my2 = A * (2 / FLT_EPSILON);
*new_y = y * (2 / FLT_EPSILON);
return;
}
*B = y / A;
*B_is_usable = 1;
if (*B > B_crossover) {
*B_is_usable = 0;
if (y == 1 && x < FLT_EPSILON / 128) {
*sqrt_A2my2 = sqrtf(x) * sqrtf((A + y) / 2);
} else if (x >= FLT_EPSILON * fabsf(y - 1)) {
Amy = f(x, y + 1, R) + f(x, y - 1, S);
*sqrt_A2my2 = sqrtf(Amy * (A + y));
} else if (y > 1) {
*sqrt_A2my2 = x * (4 / FLT_EPSILON / FLT_EPSILON) * y /
sqrtf((y + 1) * (y - 1));
*new_y = y * (4 / FLT_EPSILON / FLT_EPSILON);
} else {
*sqrt_A2my2 = sqrtf((1 - y) * (1 + y));
}
}
}
float complex
casinhf(float complex z)
{
float x, y, ax, ay, rx, ry, B, sqrt_A2my2, new_y;
int B_is_usable;
float complex w;
x = crealf(z);
y = cimagf(z);
ax = fabsf(x);
ay = fabsf(y);
if (isnan(x) || isnan(y)) {
if (isinf(x))
return (cpackf(x, y + y));
if (isinf(y))
return (cpackf(y, x + x));
if (y == 0)
return (cpackf(x + x, y));
return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
if (signbit(x) == 0)
w = clog_for_large_values(z) + m_ln2;
else
w = clog_for_large_values(-z) + m_ln2;
return (cpackf(copysignf(crealf(w), x),
copysignf(cimagf(w), y)));
}
if (x == 0 && y == 0)
return (z);
raise_inexact();
if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
return (z);
do_hard_work(ax, ay, &rx, &B_is_usable, &B, &sqrt_A2my2, &new_y);
if (B_is_usable)
ry = asinf(B);
else
ry = atan2f(new_y, sqrt_A2my2);
return (cpackf(copysignf(rx, x), copysignf(ry, y)));
}
float complex
casinf(float complex z)
{
float complex w = casinhf(cpackf(cimagf(z), crealf(z)));
return (cpackf(cimagf(w), crealf(w)));
}
float complex
cacosf(float complex z)
{
float x, y, ax, ay, rx, ry, B, sqrt_A2mx2, new_x;
int sx, sy;
int B_is_usable;
float complex w;
x = crealf(z);
y = cimagf(z);
sx = signbit(x);
sy = signbit(y);
ax = fabsf(x);
ay = fabsf(y);
if (isnan(x) || isnan(y)) {
if (isinf(x))
return (cpackf(y + y, -INFINITY));
if (isinf(y))
return (cpackf(x + x, -y));
if (x == 0) return (cpackf(pio2_hi + pio2_lo, y + y));
return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
w = clog_for_large_values(z);
rx = fabsf(cimagf(w));
ry = crealf(w) + m_ln2;
if (sy == 0)
ry = -ry;
return (cpackf(rx, ry));
}
if (x == 1 && y == 0)
return (cpackf(0, -y));
raise_inexact();
if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
return (cpackf(pio2_hi - (x - pio2_lo), -y));
do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
if (B_is_usable) {
if (sx==0)
rx = acosf(B);
else
rx = acosf(-B);
} else {
if (sx==0)
rx = atan2f(sqrt_A2mx2, new_x);
else
rx = atan2f(sqrt_A2mx2, -new_x);
}
if (sy==0)
ry = -ry;
return (cpackf(rx, ry));
}
float complex
cacoshf(float complex z)
{
float complex w;
float rx, ry;
w = cacosf(z);
rx = crealf(w);
ry = cimagf(w);
if (isnan(rx) && isnan(ry))
return (cpackf(ry, rx));
if (isnan(rx))
return (cpackf(fabsf(ry), rx));
if (isnan(ry))
return (cpackf(ry, ry));
return (cpackf(fabsf(ry), copysignf(rx, cimagf(z))));
}
static float complex
clog_for_large_values(float complex z)
{
float x, y;
float ax, ay, t;
x = crealf(z);
y = cimagf(z);
ax = fabsf(x);
ay = fabsf(y);
if (ax < ay) {
t = ax;
ax = ay;
ay = t;
}
if (ax > FLT_MAX / 2) {
return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1,
atan2f(y, x)));
}
if (ax > QUARTER_SQRT_MAX || ay < SQRT_MIN)
return (cpackf(logf(hypotf(x, y)), atan2f(y, x)));
return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
}
static inline float
sum_squares(float x, float y)
{
if (y < SQRT_MIN)
return (x*x);
return (x*x + y*y);
}
static inline float
real_part_reciprocal(float x, float y)
{
float scale;
uint32_t hx, hy;
int32_t ix, iy;
GET_FLOAT_WORD(hx, x);
ix = hx & 0x7f800000;
GET_FLOAT_WORD(hy, y);
iy = hy & 0x7f800000;
#define BIAS (FLT_MAX_EXP - 1)
#define CUTOFF (FLT_MANT_DIG / 2 + 1)
if (ix - iy >= CUTOFF << 23 || isinf(x))
return (1/x);
if (iy - ix >= CUTOFF << 23)
return (x/y/y);
if (ix <= (BIAS + FLT_MAX_EXP / 2 - CUTOFF) << 23)
return (x / (x * x + y * y));
SET_FLOAT_WORD(scale, 0x7f800000 - ix);
x *= scale;
y *= scale;
return (x / (x * x + y * y) * scale);
}
float complex
catanhf(float complex z)
{
float x, y, ax, ay, rx, ry;
x = crealf(z);
y = cimagf(z);
ax = fabsf(x);
ay = fabsf(y);
if (y == 0 && ax <= 1)
return (cpackf(atanhf(x), y));
if (x == 0)
return (cpackf(x, atanf(y)));
if (isnan(x) || isnan(y)) {
if (isinf(x))
return (cpackf(copysignf(0, x), y+y));
if (isinf(y)) {
return (cpackf(copysignf(0, x),
copysignf(pio2_hi + pio2_lo, y)));
}
return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
return (cpackf(real_part_reciprocal(x, y),
copysignf(pio2_hi + pio2_lo, y)));
}
if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
raise_inexact();
return (z);
}
if (ax == 1 && ay < FLT_EPSILON)
rx = (logf(ay) - m_ln2) / -2;
else
rx = log1pf(4 * ax / sum_squares(ax - 1, ay)) / 4;
if (ax == 1)
ry = atan2f(2, -ay) / 2;
else if (ay < FLT_EPSILON)
ry = atan2f(2 * ay, (1 - ax) * (1 + ax)) / 2;
else
ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
return (cpackf(copysignf(rx, x), copysignf(ry, y)));
}
float complex
catanf(float complex z)
{
float complex w = catanhf(cpackf(cimagf(z), crealf(z)));
return (cpackf(cimagf(w), crealf(w)));
}