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