#include <math.h>
#include <stdint.h>
typedef union{ uint64_t u; double d; } du;
double hypot( double x, double y )
{
static const double inf = __builtin_inf();
du u[3];
du *large = u;
du *small = &u[1];
u[0].d = __builtin_fabs(x);
u[1].d = __builtin_fabs(y);
if( 0x7ff0000000000000ULL == ( u[0].u & 0x7ff0000000000000ULL) ||
0x7ff0000000000000ULL == ( u[1].u & 0x7ff0000000000000ULL) )
{
if( 0x7ff0000000000000ULL == u[0].u || 0x7ff0000000000000ULL == u[1].u )
return inf;
return x + y; }
if( x == 0.0 || y == 0.0 )
return __builtin_fabs( x + y );
if( u[0].d < u[1].d )
{
large = &u[1];
small = &u[0];
}
int64_t largeExp = large->u >> 52;
int64_t smallExp = small->u >> 52;
int64_t diff = largeExp - smallExp;
if( diff >= 55L )
return large->d + small->d;
large->u &= 0x000fffffffffffffULL;
small->u &= 0x000fffffffffffffULL;
large->u |= 0x3ff0000000000000ULL;
small->u |= 0x3ff0000000000000ULL;
if( 0 == smallExp )
{
if( 0 == largeExp )
{
large->d -= 1.0;
largeExp = (large->u >> 52) - (1022);
large->u &= 0x000fffffffffffffULL;
large->u |= 0x3ff0000000000000ULL;
}
small->d -= 1.0;
smallExp = (small->u >> 52) - (1022);
small->u &= 0x000fffffffffffffULL;
small->u |= 0x3ff0000000000000ULL;
}
u[2].u = (1023ULL - largeExp + smallExp) << 52;
small->d *= u[2].d;
double r = __builtin_sqrt( large->d * large->d + small->d * small->d );
if( largeExp < 0 )
{
largeExp += 1022;
r *= 0x1.0p-1022;
}
u[2].u = largeExp << 52;
r *= u[2].d;
return r;
}