OLD | NEW |
(Empty) | |
| 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 |
OLD | NEW |