Merge pull request #76 from NuxiNL/cmplx

Change existing code to use CMPLX*() instead of cpack*() where possible.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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