complex: add C11 CMPLX macros and replace cpack with them

This commit is contained in:
Szabolcs Nagy
2012-11-13 01:31:49 +01:00
parent e2fe959fe2
commit cfbaba79a2
47 changed files with 152 additions and 156 deletions

View File

@ -112,6 +112,15 @@ long double creall(long double complex);
#define cimagf(x) __CIMAG(x, float) #define cimagf(x) __CIMAG(x, float)
#define cimagl(x) __CIMAG(x, long double) #define cimagl(x) __CIMAG(x, long double)
#define __CMPLX(x, y, t) \
((union { _Complex t __z; t __xy[2]; }){.__xy = {(x),(y)}}.__z)
#if __STDC_VERSION__ >= 201112L
#define CMPLX(x, y) __CMPLX(x, y, double)
#define CMPLXF(x, y) __CMPLX(x, y, float)
#define CMPLXL(x, y) __CMPLX(x, y, long double)
#endif
#ifdef __cplusplus #ifdef __cplusplus
} }
#endif #endif

View File

@ -83,5 +83,5 @@ double complex __ldexp_cexp(double complex z, int expt)
half_expt = expt - half_expt; half_expt = expt - half_expt;
INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0); INSERT_WORDS(scale2, (0x3ff + half_expt) << 20, 0);
return cpack(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2); return CMPLX(cos(y) * exp_x * scale1 * scale2, sin(y) * exp_x * scale1 * scale2);
} }

View File

@ -63,6 +63,6 @@ float complex __ldexp_cexpf(float complex z, int expt)
half_expt = expt - half_expt; half_expt = expt - half_expt;
SET_FLOAT_WORD(scale2, (0x7f + half_expt) << 23); 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); sinf(y) * exp_x * scale1 * scale2);
} }

View File

@ -7,5 +7,5 @@
double complex cacos(double complex z) double complex cacos(double complex z)
{ {
z = casin(z); z = casin(z);
return cpack(M_PI_2 - creal(z), -cimag(z)); return CMPLX(M_PI_2 - creal(z), -cimag(z));
} }

View File

@ -5,5 +5,5 @@
float complex cacosf(float complex z) float complex cacosf(float complex z)
{ {
z = casinf(z); z = casinf(z);
return cpackf((float)M_PI_2 - crealf(z), -cimagf(z)); return CMPLXF((float)M_PI_2 - crealf(z), -cimagf(z));
} }

View File

@ -5,5 +5,5 @@
double complex cacosh(double complex z) double complex cacosh(double complex z)
{ {
z = cacos(z); z = cacos(z);
return cpack(-cimag(z), creal(z)); return CMPLX(-cimag(z), creal(z));
} }

View File

@ -3,5 +3,5 @@
float complex cacoshf(float complex z) float complex cacoshf(float complex z)
{ {
z = cacosf(z); z = cacosf(z);
return cpackf(-cimagf(z), crealf(z)); return CMPLXF(-cimagf(z), crealf(z));
} }

View File

@ -9,6 +9,6 @@ long double complex cacoshl(long double complex z)
long double complex cacoshl(long double complex z) long double complex cacoshl(long double complex z)
{ {
z = cacosl(z); z = cacosl(z);
return cpackl(-cimagl(z), creall(z)); return CMPLXL(-cimagl(z), creall(z));
} }
#endif #endif

View File

@ -11,6 +11,6 @@ long double complex cacosl(long double complex z)
long double complex cacosl(long double complex z) long double complex cacosl(long double complex z)
{ {
z = casinl(z); z = casinl(z);
return cpackl(PI_2 - creall(z), -cimagl(z)); return CMPLXL(PI_2 - creall(z), -cimagl(z));
} }
#endif #endif

View File

@ -11,6 +11,6 @@ double complex casin(double complex z)
x = creal(z); x = creal(z);
y = cimag(z); y = cimag(z);
w = cpack(1.0 - (x - y)*(x + y), -2.0*x*y); w = CMPLX(1.0 - (x - y)*(x + y), -2.0*x*y);
return clog(cpack(-y, x) + csqrt(w)); return clog(CMPLX(-y, x) + csqrt(w));
} }

