generic-library/bionic
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Merge "Sync with upstream FreeBSD libm."

Browse Source
This commit is contained in:
Elliott Hughes 2015-08-29 03:46:21 +00:00 committed by Gerrit Code Review
commit cb58f528bf
31 changed files with 408 additions and 379 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
libm/upstream-freebsd/lib/msun/ld128/k_expl.h
View File

@ -29,7 +29,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/ld128/k_expl.h 275819 2014-12-16 09:21:56Z ed $");
/*
* ld128 version of k_expl.h. See ../ld80/s_expl.c for most comments.
@ -322,7 +322,7 @@ __ldexp_cexpl(long double complex z, int expt)
scale2 = 1;
SET_LDBL_EXPSIGN(scale1, BIAS + expt - half_expt);
return (cpackl(cos(y) * exp_x * scale1 * scale2,
return (CMPLXL(cos(y) * exp_x * scale1 * scale2,
sinl(y) * exp_x * scale1 * scale2));
}
#endif /* _COMPLEX_H */

66
libm/upstream-freebsd/lib/msun/src/catrig.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/catrig.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <float.h>
@ -286,19 +286,19 @@ casinh(double complex z)
if (isnan(x) || isnan(y)) {
/* casinh(+-Inf + I*NaN) = +-Inf + I*NaN */
if (isinf(x))
return (cpack(x, y + y));
return (CMPLX(x, y + y));
/* casinh(NaN + I*+-Inf) = opt(+-)Inf + I*NaN */
if (isinf(y))
return (cpack(y, x + x));
return (CMPLX(y, x + x));
/* casinh(NaN + I*0) = NaN + I*0 */
if (y == 0)
return (cpack(x + x, y));
return (CMPLX(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.
*/
return (cpack(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -307,7 +307,7 @@ casinh(double complex z)
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)));
return (CMPLX(copysign(creal(w), x), copysign(cimag(w), y)));
}
/* Avoid spuriously raising inexact for z = 0. */
@ -325,7 +325,7 @@ casinh(double complex z)
ry = asin(B);
else
ry = atan2(new_y, sqrt_A2my2);
return (cpack(copysign(rx, x), copysign(ry, y)));
return (CMPLX(copysign(rx, x), copysign(ry, y)));
}
/*
@ -335,9 +335,9 @@ casinh(double complex z)
double complex
casin(double complex z)
{
double complex w = casinh(cpack(cimag(z), creal(z)));
double complex w = casinh(CMPLX(cimag(z), creal(z)));
return (cpack(cimag(w), creal(w)));
return (CMPLX(cimag(w), creal(w)));
}
/*
@ -370,19 +370,19 @@ cacos(double complex z)
if (isnan(x) || isnan(y)) {
/* cacos(+-Inf + I*NaN) = NaN + I*opt(-)Inf */
if (isinf(x))
return (cpack(y + y, -INFINITY));
return (CMPLX(y + y, -INFINITY));
/* cacos(NaN + I*+-Inf) = NaN + I*-+Inf */
if (isinf(y))
return (cpack(x + x, -y));
return (CMPLX(x + x, -y));
/* cacos(0 + I*NaN) = PI/2 + I*NaN with inexact */
if (x == 0)
return (cpack(pio2_hi + pio2_lo, y + y));
return (CMPLX(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)));
return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -392,18 +392,18 @@ cacos(double complex z)
ry = creal(w) + m_ln2;
if (sy == 0)
ry = -ry;
return (cpack(rx, ry));
return (CMPLX(rx, ry));
}
/* Avoid spuriously raising inexact for z = 1. */
if (x == 1 && y == 0)
return (cpack(0, -y));
return (CMPLX(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));
return (CMPLX(pio2_hi - (x - pio2_lo), -y));
do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
if (B_is_usable) {
@ -419,7 +419,7 @@ cacos(double complex z)
}
if (sy == 0)
ry = -ry;
return (cpack(rx, ry));
return (CMPLX(rx, ry));
}
/*
@ -437,15 +437,15 @@ cacosh(double complex z)
ry = cimag(w);
/* cacosh(NaN + I*NaN) = NaN + I*NaN */
if (isnan(rx) && isnan(ry))
return (cpack(ry, rx));
return (CMPLX(ry, rx));
/* cacosh(NaN + I*+-Inf) = +Inf + I*NaN */
/* cacosh(+-Inf + I*NaN) = +Inf + I*NaN */
if (isnan(rx))
return (cpack(fabs(ry), rx));
return (CMPLX(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))));
return (CMPLX(ry, ry));
return (CMPLX(fabs(ry), copysign(rx, cimag(z))));
}
/*
@ -475,16 +475,16 @@ clog_for_large_values(double complex z)
* 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)));
return (CMPLX(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 (CMPLX(log(hypot(x, y)), atan2(y, x)));
return (cpack(log(ax * ax + ay * ay) / 2, atan2(y, x)));
return (CMPLX(log(ax * ax + ay * ay) / 2, atan2(y, x)));
}
/*
@ -575,30 +575,30 @@ catanh(double complex z)
/* This helps handle many cases. */
if (y == 0 && ax <= 1)
return (cpack(atanh(x), y));
return (CMPLX(atanh(x), y));
/* To ensure the same accuracy as atan(), and to filter out z = 0. */
if (x == 0)
return (cpack(x, atan(y)));
return (CMPLX(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));
return (CMPLX(copysign(0, x), y + y));
/* catanh(NaN + I*+-Inf) = sign(NaN)0 + I*+-PI/2 */
if (isinf(y))
return (cpack(copysign(0, x),
return (CMPLX(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)));
return (CMPLX(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
return (cpack(real_part_reciprocal(x, y),
return (CMPLX(real_part_reciprocal(x, y),
copysign(pio2_hi + pio2_lo, y)));
if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
@ -623,7 +623,7 @@ catanh(double complex z)
else
ry = atan2(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
return (cpack(copysign(rx, x), copysign(ry, y)));
return (CMPLX(copysign(rx, x), copysign(ry, y)));
}
/*
@ -633,7 +633,7 @@ catanh(double complex z)
double complex
catan(double complex z)
{
double complex w = catanh(cpack(cimag(z), creal(z)));
double complex w = catanh(CMPLX(cimag(z), creal(z)));
return (cpack(cimag(w), creal(w)));
return (CMPLX(cimag(w), creal(w)));
}

66
libm/upstream-freebsd/lib/msun/src/catrigf.c
View File

@ -39,7 +39,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/catrigf.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <float.h>
@ -156,12 +156,12 @@ casinhf(float complex z)
if (isnan(x) || isnan(y)) {
if (isinf(x))
return (cpackf(x, y + y));
return (CMPLXF(x, y + y));
if (isinf(y))
return (cpackf(y, x + x));
return (CMPLXF(y, x + x));
if (y == 0)
return (cpackf(x + x, y));
return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXF(x + x, y));
return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -169,7 +169,7 @@ casinhf(float complex z)
w = clog_for_large_values(z) + m_ln2;
else
w = clog_for_large_values(-z) + m_ln2;
return (cpackf(copysignf(crealf(w), x),
return (CMPLXF(copysignf(crealf(w), x),
copysignf(cimagf(w), y)));
}
@ -186,15 +186,15 @@ casinhf(float complex z)
ry = asinf(B);
else
ry = atan2f(new_y, sqrt_A2my2);
return (cpackf(copysignf(rx, x), copysignf(ry, y)));
return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
}
float complex
casinf(float complex z)
{
float complex w = casinhf(cpackf(cimagf(z), crealf(z)));
float complex w = casinhf(CMPLXF(cimagf(z), crealf(z)));
return (cpackf(cimagf(w), crealf(w)));
return (CMPLXF(cimagf(w), crealf(w)));
}
float complex
@ -214,12 +214,12 @@ cacosf(float complex z)
if (isnan(x) || isnan(y)) {
if (isinf(x))
return (cpackf(y + y, -INFINITY));
return (CMPLXF(y + y, -INFINITY));
if (isinf(y))
return (cpackf(x + x, -y));
return (CMPLXF(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)));
return (CMPLXF(pio2_hi + pio2_lo, y + y));
return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON) {
@ -228,16 +228,16 @@ cacosf(float complex z)
ry = crealf(w) + m_ln2;
if (sy == 0)
ry = -ry;
return (cpackf(rx, ry));
return (CMPLXF(rx, ry));
}
if (x == 1 && y == 0)
return (cpackf(0, -y));
return (CMPLXF(0, -y));
raise_inexact();
if (ax < SQRT_6_EPSILON / 4 && ay < SQRT_6_EPSILON / 4)
return (cpackf(pio2_hi - (x - pio2_lo), -y));
return (CMPLXF(pio2_hi - (x - pio2_lo), -y));
do_hard_work(ay, ax, &ry, &B_is_usable, &B, &sqrt_A2mx2, &new_x);
if (B_is_usable) {
@ -253,7 +253,7 @@ cacosf(float complex z)
}
if (sy == 0)
ry = -ry;
return (cpackf(rx, ry));
return (CMPLXF(rx, ry));
}
float complex
@ -266,12 +266,12 @@ cacoshf(float complex z)
rx = crealf(w);
ry = cimagf(w);
if (isnan(rx) && isnan(ry))
return (cpackf(ry, rx));
return (CMPLXF(ry, rx));
if (isnan(rx))
return (cpackf(fabsf(ry), rx));
return (CMPLXF(fabsf(ry), rx));
if (isnan(ry))
return (cpackf(ry, ry));
return (cpackf(fabsf(ry), copysignf(rx, cimagf(z))));
return (CMPLXF(ry, ry));
return (CMPLXF(fabsf(ry), copysignf(rx, cimagf(z))));
}
static float complex
@ -291,13 +291,13 @@ clog_for_large_values(float complex z)
}
if (ax > FLT_MAX / 2)
return (cpackf(logf(hypotf(x / m_e, y / m_e)) + 1,
return (CMPLXF(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 (CMPLXF(logf(hypotf(x, y)), atan2f(y, x)));
return (cpackf(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
return (CMPLXF(logf(ax * ax + ay * ay) / 2, atan2f(y, x)));
}
static inline float
@ -346,22 +346,22 @@ catanhf(float complex z)
ay = fabsf(y);
if (y == 0 && ax <= 1)
return (cpackf(atanhf(x), y));
return (CMPLXF(atanhf(x), y));
if (x == 0)
return (cpackf(x, atanf(y)));
return (CMPLXF(x, atanf(y)));
if (isnan(x) || isnan(y)) {
if (isinf(x))
return (cpackf(copysignf(0, x), y + y));
return (CMPLXF(copysignf(0, x), y + y));
if (isinf(y))
return (cpackf(copysignf(0, x),
return (CMPLXF(copysignf(0, x),
copysignf(pio2_hi + pio2_lo, y)));
return (cpackf(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
return (CMPLXF(x + 0.0L + (y + 0), x + 0.0L + (y + 0)));
}
if (ax > RECIP_EPSILON || ay > RECIP_EPSILON)
return (cpackf(real_part_reciprocal(x, y),
return (CMPLXF(real_part_reciprocal(x, y),
copysignf(pio2_hi + pio2_lo, y)));
if (ax < SQRT_3_EPSILON / 2 && ay < SQRT_3_EPSILON / 2) {
@ -381,13 +381,13 @@ catanhf(float complex z)
else
ry = atan2f(2 * ay, (1 - ax) * (1 + ax) - ay * ay) / 2;
return (cpackf(copysignf(rx, x), copysignf(ry, y)));
return (CMPLXF(copysignf(rx, x), copysignf(ry, y)));
}
float complex
catanf(float complex z)
{
float complex w = catanhf(cpackf(cimagf(z), crealf(z)));
float complex w = catanhf(CMPLXF(cimagf(z), crealf(z)));
return (cpackf(cimagf(w), crealf(w)));
return (CMPLXF(cimagf(w), crealf(w)));
}

32
libm/upstream-freebsd/lib/msun/src/e_j0.c
View File

@ -12,7 +12,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/e_j0.c 283032 2015-05-17 16:27:06Z kargl $");
/* __ieee754_j0(x), __ieee754_y0(x)
* Bessel function of the first and second kinds of order zero.
@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
static double pzero(double), qzero(double);
static __inline double pzero(double), qzero(double);
static const volatile double vone = 1, vzero = 0;
static const double
huge = 1e300,
@ -115,7 +117,7 @@ __ieee754_j0(double x)
if(ix<0x3f200000) { /* |x| < 2**-13 */
if(huge+x>one) { /* raise inexact if x != 0 */
if(ix<0x3e400000) return one; /* |x|<2**-27 */
else return one - 0.25*x*x;
else return one - x*x/4;
}
}
z = x*x;
@ -150,10 +152,16 @@ __ieee754_y0(double x)
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if(ix>=0x7ff00000) return one/(x+x*x);
if((ix|lx)==0) return -one/zero;
if(hx<0) return zero/zero;
/*
* y0(NaN) = NaN.
* y0(Inf) = 0.
* y0(-Inf) = NaN and raise invalid exception.
*/
if(ix>=0x7ff00000) return vone/(x+x*x);
/* y0(+-0) = -inf and raise divide-by-zero exception. */
if((ix|lx)==0) return -one/vzero;
/* y0(x<0) = NaN and raise invalid exception. */
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
@ -268,7 +276,8 @@ static const double pS2[5] = {
1.46576176948256193810e+01, /* 0x402D50B3, 0x44391809 */
};
static double pzero(double x)
static __inline double
pzero(double x)
{
const double *p,*q;
double z,r,s;
@ -278,7 +287,7 @@ static const double pS2[5] = {
if(ix>=0x40200000) {p = pR8; q= pS8;}
else if(ix>=0x40122E8B){p = pR5; q= pS5;}
else if(ix>=0x4006DB6D){p = pR3; q= pS3;}
else if(ix>=0x40000000){p = pR2; q= pS2;}
else {p = pR2; q= pS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -363,7 +372,8 @@ static const double qS2[6] = {
-5.31095493882666946917e+00, /* 0xC0153E6A, 0xF8B32931 */
};
static double qzero(double x)
static __inline double
qzero(double x)
{
const double *p,*q;
double s,r,z;
@ -373,7 +383,7 @@ static const double qS2[6] = {
if(ix>=0x40200000) {p = qR8; q= qS8;}
else if(ix>=0x40122E8B){p = qR5; q= qS5;}
else if(ix>=0x4006DB6D){p = qR3; q= qS3;}
else if(ix>=0x40000000){p = qR2; q= qS2;}
else {p = qR2; q= qS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

47
libm/upstream-freebsd/lib/msun/src/e_j0f.c
View File

@ -14,12 +14,18 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/e_j0f.c 283032 2015-05-17 16:27:06Z kargl $");
/*
* See e_j0.c for complete comments.
*/
#include "math.h"
#include "math_private.h"
static float pzerof(float), qzerof(float);
static __inline float pzerof(float), qzerof(float);
static const volatile float vone = 1, vzero = 0;
static const float
huge = 1e30,
@ -62,17 +68,17 @@ __ieee754_j0f(float x)
* j0(x) = 1/sqrt(pi) * (P(0,x)*cc - Q(0,x)*ss) / sqrt(x)
* y0(x) = 1/sqrt(pi) * (P(0,x)*ss + Q(0,x)*cc) / sqrt(x)
*/
if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(x);
if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(x); /* |x|>2**49 */
else {
u = pzerof(x); v = qzerof(x);
z = invsqrtpi*(u*cc-v*ss)/sqrtf(x);
}
return z;
}
if(ix<0x39000000) { /* |x| < 2**-13 */
if(ix<0x3b000000) { /* |x| < 2**-9 */
if(huge+x>one) { /* raise inexact if x != 0 */
if(ix<0x32000000) return one; /* |x|<2**-27 */
else return one - (float)0.25*x*x;
if(ix<0x39800000) return one; /* |x|<2**-12 */
else return one - x*x/4;
}
}
z = x*x;
@ -107,10 +113,9 @@ __ieee754_y0f(float x)
GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx;
/* Y0(NaN) is NaN, y0(-inf) is Nan, y0(inf) is 0 */
if(ix>=0x7f800000) return one/(x+x*x);
if(ix==0) return -one/zero;
if(hx<0) return zero/zero;
if(ix>=0x7f800000) return vone/(x+x*x);
if(ix==0) return -one/vzero;
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
/* y0(x) = sqrt(2/(pi*x))*(p0(x)*sin(x0)+q0(x)*cos(x0))
* where x0 = x-pi/4
@ -136,14 +141,14 @@ __ieee754_y0f(float x)
if ((s*c)<zero) cc = z/ss;
else ss = z/cc;
}
if(ix>0x80000000) z = (invsqrtpi*ss)/sqrtf(x);
if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
else {
u = pzerof(x); v = qzerof(x);
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
}
return z;
}
if(ix<=0x32000000) { /* x < 2**-27 */
if(ix<=0x39000000) { /* x < 2**-13 */
return(u00 + tpi*__ieee754_logf(x));
}
z = x*x;
@ -224,7 +229,8 @@ static const float pS2[5] = {
1.4657617569e+01, /* 0x416a859a */
};
static float pzerof(float x)
static __inline float
pzerof(float x)
{
const float *p,*q;
float z,r,s;
@ -232,9 +238,9 @@ static const float pS2[5] = {
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x41000000) {p = pR8; q= pS8;}
else if(ix>=0x40f71c58){p = pR5; q= pS5;}
else if(ix>=0x4036db68){p = pR3; q= pS3;}
else if(ix>=0x40000000){p = pR2; q= pS2;}
else if(ix>=0x409173eb){p = pR5; q= pS5;}
else if(ix>=0x4036d917){p = pR3; q= pS3;}
else {p = pR2; q= pS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -319,7 +325,8 @@ static const float qS2[6] = {
-5.3109550476e+00, /* 0xc0a9f358 */
};
static float qzerof(float x)
static __inline float
qzerof(float x)
{
const float *p,*q;
float s,r,z;
@ -327,9 +334,9 @@ static const float qS2[6] = {
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x41000000) {p = qR8; q= qS8;}
else if(ix>=0x40f71c58){p = qR5; q= qS5;}
else if(ix>=0x4036db68){p = qR3; q= qS3;}
else if(ix>=0x40000000){p = qR2; q= qS2;}
else if(ix>=0x409173eb){p = qR5; q= qS5;}
else if(ix>=0x4036d917){p = qR3; q= qS3;}
else {p = qR2; q= qS2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

30
libm/upstream-freebsd/lib/msun/src/e_j1.c
View File

@ -12,7 +12,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/e_j1.c 283032 2015-05-17 16:27:06Z kargl $");
/* __ieee754_j1(x), __ieee754_y1(x)
* Bessel function of the first and second kinds of order zero.
@ -62,7 +62,9 @@ __FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
static double pone(double), qone(double);
static __inline double pone(double), qone(double);
static const volatile double vone = 1, vzero = 0;
static const double
huge = 1e300,
@ -147,10 +149,16 @@ __ieee754_y1(double x)
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if(ix>=0x7ff00000) return one/(x+x*x);
if((ix|lx)==0) return -one/zero;
if(hx<0) return zero/zero;
/*
* y1(NaN) = NaN.
* y1(Inf) = 0.
* y1(-Inf) = NaN and raise invalid exception.
*/
if(ix>=0x7ff00000) return vone/(x+x*x);
/* y1(+-0) = -inf and raise divide-by-zero exception. */
if((ix|lx)==0) return -one/vzero;
/* y1(x<0) = NaN and raise invalid exception. */
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sin(x);
c = cos(x);
@ -262,7 +270,8 @@ static const double ps2[5] = {
8.36463893371618283368e+00, /* 0x4020BAB1, 0xF44E5192 */
};
static double pone(double x)
static __inline double
pone(double x)
{
const double *p,*q;
double z,r,s;
@ -272,7 +281,7 @@ static const double ps2[5] = {
if(ix>=0x40200000) {p = pr8; q= ps8;}
else if(ix>=0x40122E8B){p = pr5; q= ps5;}
else if(ix>=0x4006DB6D){p = pr3; q= ps3;}
else if(ix>=0x40000000){p = pr2; q= ps2;}
else {p = pr2; q= ps2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -358,7 +367,8 @@ static const double qs2[6] = {
-4.95949898822628210127e+00, /* 0xC013D686, 0xE71BE86B */
};
static double qone(double x)
static __inline double
qone(double x)
{
const double *p,*q;
double s,r,z;
@ -368,7 +378,7 @@ static const double qs2[6] = {
if(ix>=0x40200000) {p = qr8; q= qs8;}
else if(ix>=0x40122E8B){p = qr5; q= qs5;}
else if(ix>=0x4006DB6D){p = qr3; q= qs3;}
else if(ix>=0x40000000){p = qr2; q= qs2;}
else {p = qr2; q= qs2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

45
libm/upstream-freebsd/lib/msun/src/e_j1f.c
View File

@ -14,12 +14,18 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/e_j1f.c 283032 2015-05-17 16:27:06Z kargl $");
/*
* See e_j1.c for complete comments.
*/
#include "math.h"
#include "math_private.h"
static float ponef(float), qonef(float);
static __inline float ponef(float), qonef(float);
static const volatile float vone = 1, vzero = 0;
static const float
huge = 1e30,
@ -63,7 +69,7 @@ __ieee754_j1f(float x)
* j1(x) = 1/sqrt(pi) * (P(1,x)*cc - Q(1,x)*ss) / sqrt(x)
* y1(x) = 1/sqrt(pi) * (P(1,x)*ss + Q(1,x)*cc) / sqrt(x)
*/
if(ix>0x80000000) z = (invsqrtpi*cc)/sqrtf(y);
if(ix>0x58000000) z = (invsqrtpi*cc)/sqrtf(y); /* |x|>2**49 */
else {
u = ponef(y); v = qonef(y);
z = invsqrtpi*(u*cc-v*ss)/sqrtf(y);
@ -71,7 +77,7 @@ __ieee754_j1f(float x)
if(hx<0) return -z;
else return z;
}
if(ix<0x32000000) { /* |x|<2**-27 */
if(ix<0x39000000) { /* |x|<2**-13 */
if(huge+x>one) return (float)0.5*x;/* inexact if x!=0 necessary */
}
z = x*x;
@ -104,10 +110,9 @@ __ieee754_y1f(float x)
GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx;
/* if Y1(NaN) is NaN, Y1(-inf) is NaN, Y1(inf) is 0 */
if(ix>=0x7f800000) return one/(x+x*x);
if(ix==0) return -one/zero;
if(hx<0) return zero/zero;
if(ix>=0x7f800000) return vone/(x+x*x);
if(ix==0) return -one/vzero;
if(hx<0) return vzero/vzero;
if(ix >= 0x40000000) { /* |x| >= 2.0 */
s = sinf(x);
c = cosf(x);
@ -129,14 +134,14 @@ __ieee754_y1f(float x)
* sin(x) +- cos(x) = -cos(2x)/(sin(x) -+ cos(x))
* to compute the worse one.
*/
if(ix>0x48000000) z = (invsqrtpi*ss)/sqrtf(x);
if(ix>0x58000000) z = (invsqrtpi*ss)/sqrtf(x); /* |x|>2**49 */
else {
u = ponef(x); v = qonef(x);
z = invsqrtpi*(u*ss+v*cc)/sqrtf(x);
}
return z;
}
if(ix<=0x24800000) { /* x < 2**-54 */
if(ix<=0x33000000) { /* x < 2**-25 */
return(-tpi/x);
}
z = x*x;
@ -219,7 +224,8 @@ static const float ps2[5] = {
8.3646392822e+00, /* 0x4105d590 */
};
static float ponef(float x)
static __inline float
ponef(float x)
{
const float *p,*q;
float z,r,s;
@ -227,9 +233,9 @@ static const float ps2[5] = {
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x41000000) {p = pr8; q= ps8;}
else if(ix>=0x40f71c58){p = pr5; q= ps5;}
else if(ix>=0x4036db68){p = pr3; q= ps3;}
else if(ix>=0x40000000){p = pr2; q= ps2;}
else if(ix>=0x409173eb){p = pr5; q= ps5;}
else if(ix>=0x4036d917){p = pr3; q= ps3;}
else {p = pr2; q= ps2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*q[4]))));
@ -315,17 +321,18 @@ static const float qs2[6] = {
-4.9594988823e+00, /* 0xc09eb437 */
};
static float qonef(float x)
static __inline float
qonef(float x)
{
const float *p,*q;
float s,r,z;
int32_t ix;
GET_FLOAT_WORD(ix,x);
ix &= 0x7fffffff;
if(ix>=0x40200000) {p = qr8; q= qs8;}
else if(ix>=0x40f71c58){p = qr5; q= qs5;}
else if(ix>=0x4036db68){p = qr3; q= qs3;}
else if(ix>=0x40000000){p = qr2; q= qs2;}
if(ix>=0x41000000) {p = qr8; q= qs8;}
else if(ix>=0x409173eb){p = qr5; q= qs5;}
else if(ix>=0x4036d917){p = qr3; q= qs3;}
else {p = qr2; q= qs2;} /* ix>=0x40000000 */
z = one/(x*x);
r = p[0]+z*(p[1]+z*(p[2]+z*(p[3]+z*(p[4]+z*p[5]))));
s = one+z*(q[0]+z*(q[1]+z*(q[2]+z*(q[3]+z*(q[4]+z*q[5])))));

12
libm/upstream-freebsd/lib/msun/src/e_jn.c
View File

@ -12,7 +12,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/e_jn.c 279856 2015-03-10 17:10:54Z kargl $");
/*
* __ieee754_jn(n, x), __ieee754_yn(n, x)
@ -43,6 +43,8 @@ __FBSDID("$FreeBSD$");
#include "math.h"
#include "math_private.h"
static const volatile double vone = 1, vzero = 0;
static const double
invsqrtpi= 5.64189583547756279280e-01, /* 0x3FE20DD7, 0x50429B6D */
two = 2.00000000000000000000e+00, /* 0x40000000, 0x00000000 */
@ -220,10 +222,12 @@ __ieee754_yn(int n, double x)
EXTRACT_WORDS(hx,lx,x);
ix = 0x7fffffff&hx;
/* if Y(n,NaN) is NaN */
/* yn(n,NaN) = NaN */
if((ix|((u_int32_t)(lx|-lx))>>31)>0x7ff00000) return x+x;
if((ix|lx)==0) return -one/zero;
if(hx<0) return zero/zero;
/* yn(n,+-0) = -inf and raise divide-by-zero exception. */
if((ix|lx)==0) return -one/vzero;
/* yn(n,x<0) = NaN and raise invalid exception. */
if(hx<0) return vzero/vzero;
sign = 1;
if(n<0){
n = -n;

13
libm/upstream-freebsd/lib/msun/src/e_jnf.c
View File

@ -14,11 +14,17 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/e_jnf.c 279856 2015-03-10 17:10:54Z kargl $");
/*
* See e_jn.c for complete comments.
*/
#include "math.h"
#include "math_private.h"
static const volatile float vone = 1, vzero = 0;
static const float
two = 2.0000000000e+00, /* 0x40000000 */
one = 1.0000000000e+00; /* 0x3F800000 */
@ -172,10 +178,9 @@ __ieee754_ynf(int n, float x)
GET_FLOAT_WORD(hx,x);
ix = 0x7fffffff&hx;
/* if Y(n,NaN) is NaN */
if(ix>0x7f800000) return x+x;
if(ix==0) return -one/zero;
if(hx<0) return zero/zero;
if(ix==0) return -one/vzero;
if(hx<0) return vzero/vzero;
sign = 1;
if(n<0){
n = -n;

4
libm/upstream-freebsd/lib/msun/src/k_exp.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/k_exp.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
@ -103,6 +103,6 @@ __ldexp_cexp(double complex z, int expt)
half_expt = expt - half_expt;
INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
return (cpack(cos(y) * exp_x * scale1 * scale2,
return (CMPLX(cos(y) * exp_x * scale1 * scale2,
sin(y) * exp_x * scale1 * scale2));
}

4
libm/upstream-freebsd/lib/msun/src/k_expf.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/k_expf.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
@ -82,6 +82,6 @@ __ldexp_cexpf(float complex z, int expt)
half_expt = expt - half_expt;
SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23);
return (cpackf(cosf(y) * exp_x * scale1 * scale2,
return (CMPLXF(cosf(y) * exp_x * scale1 * scale2,
sinf(y) * exp_x * scale1 * scale2));
}

20
libm/upstream-freebsd/lib/msun/src/math_private.h
View File

@ -11,7 +11,7 @@
/*
* from: @(#)fdlibm.h 5.1 93/09/24
* $FreeBSD$
* $FreeBSD: head/lib/msun/src/math_private.h 276176 2014-12-24 10:13:53Z ed $
*/
#ifndef _MATH_PRIVATE_H_
@ -454,9 +454,15 @@ typedef union {
* (0.0+I)*(y+0.0*I) and laboriously computing the full complex product.
* In particular, I*Inf is corrupted to NaN+I*Inf, and I*-0 is corrupted
* to -0.0+I*0.0.
*
* The C11 standard introduced the macros CMPLX(), CMPLXF() and CMPLXL()
* to construct complex values. Compilers that conform to the C99
* standard require the following functions to avoid the above issues.
*/
#ifndef CMPLXF
static __inline float complex
cpackf(float x, float y)
CMPLXF(float x, float y)
{
float_complex z;
@ -464,9 +470,11 @@ cpackf(float x, float y)
IMAGPART(z) = y;
return (z.f);
}
#endif
#ifndef CMPLX
static __inline double complex
cpack(double x, double y)
CMPLX(double x, double y)
{
double_complex z;
@ -474,9 +482,11 @@ cpack(double x, double y)
IMAGPART(z) = y;
return (z.f);
}
#endif
#ifndef CMPLXL
static __inline long double complex
cpackl(long double x, long double y)
CMPLXL(long double x, long double y)
{
long_double_complex z;
@ -484,6 +494,8 @@ cpackl(long double x, long double y)
IMAGPART(z) = y;
return (z.f);
}
#endif
#endif /* _COMPLEX_H */
#ifdef __GNUCLIKE_ASM

61
libm/upstream-freebsd/lib/msun/src/s_ccosh.c
View File

@ -32,10 +32,12 @@
*
* Exceptional values are noted in the comments within the source code.
* These values and the return value were taken from n1124.pdf.
* The sign of the result for some exceptional values is unspecified but
* must satisfy both cosh(conj(z)) == conj(cosh(z)) and cosh(-z) == cosh(z).
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_ccosh.c 284423 2015-06-15 20:11:06Z tijl $");
#include <complex.h>
#include <math.h>
@ -62,49 +64,48 @@ ccosh(double complex z)
/* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0)
return (cpack(cosh(x), x * y));
if (ix < 0x40360000) /* small x: normal case */
return (cpack(cosh(x) * cos(y), sinh(x) * sin(y)));
return (CMPLX(cosh(x), x * y));
if (ix < 0x40360000) /* |x| < 22: normal case */
return (CMPLX(cosh(x) * cos(y), sinh(x) * sin(y)));
/* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5;
return (cpack(h * cos(y), copysign(h, x) * sin(y)));
return (CMPLX(h * cos(y), copysign(h, x) * sin(y)));
} else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */
z = __ldexp_cexp(cpack(fabs(x), y), -1);
return (cpack(creal(z), cimag(z) * copysign(1, x)));
z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
return (CMPLX(creal(z), cimag(z) * copysign(1, x)));
} else {
/* x >= 1455: the result always overflows */
h = huge * x;
return (cpack(h * h * cos(y), h * sin(y)));
return (CMPLX(h * h * cos(y), h * sin(y)));
}
}
/*
* cosh(+-0 +- I Inf) = dNaN + I sign(d(+-0, dNaN))0.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as dNaN. Raise the invalid floating-point exception.
* cosh(+-0 +- I Inf) = dNaN + I (+-)(+-)0.
* The sign of 0 in the result is unspecified. Choice = product
* of the signs of the argument. Raise the invalid floating-point
* exception.
*
* cosh(+-0 +- I NaN) = d(NaN) + I sign(d(+-0, NaN))0.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as d(NaN).
* cosh(+-0 +- I NaN) = d(NaN) + I (+-)(+-)0.
* The sign of 0 in the result is unspecified. Choice = product
* of the signs of the argument.
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
return (cpack(y - y, copysign(0, x * (y - y))));
if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */
return (CMPLX(y - y, x * copysign(0, y)));
/*
* cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
*
* cosh(NaN +- I 0) = d(NaN) + I sign(d(NaN, +-0))0.
* The sign of 0 in the result is unspecified.
* cosh(NaN +- I 0) = d(NaN) + I (+-)(+-)0.
* The sign of 0 in the result is unspecified. Choice = product
* of the signs of the argument.
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
return (cpack(x * x, copysign(0, x) * y));
return (cpack(x * x, copysign(0, (x + x) * y)));
}
if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */
return (CMPLX(x * x, copysign(0, x) * y));
/*
* cosh(x +- I Inf) = dNaN + I dNaN.
@ -114,8 +115,8 @@ ccosh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
return (cpack(y - y, x * (y - y)));
if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */
return (CMPLX(y - y, x * (y - y)));
/*
* cosh(+-Inf + I NaN) = +Inf + I d(NaN).
@ -126,10 +127,10 @@ ccosh(double complex z)
*
* cosh(+-Inf + I y) = +Inf cos(y) +- I Inf sin(y)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (ix == 0x7ff00000 && lx == 0) {
if (iy >= 0x7ff00000)
return (cpack(x * x, x * (y - y)));
return (cpack((x * x) * cos(y), x * sin(y)));
return (CMPLX(INFINITY, x * (y - y)));
return (CMPLX(INFINITY * cos(y), x * sin(y)));
}
/*
@ -143,7 +144,7 @@ ccosh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
return (cpack((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLX((x * x) * (y - y), (x + x) * (y - y)));
}
double complex
@ -151,5 +152,5 @@ ccos(double complex z)
{
/* ccos(z) = ccosh(I * z) */
return (ccosh(cpack(-cimag(z), creal(z))));
return (ccosh(CMPLX(-cimag(z), creal(z))));
}

45
libm/upstream-freebsd/lib/msun/src/s_ccoshf.c
View File

@ -25,11 +25,11 @@
*/
/*
* Hyperbolic cosine of a complex argument. See s_ccosh.c for details.
* Float version of ccosh(). See s_ccosh.c for details.
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_ccoshf.c 284423 2015-06-15 20:11:06Z tijl $");
#include <complex.h>
#include <math.h>
@ -55,50 +55,47 @@ ccoshf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
return (cpackf(coshf(x), x * y));
if (ix < 0x41100000) /* small x: normal case */
return (cpackf(coshf(x) * cosf(y), sinhf(x) * sinf(y)));
return (CMPLXF(coshf(x), x * y));
if (ix < 0x41100000) /* |x| < 9: normal case */
return (CMPLXF(coshf(x) * cosf(y), sinhf(x) * sinf(y)));
/* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f;
return (cpackf(h * cosf(y), copysignf(h, x) * sinf(y)));
h = expf(fabsf(x)) * 0.5F;
return (CMPLXF(h * cosf(y), copysignf(h, x) * sinf(y)));
} else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */
z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
return (cpackf(crealf(z), cimagf(z) * copysignf(1, x)));
z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
return (CMPLXF(crealf(z), cimagf(z) * copysignf(1, x)));
} else {
/* x >= 192.7: the result always overflows */
h = huge * x;
return (cpackf(h * h * cosf(y), h * sinf(y)));
return (CMPLXF(h * h * cosf(y), h * sinf(y)));
}
}
if (ix == 0 && iy >= 0x7f800000)
return (cpackf(y - y, copysignf(0, x * (y - y))));
if (ix == 0) /* && iy >= 0x7f800000 */
return (CMPLXF(y - y, x * copysignf(0, y)));
if (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0)
return (cpackf(x * x, copysignf(0, x) * y));
return (cpackf(x * x, copysignf(0, (x + x) * y)));
}
if (iy == 0) /* && ix >= 0x7f800000 */
return (CMPLXF(x * x, copysignf(0, x) * y));
if (ix < 0x7f800000 && iy >= 0x7f800000)
return (cpackf(y - y, x * (y - y)));
if (ix < 0x7f800000) /* && iy >= 0x7f800000 */
return (CMPLXF(y - y, x * (y - y)));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (ix == 0x7f800000) {
if (iy >= 0x7f800000)
return (cpackf(x * x, x * (y - y)));
return (cpackf((x * x) * cosf(y), x * sinf(y)));
return (CMPLXF(INFINITY, x * (y - y)));
return (CMPLXF(INFINITY * cosf(y), x * sinf(y)));
}
return (cpackf((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLXF((x * x) * (y - y), (x + x) * (y - y)));
}
float complex
ccosf(float complex z)
{
return (ccoshf(cpackf(-cimagf(z), crealf(z))));
return (ccoshf(CMPLXF(-cimagf(z), crealf(z))));
}

14
libm/upstream-freebsd/lib/msun/src/s_cexp.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_cexp.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <math.h>
@ -50,22 +50,22 @@ cexp(double complex z)
/* cexp(x + I 0) = exp(x) + I 0 */
if ((hy | ly) == 0)
return (cpack(exp(x), y));
return (CMPLX(exp(x), y));
EXTRACT_WORDS(hx, lx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if (((hx & 0x7fffffff) | lx) == 0)
return (cpack(cos(y), sin(y)));
return (CMPLX(cos(y), sin(y)));
if (hy >= 0x7ff00000) {
if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
return (cpack(y - y, y - y));
return (CMPLX(y - y, y - y));
} else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
return (cpack(0.0, 0.0));
return (CMPLX(0.0, 0.0));
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
return (cpack(x, y - y));
return (CMPLX(x, y - y));
}
}
@ -84,6 +84,6 @@ cexp(double complex z)
* - x = NaN (spurious inexact exception from y)
*/
exp_x = exp(x);
return (cpack(exp_x * cos(y), exp_x * sin(y)));
return (CMPLX(exp_x * cos(y), exp_x * sin(y)));
}
}

14
libm/upstream-freebsd/lib/msun/src/s_cexpf.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_cexpf.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <math.h>
@ -50,22 +50,22 @@ cexpf(float complex z)
/* cexp(x + I 0) = exp(x) + I 0 */
if (hy == 0)
return (cpackf(expf(x), y));
return (CMPLXF(expf(x), y));
GET_FLOAT_WORD(hx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */
if ((hx & 0x7fffffff) == 0)
return (cpackf(cosf(y), sinf(y)));
return (CMPLXF(cosf(y), sinf(y)));
if (hy >= 0x7f800000) {
if ((hx & 0x7fffffff) != 0x7f800000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */
return (cpackf(y - y, y - y));
return (CMPLXF(y - y, y - y));
} else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
return (cpackf(0.0, 0.0));
return (CMPLXF(0.0, 0.0));
} else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
return (cpackf(x, y - y));
return (CMPLXF(x, y - y));
}
}
@ -84,6 +84,6 @@ cexpf(float complex z)
* - x = NaN (spurious inexact exception from y)
*/
exp_x = expf(x);
return (cpackf(exp_x * cosf(y), exp_x * sinf(y)));
return (CMPLXF(exp_x * cosf(y), exp_x * sinf(y)));
}
}

4
libm/upstream-freebsd/lib/msun/src/s_conj.c
View File

@ -23,7 +23,7 @@
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* $FreeBSD$
* $FreeBSD: head/lib/msun/src/s_conj.c 275819 2014-12-16 09:21:56Z ed $
*/
#include <complex.h>
@ -34,5 +34,5 @@ double complex
conj(double complex z)
{
return (cpack(creal(z), -cimag(z)));
return (CMPLX(creal(z), -cimag(z)));
}

4
libm/upstream-freebsd/lib/msun/src/s_conjf.c
View File

@ -23,7 +23,7 @@
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* $FreeBSD$
* $FreeBSD: head/lib/msun/src/s_conjf.c 275819 2014-12-16 09:21:56Z ed $
*/
#include <complex.h>
@ -34,5 +34,5 @@ float complex
conjf(float complex z)
{
return (cpackf(crealf(z), -cimagf(z)));
return (CMPLXF(crealf(z), -cimagf(z)));
}

4
libm/upstream-freebsd/lib/msun/src/s_conjl.c
View File

@ -23,7 +23,7 @@
* OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
* SUCH DAMAGE.
*
* $FreeBSD$
* $FreeBSD: head/lib/msun/src/s_conjl.c 275819 2014-12-16 09:21:56Z ed $
*/
#include <complex.h>
@ -34,5 +34,5 @@ long double complex
conjl(long double complex z)
{
return (cpackl(creall(z), -cimagl(z)));
return (CMPLXL(creall(z), -cimagl(z)));
}

4
libm/upstream-freebsd/lib/msun/src/s_cproj.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_cproj.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <math.h>
@ -39,7 +39,7 @@ cproj(double complex z)
if (!isinf(creal(z)) && !isinf(cimag(z)))
return (z);
else
return (cpack(INFINITY, copysign(0.0, cimag(z))));
return (CMPLX(INFINITY, copysign(0.0, cimag(z))));
}
#if LDBL_MANT_DIG == 53

4
libm/upstream-freebsd/lib/msun/src/s_cprojf.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_cprojf.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <math.h>
@ -39,5 +39,5 @@ cprojf(float complex z)
if (!isinf(crealf(z)) && !isinf(cimagf(z)))
return (z);
else
return (cpackf(INFINITY, copysignf(0.0, cimagf(z))));
return (CMPLXF(INFINITY, copysignf(0.0, cimagf(z))));
}

4
libm/upstream-freebsd/lib/msun/src/s_cprojl.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_cprojl.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <math.h>
@ -39,5 +39,5 @@ cprojl(long double complex z)
if (!isinf(creall(z)) && !isinf(cimagl(z)))
return (z);
else
return (cpackl(INFINITY, copysignl(0.0, cimagl(z))));
return (CMPLXL(INFINITY, copysignl(0.0, cimagl(z))));
}

75
libm/upstream-freebsd/lib/msun/src/s_csinh.c
View File

@ -32,10 +32,12 @@
*
* Exceptional values are noted in the comments within the source code.
* These values and the return value were taken from n1124.pdf.
* The sign of the result for some exceptional values is unspecified but
* must satisfy both sinh(conj(z)) == conj(sinh(z)) and sinh(-z) == -sinh(z).
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_csinh.c 284426 2015-06-15 20:16:53Z tijl $");
#include <complex.h>
#include <math.h>
@ -62,48 +64,45 @@ csinh(double complex z)
/* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0)
return (cpack(sinh(x), y));
if (ix < 0x40360000) /* small x: normal case */
return (cpack(sinh(x) * cos(y), cosh(x) * sin(y)));
return (CMPLX(sinh(x), y));
if (ix < 0x40360000) /* |x| < 22: normal case */
return (CMPLX(sinh(x) * cos(y), cosh(x) * sin(y)));
/* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5;
return (cpack(copysign(h, x) * cos(y), h * sin(y)));
return (CMPLX(copysign(h, x) * cos(y), h * sin(y)));
} else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */
z = __ldexp_cexp(cpack(fabs(x), y), -1);
return (cpack(creal(z) * copysign(1, x), cimag(z)));
z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
return (CMPLX(creal(z) * copysign(1, x), cimag(z)));
} else {
/* x >= 1455: the result always overflows */
h = huge * x;
return (cpack(h * cos(y), h * h * sin(y)));
return (CMPLX(h * cos(y), h * h * sin(y)));
}
}
/*
* sinh(+-0 +- I Inf) = sign(d(+-0, dNaN))0 + I dNaN.
* The sign of 0 in the result is unspecified. Choice = normally
* the same as dNaN. Raise the invalid floating-point exception.
* sinh(+-0 +- I Inf) = +-0 + I dNaN.
* The sign of 0 in the result is unspecified. Choice = same sign
* as the argument. Raise the invalid floating-point exception.
*
* sinh(+-0 +- I NaN) = sign(d(+-0, NaN))0 + I d(NaN).
* The sign of 0 in the result is unspecified. Choice = normally
* the same as d(NaN).
* sinh(+-0 +- I NaN) = +-0 + I d(NaN).
* The sign of 0 in the result is unspecified. Choice = same sign
* as the argument.
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
return (cpack(copysign(0, x * (y - y)), y - y));
if ((ix | lx) == 0) /* && iy >= 0x7ff00000 */
return (CMPLX(x, y - y));
/*
* sinh(+-Inf +- I 0) = +-Inf + I +-0.
*
* sinh(NaN +- I 0) = d(NaN) + I +-0.
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
return (cpack(x, y));
return (cpack(x, copysign(0, y)));
}
if ((iy | ly) == 0) /* && ix >= 0x7ff00000 */
return (CMPLX(x + x, y));
/*
* sinh(x +- I Inf) = dNaN + I dNaN.
@ -113,45 +112,45 @@ csinh(double complex z)
* Optionally raises the invalid floating-point exception for finite
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
return (cpack(y - y, x * (y - y)));
if (ix < 0x7ff00000) /* && iy >= 0x7ff00000 */
return (CMPLX(y - y, y - y));
/*
* sinh(+-Inf + I NaN) = +-Inf + I d(NaN).
* The sign of Inf in the result is unspecified. Choice = normally
* the same as d(NaN).
* The sign of Inf in the result is unspecified. Choice = same sign
* as the argument.
*
* sinh(+-Inf +- I Inf) = +Inf + I dNaN.
* The sign of Inf in the result is unspecified. Choice = always +.
* Raise the invalid floating-point exception.
* sinh(+-Inf +- I Inf) = +-Inf + I dNaN.
* The sign of Inf in the result is unspecified. Choice = same sign
* as the argument. Raise the invalid floating-point exception.
*
* sinh(+-Inf + I y) = +-Inf cos(y) + I Inf sin(y)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (ix == 0x7ff00000 && lx == 0) {
if (iy >= 0x7ff00000)
return (cpack(x * x, x * (y - y)));
return (cpack(x * cos(y), INFINITY * sin(y)));
return (CMPLX(x, y - y));
return (CMPLX(x * cos(y), INFINITY * sin(y)));
}
/*
* sinh(NaN + I NaN) = d(NaN) + I d(NaN).
* sinh(NaN1 + I NaN2) = d(NaN1, NaN2) + I d(NaN1, NaN2).
*
* sinh(NaN +- I Inf) = d(NaN) + I d(NaN).
* sinh(NaN +- I Inf) = d(NaN, dNaN) + I d(NaN, dNaN).
* Optionally raises the invalid floating-point exception.
* Choice = raise.
*
* sinh(NaN + I y) = d(NaN) + I d(NaN).
* sinh(NaN + I y) = d(NaN) + I d(NaN).
* Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs).
*/
return (cpack((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLX((x + x) * (y - y), (x * x) * (y - y)));
}
double complex
csin(double complex z)
{
/* csin(z) = -I * csinh(I * z) */
z = csinh(cpack(-cimag(z), creal(z)));
return (cpack(cimag(z), -creal(z)));
/* csin(z) = -I * csinh(I * z) = I * conj(csinh(I * conj(z))). */
z = csinh(CMPLX(cimag(z), creal(z)));
return (CMPLX(cimag(z), creal(z)));
}

47
libm/upstream-freebsd/lib/msun/src/s_csinhf.c
View File

@ -25,11 +25,11 @@
*/
/*
* Hyperbolic sine of a complex argument z. See s_csinh.c for details.
* Float version of csinh(). See s_csinh.c for details.
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_csinhf.c 284426 2015-06-15 20:16:53Z tijl $");
#include <complex.h>
#include <math.h>
@ -55,51 +55,48 @@ csinhf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
return (cpackf(sinhf(x), y));
if (ix < 0x41100000) /* small x: normal case */
return (cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y)));
return (CMPLXF(sinhf(x), y));
if (ix < 0x41100000) /* |x| < 9: normal case */
return (CMPLXF(sinhf(x) * cosf(y), coshf(x) * sinf(y)));
/* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f;
return (cpackf(copysignf(h, x) * cosf(y), h * sinf(y)));
h = expf(fabsf(x)) * 0.5F;
return (CMPLXF(copysignf(h, x) * cosf(y), h * sinf(y)));
} else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */
z = __ldexp_cexpf(cpackf(fabsf(x), y), -1);
return (cpackf(crealf(z) * copysignf(1, x), cimagf(z)));
z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
return (CMPLXF(crealf(z) * copysignf(1, x), cimagf(z)));
} else {
/* x >= 192.7: the result always overflows */
h = huge * x;
return (cpackf(h * cosf(y), h * h * sinf(y)));
return (CMPLXF(h * cosf(y), h * h * sinf(y)));
}
}
if (ix == 0 && iy >= 0x7f800000)
return (cpackf(copysignf(0, x * (y - y)), y - y));
if (ix == 0) /* && iy >= 0x7f800000 */
return (CMPLXF(x, y - y));
if (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0)
return (cpackf(x, y));
return (cpackf(x, copysignf(0, y)));
}
if (iy == 0) /* && ix >= 0x7f800000 */
return (CMPLXF(x + x, y));
if (ix < 0x7f800000 && iy >= 0x7f800000)
return (cpackf(y - y, x * (y - y)));
if (ix < 0x7f800000) /* && iy >= 0x7f800000 */
return (CMPLXF(y - y, y - y));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (ix == 0x7f800000) {
if (iy >= 0x7f800000)
return (cpackf(x * x, x * (y - y)));
return (cpackf(x * cosf(y), INFINITY * sinf(y)));
return (CMPLXF(x, y - y));
return (CMPLXF(x * cosf(y), INFINITY * sinf(y)));
}
return (cpackf((x * x) * (y - y), (x + x) * (y - y)));
return (CMPLXF((x + x) * (y - y), (x * x) * (y - y)));
}
float complex
csinf(float complex z)
{
z = csinhf(cpackf(-cimagf(z), crealf(z)));
return (cpackf(cimagf(z), -crealf(z)));
z = csinhf(CMPLXF(cimagf(z), crealf(z)));
return (CMPLXF(cimagf(z), crealf(z)));
}

16
libm/upstream-freebsd/lib/msun/src/s_csqrt.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrt.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <float.h>
@ -58,12 +58,12 @@ csqrt(double complex z)
/* Handle special cases. */
if (z == 0)
return (cpack(0, b));
return (CMPLX(0, b));
if (isinf(b))
return (cpack(INFINITY, b));
return (CMPLX(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (cpack(a, t)); /* return NaN + NaN i */
return (CMPLX(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@ -73,9 +73,9 @@ csqrt(double complex z)
* csqrt(-inf + y i) = 0 + inf i
*/
if (signbit(a))
return (cpack(fabs(b - b), copysign(a, b)));
return (CMPLX(fabs(b - b), copysign(a, b)));
else
return (cpack(a, copysign(b - b, b)));
return (CMPLX(a, copysign(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
@ -94,10 +94,10 @@ csqrt(double complex z)
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
result = cpack(t, b / (2 * t));
result = CMPLX(t, b / (2 * t));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
result = cpack(fabs(b) / (2 * t), copysign(t, b));
result = CMPLX(fabs(b) / (2 * t), copysign(t, b));
}
/* Rescale. */

16
libm/upstream-freebsd/lib/msun/src/s_csqrtf.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrtf.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <math.h>
@ -49,12 +49,12 @@ csqrtf(float complex z)
/* Handle special cases. */
if (z == 0)
return (cpackf(0, b));
return (CMPLXF(0, b));
if (isinf(b))
return (cpackf(INFINITY, b));
return (CMPLXF(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (cpackf(a, t)); /* return NaN + NaN i */
return (CMPLXF(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@ -64,9 +64,9 @@ csqrtf(float complex z)
* csqrtf(-inf + y i) = 0 + inf i
*/
if (signbit(a))
return (cpackf(fabsf(b - b), copysignf(a, b)));
return (CMPLXF(fabsf(b - b), copysignf(a, b)));
else
return (cpackf(a, copysignf(b - b, b)));
return (CMPLXF(a, copysignf(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
@ -80,9 +80,9 @@ csqrtf(float complex z)
*/
if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5);
return (cpackf(t, b / (2.0 * t)));
return (CMPLXF(t, b / (2.0 * t)));
} else {
t = sqrt((-a + hypot(a, b)) * 0.5);
return (cpackf(fabsf(b) / (2.0 * t), copysignf(t, b)));
return (CMPLXF(fabsf(b) / (2.0 * t), copysignf(t, b)));
}
}

16
libm/upstream-freebsd/lib/msun/src/s_csqrtl.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_csqrtl.c 275819 2014-12-16 09:21:56Z ed $");
#include <complex.h>
#include <float.h>
@ -58,12 +58,12 @@ csqrtl(long double complex z)
/* Handle special cases. */
if (z == 0)
return (cpackl(0, b));
return (CMPLXL(0, b));
if (isinf(b))
return (cpackl(INFINITY, b));
return (CMPLXL(INFINITY, b));
if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */
return (cpackl(a, t)); /* return NaN + NaN i */
return (CMPLXL(a, t)); /* return NaN + NaN i */
}
if (isinf(a)) {
/*
@ -73,9 +73,9 @@ csqrtl(long double complex z)
* csqrt(-inf + y i) = 0 + inf i
*/
if (signbit(a))
return (cpackl(fabsl(b - b), copysignl(a, b)));
return (CMPLXL(fabsl(b - b), copysignl(a, b)));
else
return (cpackl(a, copysignl(b - b, b)));
return (CMPLXL(a, copysignl(b - b, b)));
}
/*
* The remaining special case (b is NaN) is handled just fine by
@ -94,10 +94,10 @@ csqrtl(long double complex z)
/* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) {
t = sqrtl((a + hypotl(a, b)) * 0.5);
result = cpackl(t, b / (2 * t));
result = CMPLXL(t, b / (2 * t));
} else {
t = sqrtl((-a + hypotl(a, b)) * 0.5);
result = cpackl(fabsl(b) / (2 * t), copysignl(t, b));
result = CMPLXL(fabsl(b) / (2 * t), copysignl(t, b));
}
/* Rescale. */

45
libm/upstream-freebsd/lib/msun/src/s_ctanh.c
View File

@ -25,7 +25,7 @@
*/
/*
* Hyperbolic tangent of a complex argument z = x + i y.
* Hyperbolic tangent of a complex argument z = x + I y.
*
* The algorithm is from:
*
@ -44,15 +44,15 @@
*
* tanh(z) = sinh(z) / cosh(z)
*
* sinh(x) cos(y) + i cosh(x) sin(y)
* sinh(x) cos(y) + I cosh(x) sin(y)
* = ---------------------------------
* cosh(x) cos(y) + i sinh(x) sin(y)
* cosh(x) cos(y) + I sinh(x) sin(y)
*
* cosh(x) sinh(x) / cos^2(y) + i tan(y)
* cosh(x) sinh(x) / cos^2(y) + I tan(y)
* = -------------------------------------
* 1 + sinh^2(x) / cos^2(y)
*
* beta rho s + i t
* beta rho s + I t
* = ----------------
* 1 + beta s^2
*
@ -64,7 +64,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_ctanh.c 284427 2015-06-15 20:40:44Z tijl $");
#include <complex.h>
#include <math.h>
@ -85,16 +85,16 @@ ctanh(double complex z)
ix = hx & 0x7fffffff;
/*
* ctanh(NaN + i 0) = NaN + i 0
* ctanh(NaN +- I 0) = d(NaN) +- I 0
*
* ctanh(NaN + i y) = NaN + i NaN for y != 0
* ctanh(NaN + I y) = d(NaN,y) + I d(NaN,y) for y != 0
*
* The imaginary part has the sign of x*sin(2*y), but there's no
* special effort to get this right.
*
* ctanh(+-Inf +- i Inf) = +-1 +- 0
* ctanh(+-Inf +- I Inf) = +-1 +- I 0
*
* ctanh(+-Inf + i y) = +-1 + 0 sin(2y) for y finite
* ctanh(+-Inf + I y) = +-1 + I 0 sin(2y) for y finite
*
* The imaginary part of the sign is unspecified. This special
* case is only needed to avoid a spurious invalid exception when
@ -102,26 +102,27 @@ ctanh(double complex z)
*/
if (ix >= 0x7ff00000) {
if ((ix & 0xfffff) | lx) /* x is NaN */
return (cpack(x, (y == 0 ? y : x * y)));
return (CMPLX((x + 0) * (y + 0),
y == 0 ? y : (x + 0) * (y + 0)));
SET_HIGH_WORD(x, hx - 0x40000000); /* x = copysign(1, x) */
return (cpack(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))));
return (CMPLX(x, copysign(0, isinf(y) ? y : sin(y) * cos(y))));
}
/*
* ctanh(x + i NAN) = NaN + i NaN
* ctanh(x +- i Inf) = NaN + i NaN
* ctanh(x + I NaN) = d(NaN) + I d(NaN)
* ctanh(x +- I Inf) = dNaN + I dNaN
*/
if (!isfinite(y))
return (cpack(y - y, y - y));
return (CMPLX(y - y, y - y));
/*
* ctanh(+-huge + i +-y) ~= +-1 +- i 2sin(2y)/exp(2x), using the
* ctanh(+-huge +- I y) ~= +-1 +- I 2sin(2y)/exp(2x), using the
* approximation sinh^2(huge) ~= exp(2*huge) / 4.
* We use a modified formula to avoid spurious overflow.
*/
if (ix >= 0x40360000) { /* x >= 22 */
if (ix >= 0x40360000) { /* |x| >= 22 */
double exp_mx = exp(-fabs(x));
return (cpack(copysign(1, x),
return (CMPLX(copysign(1, x),
4 * sin(y) * cos(y) * exp_mx * exp_mx));
}
@ -131,14 +132,14 @@ ctanh(double complex z)
s = sinh(x);
rho = sqrt(1 + s * s); /* = cosh(x) */
denom = 1 + beta * s * s;
return (cpack((beta * rho * s) / denom, t / denom));
return (CMPLX((beta * rho * s) / denom, t / denom));
}
double complex
ctan(double complex z)
{
/* ctan(z) = -I * ctanh(I * z) */
z = ctanh(cpack(-cimag(z), creal(z)));
return (cpack(cimag(z), -creal(z)));
/* ctan(z) = -I * ctanh(I * z) = I * conj(ctanh(I * conj(z))) */
z = ctanh(CMPLX(cimag(z), creal(z)));
return (CMPLX(cimag(z), creal(z)));
}

19
libm/upstream-freebsd/lib/msun/src/s_ctanhf.c
View File

@ -29,7 +29,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_ctanhf.c 284428 2015-06-15 20:47:26Z tijl $");
#include <complex.h>
#include <math.h>
@ -51,18 +51,19 @@ ctanhf(float complex z)
if (ix >= 0x7f800000) {
if (ix & 0x7fffff)
return (cpackf(x, (y == 0 ? y : x * y)));
return (CMPLXF((x + 0) * (y + 0),
y == 0 ? y : (x + 0) * (y + 0)));
SET_FLOAT_WORD(x, hx - 0x40000000);
return (cpackf(x,
return (CMPLXF(x,
copysignf(0, isinf(y) ? y : sinf(y) * cosf(y))));
}
if (!isfinite(y))
return (cpackf(y - y, y - y));
return (CMPLXF(y - y, y - y));
if (ix >= 0x41300000) { /* x >= 11 */
if (ix >= 0x41300000) { /* |x| >= 11 */
float exp_mx = expf(-fabsf(x));
return (cpackf(copysignf(1, x),
return (CMPLXF(copysignf(1, x),
4 * sinf(y) * cosf(y) * exp_mx * exp_mx));
}
@ -71,14 +72,14 @@ ctanhf(float complex z)
s = sinhf(x);
rho = sqrtf(1 + s * s);
denom = 1 + beta * s * s;
return (cpackf((beta * rho * s) / denom, t / denom));
return (CMPLXF((beta * rho * s) / denom, t / denom));
}
float complex
ctanf(float complex z)
{
z = ctanhf(cpackf(-cimagf(z), crealf(z)));
return (cpackf(cimagf(z), -crealf(z)));
z = ctanhf(CMPLXF(cimagf(z), crealf(z)));
return (CMPLXF(cimagf(z), crealf(z)));
}

6
libm/upstream-freebsd/lib/msun/src/s_exp2.c
View File

@ -25,7 +25,7 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_exp2.c 286515 2015-08-09 10:00:13Z dim $");
#include <float.h>
@ -376,14 +376,14 @@ exp2(double x)
/* Compute r = exp2(y) = exp2t[i0] * p(z - eps[i]). */
t = tbl[i0]; /* exp2t[i0] */
z -= tbl[i0 + 1]; /* eps[i0] */
if (k >= -1021 << 20)
if (k >= -(1021 << 20))
INSERT_WORDS(twopk, 0x3ff00000 + k, 0);
else
INSERT_WORDS(twopkp1000, 0x3ff00000 + k + (1000 << 20), 0);
r = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));
/* Scale by 2**(k>>20). */
if(k >= -1021 << 20) {
if(k >= -(1021 << 20)) {
if (k == 1024 << 20)
return (r * 2.0 * 0x1p1023);
return (r * twopk);

46
libm/upstream-freebsd/lib/msun/src/s_scalbln.c
View File

@ -25,52 +25,30 @@
*/
#include <sys/cdefs.h>
__FBSDID("$FreeBSD$");
__FBSDID("$FreeBSD: head/lib/msun/src/s_scalbln.c 278339 2015-02-07 00:38:18Z kargl $");
#include <limits.h>
#include <math.h>
double
scalbln (double x, long n)
{
int in;
#define NMAX 65536
#define NMIN -65536
in = (int)n;
if (in != n) {
if (n > 0)
in = INT_MAX;
else
in = INT_MIN;
}
return (scalbn(x, in));
double
scalbln(double x, long n)
{
return (scalbn(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
}
float
scalblnf (float x, long n)
scalblnf(float x, long n)
{
int in;
in = (int)n;
if (in != n) {
if (n > 0)
in = INT_MAX;
else
in = INT_MIN;
}
return (scalbnf(x, in));
return (scalbnf(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
}
long double
scalblnl (long double x, long n)
scalblnl(long double x, long n)
{
int in;
in = (int)n;
if (in != n) {
if (n > 0)
in = INT_MAX;
else
in = INT_MIN;
}
return (scalbnl(x, (int)n));
return (scalbnl(x, (n > NMAX) ? NMAX : (n < NMIN) ? NMIN : (int)n));
}