- Correct the range check in the double version to catch negative values

  that would overflow.
- Style fixes and improved handling of NaNs suggested by bde.

2 files changed, 22 insertions, 19 deletions

diff --git a/lib/msun/src/s_csqrt.c b/lib/msun/src/s_csqrt.c
index 3b1e187..50e24ba 100644
--- a/lib/msun/src/s_csqrt.c
+++ b/lib/msun/src/s_csqrt.c
@@ -41,17 +41,20 @@ __FBSDID("$FreeBSD$");
  */
 #pragma	STDC CX_LIMITED_RANGE	on
 
-/* We risk spurious overflow for components >= DBL_MAX/(1+sqrt(2)) */
-#define THRESH	0x1.a827999fcef32p+1022
+/* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */
+#define	THRESH	0x1.a827999fcef32p+1022
 
 double complex
 csqrt(double complex z)
 {
-	double a = creal(z), b = cimag(z);
-	double t;
 	double complex result;
+	double a, b;
+	double t;
 	int scale;
 
+	a = creal(z);
+	b = cimag(z);
+
 	/* Handle special cases. */
 	if (z == 0)
 		return (cpack(0, b));
@@ -59,13 +62,13 @@ csqrt(double complex z)
 		return (cpack(INFINITY, b));
 	if (isnan(a)) {
 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
-		return (cpack(t, t));	/* return NaN + NaN i */
+		return (cpack(a, t));	/* return NaN + NaN i */
 	}
 	if (isinf(a)) {
 		/*
-		 * csqrt(inf + nan i)  = inf +  nan i
+		 * csqrt(inf + NaN i)  = inf +  NaN i
 		 * csqrt(inf + y i)    = inf +  0 i
-		 * csqrt(-inf + nan i) = nan +- inf i
+		 * csqrt(-inf + NaN i) = NaN +- inf i
 		 * csqrt(-inf + y i)   = 0   +  inf i
 		 */
 		if (signbit(a))
@@ -79,15 +82,15 @@ csqrt(double complex z)
 	 */
 
 	/* Scale to avoid overflow. */
-	if (a >= THRESH || b >= THRESH) {
-	    a *= 0.25;
-	    b *= 0.25;
-	    scale = 1;
+	if (fabs(a) >= THRESH || fabs(b) >= THRESH) {
+		a *= 0.25;
+		b *= 0.25;
+		scale = 1;
 	} else {
-	    scale = 0;
+		scale = 0;
 	}
 
-	/* Algorithm 312, CACM vol 10, Oct 1967 */
+	/* Algorithm 312, CACM vol 10, Oct 1967. */
 	if (a >= 0) {
 		t = sqrt((a + hypot(a, b)) * 0.5);
 		result = cpack(t, b / (2 * t));
@@ -96,9 +99,9 @@ csqrt(double complex z)
 		result = cpack(fabs(b) / (2 * t), copysign(t, b));
 	}
 
-	/* Rescale */
+	/* Rescale. */
 	if (scale)
-	    return (result * 2);
+		return (result * 2);
 	else
-	    return (result);
+		return (result);
 }
diff --git a/lib/msun/src/s_csqrtf.c b/lib/msun/src/s_csqrtf.c
index cc5485a..6c70ac1 100644
--- a/lib/msun/src/s_csqrtf.c
+++ b/lib/msun/src/s_csqrtf.c
@@ -54,13 +54,13 @@ csqrtf(float complex z)
 		return (cpackf(INFINITY, b));
 	if (isnan(a)) {
 		t = (b - b) / (b - b);	/* raise invalid if b is not a NaN */
-		return (cpackf(t, t));	/* return NaN + NaN i */
+		return (cpackf(a, t));	/* return NaN + NaN i */
 	}
 	if (isinf(a)) {
 		/*
-		 * csqrtf(inf + nan i)  = inf +  nan i
+		 * csqrtf(inf + NaN i)  = inf +  NaN i
 		 * csqrtf(inf + y i)    = inf +  0 i
-		 * csqrtf(-inf + nan i) = nan +- inf i
+		 * csqrtf(-inf + NaN i) = NaN +- inf i
 		 * csqrtf(-inf + y i)   = 0   +  inf i
 		 */
 		if (signbit(a))