1 /* @(#)e_hypot.c 5.1 93/09/24 */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Developed at SunPro, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
10 * ====================================================
13 /* __ieee754_hypot(x,y)
16 * If (assume round-to-nearest) z=x*x+y*y
17 * has error less than sqrt(2)/2 ulp, than
18 * sqrt(z) has error less than 1 ulp (exercise).
20 * So, compute sqrt(x*x+y*y) with some care as
21 * follows to get the error below 1 ulp:
24 * (if possible, set rounding to round-to-nearest)
26 * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
27 * where x1 = x with lower 32 bits cleared, x2 = x-x1; else
29 * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
30 * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
31 * y1= y with lower 32 bits chopped, y2 = y-y1.
33 * NOTE: scaling may be necessary if some argument is too
37 * hypot(x,y) is INF if x or y is +INF or -INF; else
38 * hypot(x,y) is NAN if x or y is NAN.
41 * hypot(x,y) returns sqrt(x^2+y^2) with error less
42 * than 1 ulps (units in the last place)
46 #include <math_private.h>
49 __ieee754_hypot (double x
, double y
)
51 double a
, b
, t1
, t2
, y1
, y2
, w
;
54 GET_HIGH_WORD (ha
, x
);
56 GET_HIGH_WORD (hb
, y
);
60 a
= y
; b
= x
; j
= ha
; ha
= hb
; hb
= j
;
66 SET_HIGH_WORD (a
, ha
); /* a <- |a| */
67 SET_HIGH_WORD (b
, hb
); /* b <- |b| */
68 if ((ha
- hb
) > 0x3c00000)
73 if (__glibc_unlikely (ha
> 0x5f300000)) /* a>2**500 */
75 if (ha
>= 0x7ff00000) /* Inf or NaN */
78 w
= a
+ b
; /* for sNaN */
79 GET_LOW_WORD (low
, a
);
80 if (((ha
& 0xfffff) | low
) == 0)
82 GET_LOW_WORD (low
, b
);
83 if (((hb
^ 0x7ff00000) | low
) == 0)
87 /* scale a and b by 2**-600 */
88 ha
-= 0x25800000; hb
-= 0x25800000; k
+= 600;
89 SET_HIGH_WORD (a
, ha
);
90 SET_HIGH_WORD (b
, hb
);
92 if (__builtin_expect (hb
< 0x23d00000, 0)) /* b < 2**-450 */
94 if (hb
<= 0x000fffff) /* subnormal b or 0 */
97 GET_LOW_WORD (low
, b
);
101 SET_HIGH_WORD (t1
, 0x7fd00000); /* t1=2^1022 */
105 GET_HIGH_WORD (ha
, a
);
106 GET_HIGH_WORD (hb
, b
);
117 else /* scale a and b by 2^600 */
119 ha
+= 0x25800000; /* a *= 2^600 */
120 hb
+= 0x25800000; /* b *= 2^600 */
122 SET_HIGH_WORD (a
, ha
);
123 SET_HIGH_WORD (b
, hb
);
126 /* medium size a and b */
131 SET_HIGH_WORD (t1
, ha
);
133 w
= __ieee754_sqrt (t1
* t1
- (b
* (-b
) - t2
* (a
+ t1
)));
139 SET_HIGH_WORD (y1
, hb
);
142 SET_HIGH_WORD (t1
, ha
+ 0x00100000);
144 w
= __ieee754_sqrt (t1
* y1
- (w
* (-w
) - (t1
* y2
+ t2
* b
)));
150 GET_HIGH_WORD (high
, t1
);
151 SET_HIGH_WORD (t1
, high
+ (k
<< 20));
153 math_check_force_underflow_nonneg (w
);
159 strong_alias (__ieee754_hypot
, __hypot_finite
)