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