#include "libm.h"

#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
long double hypotl(long double x, long double y)
{
	return hypot(x, y);
}
#elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
#if LDBL_MANT_DIG == 64
#define SPLIT (0x1p32L + 1)
#elif LDBL_MANT_DIG == 113
#define SPLIT (0x1p57L + 1)
#endif

static void sq(long double *hi, long double *lo, long double x)
{
	long double xh, xl, xc;
	xc = x * SPLIT;
	xh = x - xc + xc;
	xl = x - xh;
	*hi = x * x;
	*lo = xh * xh - *hi + 2 * xh * xl + xl * xl;
}

long double hypotl(long double x, long double y)
{
	union ldshape ux = { x }, uy = { y };
	int ex, ey;
	long double hx, lx, hy, ly, z;

	ux.i.se &= 0x7fff;
	uy.i.se &= 0x7fff;
	if (ux.i.se < uy.i.se) {
		ex = uy.i.se;
		ey = ux.i.se;
		x = uy.f;
		y = ux.f;
	} else {
		ex = ux.i.se;
		ey = uy.i.se;
		x = ux.f;
		y = uy.f;
	}

	if (ex == 0x7fff && isinf(y))
		return y;
	if (ex == 0x7fff || y == 0)
		return x;
	if (ex - ey > LDBL_MANT_DIG)
		return x + y;

	z = 1;
	if (ex > 0x3fff + 8000) {
		z = 0x1p10000L;
		x *= 0x1p-10000L;
		y *= 0x1p-10000L;
	} else if (ey < 0x3fff - 8000) {
		z = 0x1p-10000L;
		x *= 0x1p10000L;
		y *= 0x1p10000L;
	}
	sq(&hx, &lx, x);
	sq(&hy, &ly, y);
	return z * sqrtl(ly + lx + hy + hx);
}
#endif