Rename cpack*() to CMPLX*().

The C11 standard introduced a set of macros (CMPLX, CMPLXF, CMPLXL) that
can be used to construct complex numbers from a pair of real and
imaginary numbers. Unfortunately, they require some compiler support,
which is why we only define them for Clang and GCC>=4.7.

The cpack() function in libm performs the same task as CMPLX(), but
cannot be used to generate compile-time constants. This means that all
invocations of cpack() can safely be replaced by C11's CMPLX(). To keep
the code building with GCC 4.2, provide copies of CMPLX() that can at
least be used to generate run-time complex numbers.

This makes it easier to build some of the functions outside of libm.
This commit is contained in:
Ed Schouten 2014-12-16 09:21:56 +00:00
parent c7d73a4d23
commit 2cec876a59
24 changed files with 195 additions and 182 deletions

View File

@ -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 */

View File

@ -299,7 +299,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 */

View File

@ -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)));
}

View File

@ -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)));
}

View File

@ -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));
}

View File

@ -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));
}

View File

@ -454,9 +454,16 @@ 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. The functions below are modelled after
* these macros, with the exception that they cannot be used to
* construct compile-time complex values.
*/
#ifndef CMPLXF
static __inline float complex
cpackf(float x, float y)
CMPLXF(float x, float y)
{
float_complex z;
@ -464,9 +471,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 +483,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 +495,8 @@ cpackl(long double x, long double y)
IMAGPART(z) = y;
return (z.f);
}
#endif
#endif /* _COMPLEX_H */
#ifdef __GNUCLIKE_ASM

View File

@ -62,23 +62,23 @@ 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));
return (CMPLX(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) * 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)));
}
}
@ -92,7 +92,7 @@ ccosh(double complex z)
* the same as d(NaN).
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
return (cpack(y - y, copysign(0, x * (y - y))));
return (CMPLX(y - y, copysign(0, x * (y - y))));
/*
* cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
@ -102,8 +102,8 @@ ccosh(double complex z)
*/
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)));
return (CMPLX(x * x, copysign(0, x) * y));
return (CMPLX(x * x, copysign(0, (x + x) * y)));
}
/*
@ -115,7 +115,7 @@ ccosh(double complex z)
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
return (cpack(y - y, x * (y - y)));
return (CMPLX(y - y, x * (y - y)));
/*
* cosh(+-Inf + I NaN) = +Inf + I d(NaN).
@ -128,8 +128,8 @@ ccosh(double complex z)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (iy >= 0x7ff00000)
return (cpack(x * x, x * (y - y)));
return (cpack((x * x) * cos(y), x * sin(y)));
return (CMPLX(x * x, x * (y - y)));
return (CMPLX((x * x) * cos(y), x * sin(y)));
}
/*
@ -143,7 +143,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 +151,5 @@ ccos(double complex z)
{
/* ccos(z) = ccosh(I * z) */
return (ccosh(cpack(-cimag(z), creal(z))));
return (ccosh(CMPLX(-cimag(z), creal(z))));
}

View File

@ -55,50 +55,50 @@ ccoshf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
return (cpackf(coshf(x), x * y));
return (CMPLXF(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) * 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)));
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))));
return (CMPLXF(y - y, copysignf(0, x * (y - 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)));
return (CMPLXF(x * x, copysignf(0, x) * y));
return (CMPLXF(x * x, copysignf(0, (x + x) * y)));
}
if (ix < 0x7f800000 && iy >= 0x7f800000)
return (cpackf(y - y, x * (y - y)));
return (CMPLXF(y - y, x * (y - y)));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (iy >= 0x7f800000)
return (cpackf(x * x, x * (y - y)));
return (cpackf((x * x) * cosf(y), x * sinf(y)));
return (CMPLXF(x * x, x * (y - y)));
return (CMPLXF((x * x) * 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))));
}

View File

@ -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)));
}
}

View File

@ -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)));
}
}

View File

@ -34,5 +34,5 @@ double complex
conj(double complex z)
{
return (cpack(creal(z), -cimag(z)));
return (CMPLX(creal(z), -cimag(z)));
}

View File

@ -34,5 +34,5 @@ float complex
conjf(float complex z)
{
return (cpackf(crealf(z), -cimagf(z)));
return (CMPLXF(crealf(z), -cimagf(z)));
}

View File

