1 #include "libm.h"
2
3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
hypotl(long double x,long double y)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
sq(long double * hi,long double * lo,long double x)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
hypotl(long double x,long double y)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
67