Index: fusl/src/complex/csqrt.c |
diff --git a/fusl/src/complex/csqrt.c b/fusl/src/complex/csqrt.c |
index 8a2ba608012230d79a16ce788feb96e4b4d53a0f..951ab4f53bb032ef907031bccfa4a4505675b697 100644 |
--- a/fusl/src/complex/csqrt.c |
+++ b/fusl/src/complex/csqrt.c |
@@ -37,64 +37,63 @@ |
#pragma STDC CX_LIMITED_RANGE ON |
/* We risk spurious overflow for components >= DBL_MAX / (1 + sqrt(2)). */ |
-#define THRESH 0x1.a827999fcef32p+1022 |
+#define THRESH 0x1.a827999fcef32p+1022 |
-double complex csqrt(double complex z) |
-{ |
- double complex result; |
- double a, b; |
- double t; |
- int scale; |
+double complex csqrt(double complex z) { |
+ double complex result; |
+ double a, b; |
+ double t; |
+ int scale; |
- a = creal(z); |
- b = cimag(z); |
+ a = creal(z); |
+ b = cimag(z); |
- /* Handle special cases. */ |
- if (z == 0) |
- return CMPLX(0, b); |
- if (isinf(b)) |
- return CMPLX(INFINITY, b); |
- if (isnan(a)) { |
- t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ |
- return CMPLX(a, t); /* return NaN + NaN i */ |
- } |
- if (isinf(a)) { |
- /* |
- * csqrt(inf + NaN i) = inf + NaN i |
- * csqrt(inf + y i) = inf + 0 i |
- * csqrt(-inf + NaN i) = NaN +- inf i |
- * csqrt(-inf + y i) = 0 + inf i |
- */ |
- if (signbit(a)) |
- return CMPLX(fabs(b - b), copysign(a, b)); |
- else |
- return CMPLX(a, copysign(b - b, b)); |
- } |
- /* |
- * The remaining special case (b is NaN) is handled just fine by |
- * the normal code path below. |
- */ |
+ /* Handle special cases. */ |
+ if (z == 0) |
+ return CMPLX(0, b); |
+ if (isinf(b)) |
+ return CMPLX(INFINITY, b); |
+ if (isnan(a)) { |
+ t = (b - b) / (b - b); /* raise invalid if b is not a NaN */ |
+ return CMPLX(a, t); /* return NaN + NaN i */ |
+ } |
+ if (isinf(a)) { |
+ /* |
+ * csqrt(inf + NaN i) = inf + NaN i |
+ * csqrt(inf + y i) = inf + 0 i |
+ * csqrt(-inf + NaN i) = NaN +- inf i |
+ * csqrt(-inf + y i) = 0 + inf i |
+ */ |
+ if (signbit(a)) |
+ return CMPLX(fabs(b - b), copysign(a, b)); |
+ else |
+ return CMPLX(a, copysign(b - b, b)); |
+ } |
+ /* |
+ * The remaining special case (b is NaN) is handled just fine by |
+ * the normal code path below. |
+ */ |
- /* Scale to avoid overflow. */ |
- if (fabs(a) >= THRESH || fabs(b) >= THRESH) { |
- a *= 0.25; |
- b *= 0.25; |
- scale = 1; |
- } else { |
- scale = 0; |
- } |
+ /* Scale to avoid overflow. */ |
+ if (fabs(a) >= THRESH || fabs(b) >= THRESH) { |
+ a *= 0.25; |
+ b *= 0.25; |
+ scale = 1; |
+ } else { |
+ scale = 0; |
+ } |
- /* Algorithm 312, CACM vol 10, Oct 1967. */ |
- if (a >= 0) { |
- t = sqrt((a + hypot(a, b)) * 0.5); |
- result = CMPLX(t, b / (2 * t)); |
- } else { |
- t = sqrt((-a + hypot(a, b)) * 0.5); |
- result = CMPLX(fabs(b) / (2 * t), copysign(t, b)); |
- } |
+ /* Algorithm 312, CACM vol 10, Oct 1967. */ |
+ if (a >= 0) { |
+ t = sqrt((a + hypot(a, b)) * 0.5); |
+ result = CMPLX(t, b / (2 * t)); |
+ } else { |
+ t = sqrt((-a + hypot(a, b)) * 0.5); |
+ result = CMPLX(fabs(b) / (2 * t), copysign(t, b)); |
+ } |
- /* Rescale. */ |
- if (scale) |
- result *= 2; |
- return result; |
+ /* Rescale. */ |
+ if (scale) |
+ result *= 2; |
+ return result; |
} |