From 5d6cb09b963ec564ad6244d357294ac1444cfc34 Mon Sep 17 00:00:00 2001 From: Ed Schouten Date: Thu, 8 Jan 2015 12:02:56 +0100 Subject: [PATCH] Change existing code to use CMPLX*() instead of cpack*() where possible. I am currently working on building openlibm against stock copies of , instead of the openlibm.h header. It seems that a C compliant header can be used as a drop-in replacement for openlibm.h, with the exception that it lacks cpack*(). In FreeBSD SVN r275819 I patched up the math library by replacing cpack*() by CMPLX*(). That way many functions become less tied to the intrinsics of the math library. Make the same change to openlibm. --- src/k_exp.c | 2 +- src/k_expf.c | 2 +- src/openlibm.h | 44 ++++++++++++++++++++++++++++++-------------- src/s_ccosh.c | 28 ++++++++++++++-------------- src/s_ccoshf.c | 28 ++++++++++++++-------------- src/s_cexp.c | 12 ++++++------ src/s_cexpf.c | 12 ++++++------ src/s_conj.c | 2 +- src/s_conjf.c | 2 +- src/s_conjl.c | 2 +- src/s_cproj.c | 2 +- src/s_cprojf.c | 2 +- src/s_cprojl.c | 2 +- src/s_csinh.c | 30 +++++++++++++++--------------- src/s_csinhf.c | 30 +++++++++++++++--------------- src/s_csqrt.c | 14 +++++++------- src/s_csqrtf.c | 14 +++++++------- src/s_csqrtl.c | 14 +++++++------- src/s_ctanh.c | 14 +++++++------- src/s_ctanhf.c | 14 +++++++------- 20 files changed, 143 insertions(+), 127 deletions(-) diff --git a/src/k_exp.c b/src/k_exp.c index 8735088..7bed88b 100644 --- a/src/k_exp.c +++ b/src/k_exp.c @@ -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)); } diff --git a/src/k_expf.c b/src/k_expf.c index c94da5c..9df0b0f 100644 --- a/src/k_expf.c +++ b/src/k_expf.c @@ -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)); } diff --git a/src/openlibm.h b/src/openlibm.h index 5635612..991dff6 100644 --- a/src/openlibm.h +++ b/src/openlibm.h @@ -204,7 +204,7 @@ typedef union { #define IMAGPART(z) ((z).a[1]) /* - * Inline functions that can be used to construct complex values. + * Macros that can be used to construct complex values. * * The C99 standard intends x+I*y to be used for this, but x+I*y is * currently unusable in general since gcc introduces many overflow, @@ -217,18 +217,20 @@ typedef union { * and gcc 4.7 added a __builtin_complex feature to simplify implementation * of CMPLX in libc, so we can take advantage of these features if they * are available. + * + * If __builtin_complex is not available, resort to using inline + * functions instead. These can unfortunately not be used to construct + * compile-time constants. */ -#if defined(CMPLXF) && defined(CMPLX) && defined(CMPLXL) /* C11 */ -# define cpackf(x,y) CMPLXF(x,y) -# define cpack(x,y) CMPLX(x,y) -# define cpackl(x,y) CMPLXL(x,y) -#elif (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && !defined(__INTEL_COMPILER) -# define cpackf(x,y) __builtin_complex ((float) (x), (float) (y)) -# define cpack(x,y) __builtin_complex ((double) (x), (double) (y)) -# define cpackl(x,y) __builtin_complex ((long double) (x), (long double) (y)) -#else /* define our own cpack functions */ + +#define HAVE_BUILTIN_COMPLEX (__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ >= 7)) && !defined(__INTEL_COMPILER) + +#ifndef CMPLXF +#if HAVE_BUILTIN_COMPLEX +# define CMPLXF(x,y) __builtin_complex ((float) (x), (float) (y)) +#else static __inline float complex -cpackf(float x, float y) +CMPLXF(float x, float y) { float_complex z; @@ -236,9 +238,15 @@ cpackf(float x, float y) IMAGPART(z) = y; return (z.f); } +#endif +#endif +#ifndef CMPLX +#if HAVE_BUILTIN_COMPLEX +# define CMPLX(x,y) __builtin_complex ((double) (x), (double) (y)) +#else static __inline double complex -cpack(double x, double y) +CMPLX(double x, double y) { double_complex z; @@ -246,9 +254,15 @@ cpack(double x, double y) IMAGPART(z) = y; return (z.f); } +#endif +#endif +#ifndef CMPLXL +#if HAVE_BUILTIN_COMPLEX +# define CMPLXL(x,y) __builtin_complex ((long double) (x), (long double) (y)) +#else static __inline long double complex -cpackl(long double x, long double y) +CMPLXL(long double x, long double y) { long_double_complex z; @@ -256,7 +270,9 @@ cpackl(long double x, long double y) IMAGPART(z) = y; return (z.f); } -#endif /* define our own cpack functions */ +#endif +#endif + //VBS //#endif /* _COMPLEX_H */ diff --git a/src/s_ccosh.c b/src/s_ccosh.c index 8f55f42..32b780a 100644 --- a/src/s_ccosh.c +++ b/src/s_ccosh.c @@ -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))); } DLLEXPORT 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)))); } diff --git a/src/s_ccoshf.c b/src/s_ccoshf.c index 6fcf836..6d79ac6 100644 --- a/src/s_ccoshf.c +++ b/src/s_ccoshf.c @@ -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))); } DLLEXPORT float complex ccosf(float complex z) { - return (ccoshf(cpackf(-cimagf(z), crealf(z)))); + return (ccoshf(CMPLXF(-cimagf(z), crealf(z)))); } diff --git a/src/s_cexp.c b/src/s_cexp.c index 38e6653..f247cfe 100644 --- a/src/s_cexp.c +++ b/src/s_cexp.c @@ -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))); } } diff --git a/src/s_cexpf.c b/src/s_cexpf.c index 90eab17..26e963a 100644 --- a/src/s_cexpf.c +++ b/src/s_cexpf.c @@ -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))); } } diff --git a/src/s_conj.c b/src/s_conj.c index 4dbf1bc..b383f69 100644 --- a/src/s_conj.c +++ b/src/s_conj.c @@ -35,5 +35,5 @@ DLLEXPORT double complex conj(double complex z) { - return (cpack(creal(z), -cimag(z))); + return (CMPLX(creal(z), -cimag(z))); } diff --git a/src/s_conjf.c b/src/s_conjf.c index 9b1a37a..2b24486 100644 --- a/src/s_conjf.c +++ b/src/s_conjf.c @@ -35,5 +35,5 @@ DLLEXPORT float complex conjf(float complex z) { - return (cpackf(crealf(z), -cimagf(z))); + return (CMPLXF(crealf(z), -cimagf(z))); } diff --git a/src/s_conjl.c b/src/s_conjl.c index 477ccc1..4b73472 100644 --- a/src/s_conjl.c +++ b/src/s_conjl.c @@ -35,5 +35,5 @@ DLLEXPORT long double complex conjl(long double complex z) { - return (cpackl(creall(z), -cimagl(z))); + return (CMPLXL(creall(z), -cimagl(z))); } diff --git a/src/s_cproj.c b/src/s_cproj.c index 67f0bd0..303812c 100644 --- a/src/s_cproj.c +++ b/src/s_cproj.c @@ -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 diff --git a/src/s_cprojf.c b/src/s_cprojf.c index 4750b6b..30c7f47 100644 --- a/src/s_cprojf.c +++ b/src/s_cprojf.c @@ -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)))); } diff --git a/src/s_cprojl.c b/src/s_cprojl.c index a048107..a10553b 100644 --- a/src/s_cprojl.c +++ b/src/s_cprojl.c @@ -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)))); } diff --git a/src/s_csinh.c b/src/s_csinh.c index 8402c49..9f2fd95 100644 --- a/src/s_csinh.c +++ b/src/s_csinh.c @@ -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))); } DLLEXPORT 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))); } diff --git a/src/s_csinhf.c b/src/s_csinhf.c index 78ed180..e5545ab 100644 --- a/src/s_csinhf.c +++ b/src/s_csinhf.c @@ -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))); } DLLEXPORT 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))); } diff --git a/src/s_csqrt.c b/src/s_csqrt.c index efd0701..7ddf944 100644 --- a/src/s_csqrt.c +++ b/src/s_csqrt.c @@ -60,12 +60,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)) { /* @@ -75,9 +75,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 @@ -96,10 +96,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. */ diff --git a/src/s_csqrtf.c b/src/s_csqrtf.c index 7a79c88..58c4b05 100644 --- a/src/s_csqrtf.c +++ b/src/s_csqrtf.c @@ -51,12 +51,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)) { /* @@ -66,9 +66,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 @@ -82,9 +82,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))); } } diff --git a/src/s_csqrtl.c b/src/s_csqrtl.c index c542a13..416e6b6 100644 --- a/src/s_csqrtl.c +++ b/src/s_csqrtl.c @@ -59,12 +59,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)) { /* @@ -74,9 +74,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 @@ -95,10 +95,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. */ diff --git a/src/s_ctanh.c b/src/s_ctanh.c index d26d128..0042bbe 100644 --- a/src/s_ctanh.c +++ b/src/s_ctanh.c @@ -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)); } DLLEXPORT 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))); } diff --git a/src/s_ctanhf.c b/src/s_ctanhf.c index c123b54..cc70c5e 100644 --- a/src/s_ctanhf.c +++ b/src/s_ctanhf.c @@ -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)); } DLLEXPORT 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))); }