View File

@ -9,6 +9,6 @@ float complex casinf(float complex z)
x = crealf(z); x = crealf(z);
y = cimagf(z); y = cimagf(z);
w = cpackf(1.0 - (x - y)*(x + y), -2.0*x*y); w = CMPLXF(1.0 - (x - y)*(x + y), -2.0*x*y);
return clogf(cpackf(-y, x) + csqrtf(w)); return clogf(CMPLXF(-y, x) + csqrtf(w));
} }

View File

@ -4,6 +4,6 @@
double complex casinh(double complex z) double complex casinh(double complex z)
{ {
z = casin(cpack(-cimag(z), creal(z))); z = casin(CMPLX(-cimag(z), creal(z)));
return cpack(cimag(z), -creal(z)); return CMPLX(cimag(z), -creal(z));
} }

View File

@ -2,6 +2,6 @@
float complex casinhf(float complex z) float complex casinhf(float complex z)
{ {
z = casinf(cpackf(-cimagf(z), crealf(z))); z = casinf(CMPLXF(-cimagf(z), crealf(z)));
return cpackf(cimagf(z), -crealf(z)); return CMPLXF(cimagf(z), -crealf(z));
} }

View File

@ -8,7 +8,7 @@ long double complex casinhl(long double complex z)
#else #else
long double complex casinhl(long double complex z) long double complex casinhl(long double complex z)
{ {
z = casinl(cpackl(-cimagl(z), creall(z))); z = casinl(CMPLXL(-cimagl(z), creall(z)));
return cpackl(cimagl(z), -creall(z)); return CMPLXL(cimagl(z), -creall(z));
} }
#endif #endif

View File

@ -14,7 +14,7 @@ long double complex casinl(long double complex z)
x = creall(z); x = creall(z);
y = cimagl(z); y = cimagl(z);
w = cpackl(1.0 - (x - y)*(x + y), -2.0*x*y); w = CMPLXL(1.0 - (x - y)*(x + y), -2.0*x*y);
return clogl(cpackl(-y, x) + csqrtl(w)); return clogl(CMPLXL(-y, x) + csqrtl(w));
} }
#endif #endif

View File

@ -4,6 +4,6 @@
double complex catanh(double complex z) double complex catanh(double complex z)
{ {
z = catan(cpack(-cimag(z), creal(z))); z = catan(CMPLX(-cimag(z), creal(z)));
return cpack(cimag(z), -creal(z)); return CMPLX(cimag(z), -creal(z));
} }

View File

@ -2,6 +2,6 @@
float complex catanhf(float complex z) float complex catanhf(float complex z)
{ {
z = catanf(cpackf(-cimagf(z), crealf(z))); z = catanf(CMPLXF(-cimagf(z), crealf(z)));
return cpackf(cimagf(z), -crealf(z)); return CMPLXF(cimagf(z), -crealf(z));
} }

View File

@ -8,7 +8,7 @@ long double complex catanhl(long double complex z)
#else #else
long double complex catanhl(long double complex z) long double complex catanhl(long double complex z)
{ {
z = catanl(cpackl(-cimagl(z), creall(z))); z = catanl(CMPLXL(-cimagl(z), creall(z)));
return cpackl(cimagl(z), -creall(z)); return CMPLXL(cimagl(z), -creall(z));
} }
#endif #endif

View File

@ -4,5 +4,5 @@
double complex ccos(double complex z) double complex ccos(double complex z)
{ {
return ccosh(cpack(-cimag(z), creal(z))); return ccosh(CMPLX(-cimag(z), creal(z)));
} }

View File

@ -2,5 +2,5 @@
float complex ccosf(float complex z) float complex ccosf(float complex z)
{ {
return ccoshf(cpackf(-cimagf(z), crealf(z))); return ccoshf(CMPLXF(-cimagf(z), crealf(z)));
} }

