summaryrefslogtreecommitdiff log msg author committer range
path: root/src/math/hypot.c
blob: 6071bf1e284f376f417dbab5bc306964a7491f3f (plain)
 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 #include #include #include #if FLT_EVAL_METHOD > 1U && LDBL_MANT_DIG == 64 #define SPLIT (0x1p32 + 1) #else #define SPLIT (0x1p27 + 1) #endif 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; } 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; } /* 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); }