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