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{
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
15static 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
25long 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