1 #include <math.h> 2 #include <stdint.h> 3 #include <float.h> 4 5 #if FLT_EVAL_METHOD > 1U && LDBL_MANT_DIG == 64 6 #define SPLIT (0x1p32 + 1) 7 #else 8 #define SPLIT (0x1p27 + 1) 9 #endif 10 11 static void sq(double_t *hi, double_t *lo, double x) 12 { 13 double_t xh, xl, xc; 14 15 xc = (double_t)x*SPLIT; 16 xh = x - xc + xc; 17 xl = x - xh; 18 *hi = (double_t)x*x; 19 *lo = xh*xh - *hi + 2*xh*xl + xl*xl; 20 } 21 22 double hypot(double x, double y) 23 { 24 union {double f; uint64_t i;} ux = {x}, uy = {y}, ut; 25 int ex, ey; 26 double_t hx, lx, hy, ly, z; 27 28 /* arrange |x| >= |y| */ 29 ux.i &= -1ULL>>1; 30 uy.i &= -1ULL>>1; 31 if (ux.i < uy.i) { 32 ut = ux; 33 ux = uy; 34 uy = ut; 35 } 36 37 /* special cases */ 38 ex = ux.i>>52; 39 ey = uy.i>>52; 40 x = ux.f; 41 y = uy.f; 42 /* note: hypot(inf,nan) == inf */ 43 if (ey == 0x7ff) 44 return y; 45 if (ex == 0x7ff || uy.i == 0) 46 return x; 47 /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */ 48 /* 64 difference is enough for ld80 double_t */ 49 if (ex - ey > 64) 50 return x + y; 51 52 /* precise sqrt argument in nearest rounding mode without overflow */ 53 /* xh*xh must not overflow and xl*xl must not underflow in sq */ 54 z = 1; 55 if (ex > 0x3ff+510) { 56 z = 0x1p700; 57 x *= 0x1p-700; 58 y *= 0x1p-700; 59 } else if (ey < 0x3ff-450) { 60 z = 0x1p-700; 61 x *= 0x1p700; 62 y *= 0x1p700; 63 } 64 sq(&hx, &lx, x); 65 sq(&hy, &ly, y); 66 return z*sqrt(ly+lx+hy+hx); 67 } 68