Index: fusl/src/math/hypot.c |
diff --git a/fusl/src/math/hypot.c b/fusl/src/math/hypot.c |
index 6071bf1e284f376f417dbab5bc306964a7491f3f..5f9888d88433e0bd81366e5fafde9ce93c0b2cec 100644 |
--- a/fusl/src/math/hypot.c |
+++ b/fusl/src/math/hypot.c |
@@ -8,60 +8,61 @@ |
#define SPLIT (0x1p27 + 1) |
#endif |
-static void sq(double_t *hi, double_t *lo, double x) |
-{ |
- double_t xh, xl, xc; |
+static void sq(double_t* hi, double_t* lo, double x) { |
+ double_t xh, xl, xc; |
- xc = (double_t)x*SPLIT; |
- xh = x - xc + xc; |
- xl = x - xh; |
- *hi = (double_t)x*x; |
- *lo = xh*xh - *hi + 2*xh*xl + xl*xl; |
+ xc = (double_t)x * SPLIT; |
+ xh = x - xc + xc; |
+ xl = x - xh; |
+ *hi = (double_t)x * x; |
+ *lo = xh * xh - *hi + 2 * xh * xl + xl * xl; |
} |
-double hypot(double x, double y) |
-{ |
- union {double f; uint64_t i;} ux = {x}, uy = {y}, ut; |
- int ex, ey; |
- double_t hx, lx, hy, ly, z; |
+double hypot(double x, double y) { |
+ union { |
+ double f; |
+ uint64_t i; |
+ } ux = {x}, uy = {y}, ut; |
+ int ex, ey; |
+ double_t hx, lx, hy, ly, z; |
- /* arrange |x| >= |y| */ |
- ux.i &= -1ULL>>1; |
- uy.i &= -1ULL>>1; |
- if (ux.i < uy.i) { |
- ut = ux; |
- ux = uy; |
- uy = ut; |
- } |
+ /* arrange |x| >= |y| */ |
+ ux.i &= -1ULL >> 1; |
+ uy.i &= -1ULL >> 1; |
+ if (ux.i < uy.i) { |
+ ut = ux; |
+ ux = uy; |
+ uy = ut; |
+ } |
- /* special cases */ |
- ex = ux.i>>52; |
- ey = uy.i>>52; |
- x = ux.f; |
- y = uy.f; |
- /* note: hypot(inf,nan) == inf */ |
- if (ey == 0x7ff) |
- return y; |
- if (ex == 0x7ff || uy.i == 0) |
- return x; |
- /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */ |
- /* 64 difference is enough for ld80 double_t */ |
- if (ex - ey > 64) |
- return x + y; |
+ /* special cases */ |
+ ex = ux.i >> 52; |
+ ey = uy.i >> 52; |
+ x = ux.f; |
+ y = uy.f; |
+ /* note: hypot(inf,nan) == inf */ |
+ if (ey == 0x7ff) |
+ return y; |
+ if (ex == 0x7ff || uy.i == 0) |
+ return x; |
+ /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */ |
+ /* 64 difference is enough for ld80 double_t */ |
+ if (ex - ey > 64) |
+ return x + y; |
- /* precise sqrt argument in nearest rounding mode without overflow */ |
- /* xh*xh must not overflow and xl*xl must not underflow in sq */ |
- z = 1; |
- if (ex > 0x3ff+510) { |
- z = 0x1p700; |
- x *= 0x1p-700; |
- y *= 0x1p-700; |
- } else if (ey < 0x3ff-450) { |
- z = 0x1p-700; |
- x *= 0x1p700; |
- y *= 0x1p700; |
- } |
- sq(&hx, &lx, x); |
- sq(&hy, &ly, y); |
- return z*sqrt(ly+lx+hy+hx); |
+ /* precise sqrt argument in nearest rounding mode without overflow */ |
+ /* xh*xh must not overflow and xl*xl must not underflow in sq */ |
+ z = 1; |
+ if (ex > 0x3ff + 510) { |
+ z = 0x1p700; |
+ x *= 0x1p-700; |
+ y *= 0x1p-700; |
+ } else if (ey < 0x3ff - 450) { |
+ z = 0x1p-700; |
+ x *= 0x1p700; |
+ y *= 0x1p700; |
+ } |
+ sq(&hx, &lx, x); |
+ sq(&hy, &ly, y); |
+ return z * sqrt(ly + lx + hy + hx); |
} |