View File

@ -55,23 +55,23 @@ double complex ccosh(double complex z)
/* Handle the nearly-non-exceptional cases where x and y are finite. */ /* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) { if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0) if ((iy | ly) == 0)
return cpack(cosh(x), x * y); return CMPLX(cosh(x), x * y);
if (ix < 0x40360000) /* small x: normal case */ 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|) */ /* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) { if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */ /* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5; 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) { } else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */ /* x < 1455: scale to avoid overflow */
z = __ldexp_cexp(cpack(fabs(x), y), -1); z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
return cpack(creal(z), cimag(z) * copysign(1, x)); return CMPLX(creal(z), cimag(z) * copysign(1, x));
} else { } else {
/* x >= 1455: the result always overflows */ /* x >= 1455: the result always overflows */
h = huge * x; h = huge * x;
return cpack(h * h * cos(y), h * sin(y)); return CMPLX(h * h * cos(y), h * sin(y));
} }
} }
@ -85,7 +85,7 @@ double complex ccosh(double complex z)
* the same as d(NaN). * the same as d(NaN).
*/ */
if ((ix | lx) == 0 && iy >= 0x7ff00000) 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. * cosh(+-Inf +- I 0) = +Inf + I (+-)(+-)0.
@ -95,8 +95,8 @@ double complex ccosh(double complex z)
*/ */
if ((iy | ly) == 0 && ix >= 0x7ff00000) { if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0) if (((hx & 0xfffff) | lx) == 0)
return cpack(x * x, copysign(0, x) * y); return CMPLX(x * x, copysign(0, x) * y);
return cpack(x * x, copysign(0, (x + x) * y)); return CMPLX(x * x, copysign(0, (x + x) * y));
} }
/* /*
@ -108,7 +108,7 @@ double complex ccosh(double complex z)
* nonzero x. Choice = don't raise (except for signaling NaNs). * nonzero x. Choice = don't raise (except for signaling NaNs).
*/ */
if (ix < 0x7ff00000 && iy >= 0x7ff00000) 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). * cosh(+-Inf + I NaN) = +Inf + I d(NaN).
@ -121,8 +121,8 @@ double complex ccosh(double complex z)
*/ */
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (iy >= 0x7ff00000) if (iy >= 0x7ff00000)
return cpack(x * x, x * (y - y)); return CMPLX(x * x, x * (y - y));
return cpack((x * x) * cos(y), x * sin(y)); return CMPLX((x * x) * cos(y), x * sin(y));
} }
/* /*
@ -136,5 +136,5 @@ double complex ccosh(double complex z)
* Optionally raises the invalid floating-point exception for finite * Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs). * 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));
} }

View File

@ -48,43 +48,43 @@ float complex ccoshf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) { if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0) if (iy == 0)
return cpackf(coshf(x), x * y); return CMPLXF(coshf(x), x * y);
if (ix < 0x41100000) /* small x: normal case */ 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|) */ /* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) { if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */ /* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f; 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) { } else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */ /* x < 192.7: scale to avoid overflow */
z = __ldexp_cexpf(cpackf(fabsf(x), y), -1); z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
return cpackf(crealf(z), cimagf(z) * copysignf(1, x)); return CMPLXF(crealf(z), cimagf(z) * copysignf(1, x));
} else { } else {
/* x >= 192.7: the result always overflows */ /* x >= 192.7: the result always overflows */
h = huge * x; 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) 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 (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0) if ((hx & 0x7fffff) == 0)
return cpackf(x * x, copysignf(0, x) * y); return CMPLXF(x * x, copysignf(0, x) * y);
return cpackf(x * x, copysignf(0, (x + x) * y)); return CMPLXF(x * x, copysignf(0, (x + x) * y));
} }
if (ix < 0x7f800000 && iy >= 0x7f800000) 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 (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (iy >= 0x7f800000) if (iy >= 0x7f800000)
return cpackf(x * x, x * (y - y)); return CMPLXF(x * x, x * (y - y));
return cpackf((x * x) * cosf(y), x * sinf(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));
} }

