Loading...
#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 |