7  * Developed at SunPro, a Sun Microsystems, Inc. business.
17  *	If (assume round-to-nearest) z=x*x+y*y
25  *	(if possible, set rounding to round-to-nearest)
28  *	where x1 = x with lower 32 bits cleared, x2 = x-x1; else
30  *		t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
31  *	where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
32  *	y1= y with lower 32 bits chopped, y2 = y-y1.
38  *	hypot(x,y) is INF if x or y is +INF or -INF; else
53     __float64 a = x, b = y, t1, t2, y1, y2, w;  in hypot64()  local
61         a = y;  in hypot64()
67         a = x;  in hypot64()
70     SET_HIGH_WORD(a, ha); /* a <- |a| */  in hypot64()
71     SET_HIGH_WORD(b, hb); /* b <- |b| */  in hypot64()
72     if ((ha - hb) > 0x3c00000) {  in hypot64()
73         return a + b;  in hypot64()
76     if (ha > 0x5f300000) { /* a>2**500 */  in hypot64()
79             w = a + b; /* for sNaN */  in hypot64()
80             GET_LOW_WORD(low, a);  in hypot64()
82                 w = a;  in hypot64()
84             if (((hb ^ 0x7ff00000) | low) == 0 && !issignaling(a))  in hypot64()
88         /* scale a and b by 2**-600 */  in hypot64()
89         ha -= 0x25800000;  in hypot64()
90         hb -= 0x25800000;  in hypot64()
92         SET_HIGH_WORD(a, ha);  in hypot64()
95     if (hb < 0x20b00000) { /* b < 2**-500 */  in hypot64()
100                 return a;  in hypot64()
104             a *= t1;  in hypot64()
105             k -= 1022;  in hypot64()
106         } else { /* scale a and b by 2^600 */  in hypot64()
107             ha += 0x25800000; /* a *= 2^600 */  in hypot64()
109             k -= 600;  in hypot64()
110             SET_HIGH_WORD(a, ha);  in hypot64()
114     /* medium size a and b */  in hypot64()
115     w = a - b;  in hypot64()
119         t2 = a - t1;  in hypot64()
120         w = sqrt64(t1 * t1 - (b * (-b) - t2 * (a + t1)));  in hypot64()
122         a = a + a;  in hypot64()
125         y2 = b - y1;  in hypot64()
128         t2 = a - t1;  in hypot64()
129         w = sqrt64(t1 * y1 - (w * (-w) - (t1 * y2 + t2 * b)));  in hypot64()