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
11static 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
21double 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