View File

@ -8,6 +8,6 @@ long double complex ccosl(long double complex z)
#else #else
long double complex ccosl(long double complex z) long double complex ccosl(long double complex z)
{ {
return ccoshl(cpackl(-cimagl(z), creall(z))); return ccoshl(CMPLXL(-cimagl(z), creall(z)));
} }
#endif #endif

View File

@ -44,22 +44,22 @@ double complex cexp(double complex z)
/* cexp(x + I 0) = exp(x) + I 0 */ /* cexp(x + I 0) = exp(x) + I 0 */
if ((hy | ly) == 0) if ((hy | ly) == 0)
return cpack(exp(x), y); return CMPLX(exp(x), y);
EXTRACT_WORDS(hx, lx, x); EXTRACT_WORDS(hx, lx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */ /* cexp(0 + I y) = cos(y) + I sin(y) */
if (((hx & 0x7fffffff) | lx) == 0) if (((hx & 0x7fffffff) | lx) == 0)
return cpack(cos(y), sin(y)); return CMPLX(cos(y), sin(y));
if (hy >= 0x7ff00000) { if (hy >= 0x7ff00000) {
if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) { if (lx != 0 || (hx & 0x7fffffff) != 0x7ff00000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ /* 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) { } else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
return cpack(0.0, 0.0); return CMPLX(0.0, 0.0);
} else { } else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
return cpack(x, y - y); return CMPLX(x, y - y);
} }
} }
@ -78,6 +78,6 @@ double complex cexp(double complex z)
* - x = NaN (spurious inexact exception from y) * - x = NaN (spurious inexact exception from y)
*/ */
exp_x = exp(x); 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

@ -44,22 +44,22 @@ float complex cexpf(float complex z)
/* cexp(x + I 0) = exp(x) + I 0 */ /* cexp(x + I 0) = exp(x) + I 0 */
if (hy == 0) if (hy == 0)
return cpackf(expf(x), y); return CMPLXF(expf(x), y);
GET_FLOAT_WORD(hx, x); GET_FLOAT_WORD(hx, x);
/* cexp(0 + I y) = cos(y) + I sin(y) */ /* cexp(0 + I y) = cos(y) + I sin(y) */
if ((hx & 0x7fffffff) == 0) if ((hx & 0x7fffffff) == 0)
return cpackf(cosf(y), sinf(y)); return CMPLXF(cosf(y), sinf(y));
if (hy >= 0x7f800000) { if (hy >= 0x7f800000) {
if ((hx & 0x7fffffff) != 0x7f800000) { if ((hx & 0x7fffffff) != 0x7f800000) {
/* cexp(finite|NaN +- I Inf|NaN) = NaN + I NaN */ /* 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) { } else if (hx & 0x80000000) {
/* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */ /* cexp(-Inf +- I Inf|NaN) = 0 + I 0 */
return cpackf(0.0, 0.0); return CMPLXF(0.0, 0.0);
} else { } else {
/* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */ /* cexp(+Inf +- I Inf|NaN) = Inf + I NaN */
return cpackf(x, y - y); return CMPLXF(x, y - y);
} }
} }
@ -78,6 +78,6 @@ float complex cexpf(float complex z)
* - x = NaN (spurious inexact exception from y) * - x = NaN (spurious inexact exception from y)
*/ */
exp_x = expf(x); 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

@ -10,5 +10,5 @@ double complex clog(double complex z)
r = cabs(z); r = cabs(z);
phi = carg(z); phi = carg(z);
return cpack(log(r), phi); return CMPLX(log(r), phi);
} }

View File

@ -8,5 +8,5 @@ float complex clogf(float complex z)
r = cabsf(z); r = cabsf(z);
phi = cargf(z); phi = cargf(z);
return cpackf(logf(r), phi); return CMPLXF(logf(r), phi);
} }

