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