@ -34,5 +34,5 @@ long double complex
conjl(long double complex z)
{
return (cpackl(creall(z), -cimagl(z)));
return (CMPLXL(creall(z), -cimagl(z)));
}

View File

@ -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

View File

@ -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))));
}

View File

@ -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))));
}

View File

@ -62,23 +62,23 @@ 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));
return (CMPLX(sinh(x), y));
if (ix < 0x40360000) /* small x: normal case */
return (cpack(sinh(x) * cos(y), cosh(x) * sin(y)));
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)));
}
}
@ -92,7 +92,7 @@ csinh(double complex z)
* the same as d(NaN).
*/
if ((ix | lx) == 0 && iy >= 0x7ff00000)
return (cpack(copysign(0, x * (y - y)), y - y));
return (CMPLX(copysign(0, x * (y - y)), y - y));
/*
* sinh(+-Inf +- I 0) = +-Inf + I +-0.
@ -101,8 +101,8 @@ csinh(double complex z)
*/
if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0)
return (cpack(x, y));
return (cpack(x, copysign(0, y)));
return (CMPLX(x, y));
return (CMPLX(x, copysign(0, y)));
}
/*
@ -114,7 +114,7 @@ csinh(double complex z)
* nonzero x. Choice = don't raise (except for signaling NaNs).
*/
if (ix < 0x7ff00000 && iy >= 0x7ff00000)
return (cpack(y - y, x * (y - y)));
return (CMPLX(y - y, x * (y - y)));
/*
* sinh(+-Inf + I NaN) = +-Inf + I d(NaN).
@ -129,8 +129,8 @@ csinh(double complex z)
*/
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (iy >= 0x7ff00000)
return (cpack(x * x, x * (y - y)));
return (cpack(x * cos(y), INFINITY * sin(y)));
return (CMPLX(x * x, x * (y - y)));
return (CMPLX(x * cos(y), INFINITY * sin(y)));
}
/*
@ -144,7 +144,7 @@ csinh(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
@ -152,6 +152,6 @@ csin(double complex z)
{
/* csin(z) = -I * csinh(I * z) */
z = csinh(cpack(-cimag(z), creal(z)));
return (cpack(cimag(z), -creal(z)));
z = csinh(CMPLX(-cimag(z), creal(z)));
return (CMPLX(cimag(z), -creal(z)));
}

View File

@ -55,51 +55,51 @@ csinhf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0)
return (cpackf(sinhf(x), y));
return (CMPLXF(sinhf(x), y));
if (ix < 0x41100000) /* small x: normal case */
return (cpackf(sinhf(x) * cosf(y), coshf(x) * sinf(y)));
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)));
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));
return (CMPLXF(copysignf(0, x * (y - y)), y - y));
if (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0)
return (cpackf(x, y));
return (cpackf(x, copysignf(0, y)));
return (CMPLXF(x, y));
return (CMPLXF(x, copysignf(0, y)));
}
if (ix < 0x7f800000 && iy >= 0x7f800000)
return (cpackf(y - y, x * (y - y)));
return (CMPLXF(y - y, x * (y - y)));
if (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (iy >= 0x7f800000)
return (cpackf(x * x, x * (y - y)));
return (cpackf(x * cosf(y), INFINITY * sinf(y)));
return (CMPLXF(x * x, 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)));
}

View File

@ -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. */

View File

@ -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)));
}
}

View File

@ -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. */

View File

@ -102,9 +102,9 @@ 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, (y == 0 ? y : x * y)));
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))));
}
/*
@ -112,7 +112,7 @@ ctanh(double complex z)
* ctanh(x +- i Inf) = NaN + i NaN
*/
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
@ -121,7 +121,7 @@ ctanh(double complex z)
*/
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,7 +131,7 @@ 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
@ -139,6 +139,6 @@ ctan(double complex z)
{
/* ctan(z) = -I * ctanh(I * z) */
z = ctanh(cpack(-cimag(z), creal(z)));
return (cpack(cimag(z), -creal(z)));
z = ctanh(CMPLX(-cimag(z), creal(z)));
return (CMPLX(cimag(z), -creal(z)));
}

View File

@ -51,18 +51,18 @@ ctanhf(float complex z)
if (ix >= 0x7f800000) {
if (ix & 0x7fffff)
return (cpackf(x, (y == 0 ? y : x * y)));
return (CMPLXF(x, (y == 0 ? y : x * y)));
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 */
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 +71,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)));
}