View File

@ -13,6 +13,6 @@ long double complex clogl(long double complex z)
r = cabsl(z); r = cabsl(z);
phi = cargl(z); phi = cargl(z);
return cpackl(logl(r), phi); return CMPLXL(logl(r), phi);
} }
#endif #endif

View File

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

View File

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

View File

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

View File

@ -3,6 +3,6 @@
double complex cproj(double complex z) double complex cproj(double complex z)
{ {
if (isinf(creal(z)) || isinf(cimag(z))) if (isinf(creal(z)) || isinf(cimag(z)))
return cpack(INFINITY, copysign(0.0, creal(z))); return CMPLX(INFINITY, copysign(0.0, creal(z)));
return z; return z;
} }

View File

@ -3,6 +3,6 @@
float complex cprojf(float complex z) float complex cprojf(float complex z)
{ {
if (isinf(crealf(z)) || isinf(cimagf(z))) if (isinf(crealf(z)) || isinf(cimagf(z)))
return cpackf(INFINITY, copysignf(0.0, crealf(z))); return CMPLXF(INFINITY, copysignf(0.0, crealf(z)));
return z; return z;
} }

View File

@ -9,7 +9,7 @@ long double complex cprojl(long double complex z)
long double complex cprojl(long double complex z) long double complex cprojl(long double complex z)
{ {
if (isinf(creall(z)) || isinf(cimagl(z))) if (isinf(creall(z)) || isinf(cimagl(z)))
return cpackl(INFINITY, copysignl(0.0, creall(z))); return CMPLXL(INFINITY, copysignl(0.0, creall(z)));
return z; return z;
} }
#endif #endif

View File

@ -4,6 +4,6 @@
double complex csin(double complex z) double complex csin(double complex z)
{ {
z = csinh(cpack(-cimag(z), creal(z))); z = csinh(CMPLX(-cimag(z), creal(z)));
return cpack(cimag(z), -creal(z)); return CMPLX(cimag(z), -creal(z));
} }

View File

@ -2,6 +2,6 @@
float complex csinf(float complex z) float complex csinf(float complex z)
{ {
z = csinhf(cpackf(-cimagf(z), crealf(z))); z = csinhf(CMPLXF(-cimagf(z), crealf(z)));
return cpackf(cimagf(z), -crealf(z)); return CMPLXF(cimagf(z), -crealf(z));
} }

View File

