| #include "libm.h" |
| |
| #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 |
| long double hypotl(long double x, long double y) { |
| return hypot(x, y); |
| } |
| #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 |
| #if LDBL_MANT_DIG == 64 |
| #define SPLIT (0x1p32L + 1) |
| #elif LDBL_MANT_DIG == 113 |
| #define SPLIT (0x1p57L + 1) |
| #endif |
| |
| static void sq(long double* hi, long double* lo, long double x) { |
| long double xh, xl, xc; |
| xc = x * SPLIT; |
| xh = x - xc + xc; |
| xl = x - xh; |
| *hi = x * x; |
| *lo = xh * xh - *hi + 2 * xh * xl + xl * xl; |
| } |
| |
| long double hypotl(long double x, long double y) { |
| union ldshape ux = {x}, uy = {y}; |
| int ex, ey; |
| long double hx, lx, hy, ly, z; |
| |
| ux.i.se &= 0x7fff; |
| uy.i.se &= 0x7fff; |
| if (ux.i.se < uy.i.se) { |
| ex = uy.i.se; |
| ey = ux.i.se; |
| x = uy.f; |
| y = ux.f; |
| } else { |
| ex = ux.i.se; |
| ey = uy.i.se; |
| x = ux.f; |
| y = uy.f; |
| } |
| |
| if (ex == 0x7fff && isinf(y)) |
| return y; |
| if (ex == 0x7fff || y == 0) |
| return x; |
| if (ex - ey > LDBL_MANT_DIG) |
| return x + y; |
| |
| z = 1; |
| if (ex > 0x3fff + 8000) { |
| z = 0x1p10000L; |
| x *= 0x1p-10000L; |
| y *= 0x1p-10000L; |
| } else if (ey < 0x3fff - 8000) { |
| z = 0x1p-10000L; |
| x *= 0x1p10000L; |
| y *= 0x1p10000L; |
| } |
| sq(&hx, &lx, x); |
| sq(&hy, &ly, y); |
| return z * sqrtl(ly + lx + hy + hx); |
| } |
| #endif |
