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