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