@ -55,23 +55,23 @@ double complex csinh(double complex z)
/* Handle the nearly-non-exceptional cases where x and y are finite. */ /* Handle the nearly-non-exceptional cases where x and y are finite. */
if (ix < 0x7ff00000 && iy < 0x7ff00000) { if (ix < 0x7ff00000 && iy < 0x7ff00000) {
if ((iy | ly) == 0) if ((iy | ly) == 0)
return cpack(sinh(x), y); return CMPLX(sinh(x), y);
if (ix < 0x40360000) /* small x: normal case */ 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|) */ /* |x| >= 22, so cosh(x) ~= exp(|x|) */
if (ix < 0x40862e42) { if (ix < 0x40862e42) {
/* x < 710: exp(|x|) won't overflow */ /* x < 710: exp(|x|) won't overflow */
h = exp(fabs(x)) * 0.5; 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) { } else if (ix < 0x4096bbaa) {
/* x < 1455: scale to avoid overflow */ /* x < 1455: scale to avoid overflow */
z = __ldexp_cexp(cpack(fabs(x), y), -1); z = __ldexp_cexp(CMPLX(fabs(x), y), -1);
return cpack(creal(z) * copysign(1, x), cimag(z)); return CMPLX(creal(z) * copysign(1, x), cimag(z));
} else { } else {
/* x >= 1455: the result always overflows */ /* x >= 1455: the result always overflows */
h = huge * x; h = huge * x;
return cpack(h * cos(y), h * h * sin(y)); return CMPLX(h * cos(y), h * h * sin(y));
} }
} }
@ -85,7 +85,7 @@ double complex csinh(double complex z)
* the same as d(NaN). * the same as d(NaN).
*/ */
if ((ix | lx) == 0 && iy >= 0x7ff00000) 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. * sinh(+-Inf +- I 0) = +-Inf + I +-0.
@ -94,8 +94,8 @@ double complex csinh(double complex z)
*/ */
if ((iy | ly) == 0 && ix >= 0x7ff00000) { if ((iy | ly) == 0 && ix >= 0x7ff00000) {
if (((hx & 0xfffff) | lx) == 0) if (((hx & 0xfffff) | lx) == 0)
return cpack(x, y); return CMPLX(x, y);
return cpack(x, copysign(0, y)); return CMPLX(x, copysign(0, y));
} }
/* /*
@ -107,7 +107,7 @@ double complex csinh(double complex z)
* nonzero x. Choice = don't raise (except for signaling NaNs). * nonzero x. Choice = don't raise (except for signaling NaNs).
*/ */
if (ix < 0x7ff00000 && iy >= 0x7ff00000) 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). * sinh(+-Inf + I NaN) = +-Inf + I d(NaN).
@ -122,8 +122,8 @@ double complex csinh(double complex z)
*/ */
if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) { if (ix >= 0x7ff00000 && ((hx & 0xfffff) | lx) == 0) {
if (iy >= 0x7ff00000) if (iy >= 0x7ff00000)
return cpack(x * x, x * (y - y)); return CMPLX(x * x, x * (y - y));
return cpack(x * cos(y), INFINITY * sin(y)); return CMPLX(x * cos(y), INFINITY * sin(y));
} }
/* /*
@ -137,5 +137,5 @@ double complex csinh(double complex z)
* Optionally raises the invalid floating-point exception for finite * Optionally raises the invalid floating-point exception for finite
* nonzero y. Choice = don't raise (except for signaling NaNs). * 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));
} }

View File

@ -48,43 +48,43 @@ float complex csinhf(float complex z)
if (ix < 0x7f800000 && iy < 0x7f800000) { if (ix < 0x7f800000 && iy < 0x7f800000) {
if (iy == 0) if (iy == 0)
return cpackf(sinhf(x), y); return CMPLXF(sinhf(x), y);
if (ix < 0x41100000) /* small x: normal case */ 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|) */ /* |x| >= 9, so cosh(x) ~= exp(|x|) */
if (ix < 0x42b17218) { if (ix < 0x42b17218) {
/* x < 88.7: expf(|x|) won't overflow */ /* x < 88.7: expf(|x|) won't overflow */
h = expf(fabsf(x)) * 0.5f; 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) { } else if (ix < 0x4340b1e7) {
/* x < 192.7: scale to avoid overflow */ /* x < 192.7: scale to avoid overflow */
z = __ldexp_cexpf(cpackf(fabsf(x), y), -1); z = __ldexp_cexpf(CMPLXF(fabsf(x), y), -1);
return cpackf(crealf(z) * copysignf(1, x), cimagf(z)); return CMPLXF(crealf(z) * copysignf(1, x), cimagf(z));
} else { } else {
/* x >= 192.7: the result always overflows */ /* x >= 192.7: the result always overflows */
h = huge * x; 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) 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 (iy == 0 && ix >= 0x7f800000) {
if ((hx & 0x7fffff) == 0) if ((hx & 0x7fffff) == 0)
return cpackf(x, y); return CMPLXF(x, y);
return cpackf(x, copysignf(0, y)); return CMPLXF(x, copysignf(0, y));
} }
if (ix < 0x7f800000 && iy >= 0x7f800000) 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 (ix >= 0x7f800000 && (hx & 0x7fffff) == 0) {
if (iy >= 0x7f800000) if (iy >= 0x7f800000)
return cpackf(x * x, x * (y - y)); return CMPLXF(x * x, x * (y - y));
return cpackf(x * cosf(y), INFINITY * sinf(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));
} }

View File

