Szabolcs Nagy
c2a0dfea62
method: if there is a large difference between the scale of x and y then the larger magnitude dominates, otherwise reduce x,y so the argument of sqrt (x*x+y*y) does not overflow or underflow and calculate the argument precisely using exact multiplication. If the argument has less error than 1/sqrt(2) ~ 0.7 ulp, then the result has less error than 1 ulp in nearest rounding mode. the original fdlibm method was the same, except it used bit hacks instead of dekker-veltkamp algorithm, which is problematic for long double where different representations are supported. (the new hypot and hypotl code should be smaller and faster on 32bit cpu archs with fast fpu), the new code behaves differently in non-nearest rounding, but the error should be still less than 2ulps. ld80 and ld128 are supported
67 lines
1.2 KiB
C
67 lines
1.2 KiB
C
#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
|