1#include "libm.h" 2 3#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 4long double hypotl(long double x, long double y) { 5 return hypot(x, y); 6} 7#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 8#if LDBL_MANT_DIG == 64 9#define SPLIT (0x1p32L + 1) 10#elif LDBL_MANT_DIG == 113 11#define SPLIT (0x1p57L + 1) 12#endif 13 14static void sq(long double* hi, long double* lo, long double x) { 15 long double xh, xl, xc; 16 xc = x * SPLIT; 17 xh = x - xc + xc; 18 xl = x - xh; 19 *hi = x * x; 20 *lo = xh * xh - *hi + 2 * xh * xl + xl * xl; 21} 22 23long double hypotl(long double x, long double y) { 24 union ldshape ux = {x}, uy = {y}; 25 int ex, ey; 26 long double hx, lx, hy, ly, z; 27 28 ux.i.se &= 0x7fff; 29 uy.i.se &= 0x7fff; 30 if (ux.i.se < uy.i.se) { 31 ex = uy.i.se; 32 ey = ux.i.se; 33 x = uy.f; 34 y = ux.f; 35 } else { 36 ex = ux.i.se; 37 ey = uy.i.se; 38 x = ux.f; 39 y = uy.f; 40 } 41 42 if (ex == 0x7fff && isinf(y)) 43 return y; 44 if (ex == 0x7fff || y == 0) 45 return x; 46 if (ex - ey > LDBL_MANT_DIG) 47 return x + y; 48 49 z = 1; 50 if (ex > 0x3fff + 8000) { 51 z = 0x1p10000L; 52 x *= 0x1p-10000L; 53 y *= 0x1p-10000L; 54 } else if (ey < 0x3fff - 8000) { 55 z = 0x1p-10000L; 56 x *= 0x1p10000L; 57 y *= 0x1p10000L; 58 } 59 sq(&hx, &lx, x); 60 sq(&hy, &ly, y); 61 return z * sqrtl(ly + lx + hy + hx); 62} 63#endif 64