@ -8,7 +8,7 @@ long double complex csinl(long double complex z)
#else #else
long double complex csinl(long double complex z) long double complex csinl(long double complex z)
{ {
z = csinhl(cpackl(-cimagl(z), creall(z))); z = csinhl(CMPLXL(-cimagl(z), creall(z)));
return cpackl(cimagl(z), -creall(z)); return CMPLXL(cimagl(z), -creall(z));
} }
#endif #endif

View File

@ -51,12 +51,12 @@ double complex csqrt(double complex z)
/* Handle special cases. */ /* Handle special cases. */
if (z == 0) if (z == 0)
return cpack(0, b); return CMPLX(0, b);
if (isinf(b)) if (isinf(b))
return cpack(INFINITY, b); return CMPLX(INFINITY, b);
if (isnan(a)) { if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ 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)) { if (isinf(a)) {
/* /*
@ -66,9 +66,9 @@ double complex csqrt(double complex z)
* csqrt(-inf + y i) = 0 + inf i * csqrt(-inf + y i) = 0 + inf i
*/ */
if (signbit(a)) if (signbit(a))
return cpack(fabs(b - b), copysign(a, b)); return CMPLX(fabs(b - b), copysign(a, b));
else 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 * The remaining special case (b is NaN) is handled just fine by
@ -87,10 +87,10 @@ double complex csqrt(double complex z)
/* Algorithm 312, CACM vol 10, Oct 1967. */ /* Algorithm 312, CACM vol 10, Oct 1967. */
if (a >= 0) { if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5); t = sqrt((a + hypot(a, b)) * 0.5);
result = cpack(t, b / (2 * t)); result = CMPLX(t, b / (2 * t));
} else { } else {
t = sqrt((-a + hypot(a, b)) * 0.5); 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. */ /* Rescale. */

View File

@ -43,12 +43,12 @@ float complex csqrtf(float complex z)
/* Handle special cases. */ /* Handle special cases. */
if (z == 0) if (z == 0)
return cpackf(0, b); return CMPLXF(0, b);
if (isinf(b)) if (isinf(b))
return cpackf(INFINITY, b); return CMPLXF(INFINITY, b);
if (isnan(a)) { if (isnan(a)) {
t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ 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)) { if (isinf(a)) {
/* /*
@ -58,9 +58,9 @@ float complex csqrtf(float complex z)
* csqrtf(-inf + y i) = 0 + inf i * csqrtf(-inf + y i) = 0 + inf i
*/ */
if (signbit(a)) if (signbit(a))
return cpackf(fabsf(b - b), copysignf(a, b)); return CMPLXF(fabsf(b - b), copysignf(a, b));
else 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 * The remaining special case (b is NaN) is handled just fine by
@ -74,9 +74,9 @@ float complex csqrtf(float complex z)
*/ */
if (a >= 0) { if (a >= 0) {
t = sqrt((a + hypot(a, b)) * 0.5); t = sqrt((a + hypot(a, b)) * 0.5);
return cpackf(t, b / (2.0 * t)); return CMPLXF(t, b / (2.0 * t));
} else { } else {
t = sqrt((-a + hypot(a, b)) * 0.5); 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

@ -4,6 +4,6 @@
double complex ctan(double complex z) double complex ctan(double complex z)
{ {
z = ctanh(cpack(-cimag(z), creal(z))); z = ctanh(CMPLX(-cimag(z), creal(z)));
return cpack(cimag(z), -creal(z)); return CMPLX(cimag(z), -creal(z));
} }

View File

@ -2,6 +2,6 @@
float complex ctanf(float complex z) float complex ctanf(float complex z)
{ {
z = ctanhf(cpackf(-cimagf(z), crealf(z))); z = ctanhf(CMPLXF(-cimagf(z), crealf(z)));
return cpackf(cimagf(z), -crealf(z)); return CMPLXF(cimagf(z), -crealf(z));
} }

View File

@ -95,9 +95,9 @@ double complex ctanh(double complex z)
*/ */
if (ix >= 0x7ff00000) { if (ix >= 0x7ff00000) {
if ((ix & 0xfffff) | lx) /* x is NaN */ 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) */ 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)));
} }
/* /*
@ -105,7 +105,7 @@ double complex ctanh(double complex z)
* ctanh(x +- i Inf) = NaN + i NaN * ctanh(x +- i Inf) = NaN + i NaN
*/ */
if (!isfinite(y)) 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
@ -114,7 +114,7 @@ double complex ctanh(double complex z)
*/ */
if (ix >= 0x40360000) { /* x >= 22 */ if (ix >= 0x40360000) { /* x >= 22 */
double exp_mx = exp(-fabs(x)); double exp_mx = exp(-fabs(x));
return cpack(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx); return CMPLX(copysign(1, x), 4 * sin(y) * cos(y) * exp_mx * exp_mx);
} }
/* Kahan's algorithm */ /* Kahan's algorithm */
@ -123,5 +123,5 @@ double complex ctanh(double complex z)
s = sinh(x); s = sinh(x);
rho = sqrt(1 + s * s); /* = cosh(x) */ rho = sqrt(1 + s * s); /* = cosh(x) */
denom = 1 + beta * s * s; denom = 1 + beta * s * s;
return cpack((beta * rho * s) / denom, t / denom); return CMPLX((beta * rho * s) / denom, t / denom);
} }

