1 #include <float.h>
2 #include <math.h>
3 #include <stdint.h>
4
5 #if FLT_EVAL_METHOD > 1U && LDBL_MANT_DIG == 64
6 #define SPLIT (0x1p32 + 1)
7 #else
8 #define SPLIT (0x1p27 + 1)
9 #endif
10
sq(double_t * hi,double_t * lo,double x)11 static void sq(double_t* hi, double_t* lo, double x) {
12 double_t xh, xl, xc;
13
14 xc = (double_t)x * SPLIT;
15 xh = x - xc + xc;
16 xl = x - xh;
17 *hi = (double_t)x * x;
18 *lo = xh * xh - *hi + 2 * xh * xl + xl * xl;
19 }
20
hypot(double x,double y)21 double hypot(double x, double y) {
22 union {
23 double f;
24 uint64_t i;
25 } ux = {x}, uy = {y}, ut;
26 int ex, ey;
27 double_t hx, lx, hy, ly, z;
28
29 /* arrange |x| >= |y| */
30 ux.i &= -1ULL >> 1;
31 uy.i &= -1ULL >> 1;
32 if (ux.i < uy.i) {
33 ut = ux;
34 ux = uy;
35 uy = ut;
36 }
37
38 /* special cases */
39 ex = ux.i >> 52;
40 ey = uy.i >> 52;
41 x = ux.f;
42 y = uy.f;
43 /* note: hypot(inf,nan) == inf */
44 if (ey == 0x7ff)
45 return y;
46 if (ex == 0x7ff || uy.i == 0)
47 return x;
48 /* note: hypot(x,y) ~= x + y*y/x/2 with inexact for small y/x */
49 /* 64 difference is enough for ld80 double_t */
50 if (ex - ey > 64)
51 return x + y;
52
53 /* precise sqrt argument in nearest rounding mode without overflow */
54 /* xh*xh must not overflow and xl*xl must not underflow in sq */
55 z = 1;
56 if (ex > 0x3ff + 510) {
57 z = 0x1p700;
58 x *= 0x1p-700;
59 y *= 0x1p-700;
60 } else if (ey < 0x3ff - 450) {
61 z = 0x1p-700;
62 x *= 0x1p700;
63 y *= 0x1p700;
64 }
65 sq(&hx, &lx, x);
66 sq(&hy, &ly, y);
67 return z * sqrt(ly + lx + hy + hx);
68 }
69