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 |