View File

@ -44,17 +44,17 @@ float complex ctanhf(float complex z)
if (ix >= 0x7f800000) { if (ix >= 0x7f800000) {
if (ix & 0x7fffff) 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); SET_FLOAT_WORD(x, hx - 0x40000000);
return cpackf(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y))); return CMPLXF(x, copysignf(0, isinf(y) ? y : sinf(y) * cosf(y)));
} }
if (!isfinite(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)); float exp_mx = expf(-fabsf(x));
return cpackf(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx); return CMPLXF(copysignf(1, x), 4 * sinf(y) * cosf(y) * exp_mx * exp_mx);
} }
t = tanf(y); t = tanf(y);
@ -62,5 +62,5 @@ float complex ctanhf(float complex z)
s = sinhf(x); s = sinhf(x);
rho = sqrtf(1 + s * s); rho = sqrtf(1 + s * s);
denom = 1 + beta * s * s; denom = 1 + beta * s * s;
return cpackf((beta * rho * s) / denom, t / denom); return CMPLXF((beta * rho * s) / denom, t / denom);
} }

View File

@ -8,7 +8,7 @@ long double complex ctanl(long double complex z)
#else #else
long double complex ctanl(long double complex z) long double complex ctanl(long double complex z)
{ {
z = ctanhl(cpackl(-cimagl(z), creall(z))); z = ctanhl(CMPLXL(-cimagl(z), creall(z)));
return cpackl(cimagl(z), -creall(z)); return CMPLXL(cimagl(z), -creall(z));
} }
#endif #endif

View File

@ -170,25 +170,12 @@ long double __p1evll(long double, const long double *, int);
#define STRICT_ASSIGN(type, lval, rval) ((lval) = (type)(rval)) #define STRICT_ASSIGN(type, lval, rval) ((lval) = (type)(rval))
#endif #endif
/* complex */ /* complex */
union dcomplex { #ifndef CMPLX
double complex z; #define CMPLX(x, y) __CMPLX(x, y, double)
double a[2]; #define CMPLXF(x, y) __CMPLX(x, y, float)
}; #define CMPLXL(x, y) __CMPLX(x, y, long double)
union fcomplex { #endif
float complex z;
float a[2];
};
union lcomplex {
long double complex z;
long double a[2];
};
/* x + y*I is not supported properly by gcc */
#define cpack(x,y) ((union dcomplex){.a={(x),(y)}}.z)
#define cpackf(x,y) ((union fcomplex){.a={(x),(y)}}.z)
#define cpackl(x,y) ((union lcomplex){.a={(x),(y)}}.z)
#endif #endif