1 /* @(#)e_hypotl.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_hypotl(x,y)
16 * If (assume round-to-nearest) z=x*x+y*y
17 * has error less than sqrtl(2)/2 ulp, than
18 * sqrtl(z) has error less than 1 ulp (exercise).
20 * So, compute sqrtl(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 53 bits cleared, x2 = x-x1; else
29 * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
30 * where t1 = 2x with lower 53 bits cleared, t2 = 2x-t1,
31 * y1= y with lower 53 bits chopped, y2 = y-y1.
33 * NOTE: scaling may be necessary if some argument is too
37 * hypotl(x,y) is INF if x or y is +INF or -INF; else
38 * hypotl(x,y) is NAN if x or y is NAN.
41 * hypotl(x,y) returns sqrtl(x^2+y^2) with error less
42 * than 1 ulps (units in the last place)
46 #include <math_private.h>
47 #include <math-underflow.h>
48 #include <libm-alias-finite.h>
51 __ieee754_hypotl(long double x
, long double y
)
53 long double a
,b
,a1
,a2
,b1
,b2
,w
,kld
;
55 double xhi
, yhi
, hi
, lo
;
58 EXTRACT_WORDS64 (ha
, xhi
);
60 EXTRACT_WORDS64 (hb
, yhi
);
61 ha
&= 0x7fffffffffffffffLL
;
62 hb
&= 0x7fffffffffffffffLL
;
63 if(hb
> ha
) {a
=y
;b
=x
;j
=ha
; ha
=hb
;hb
=j
;} else {a
=x
;b
=y
;}
64 a
= fabsl(a
); /* a <- |a| */
65 b
= fabsl(b
); /* b <- |b| */
66 if((ha
-hb
)>0x0780000000000000LL
) {return a
+b
;} /* x/y > 2**120 */
69 if(ha
> 0x5f30000000000000LL
) { /* a>2**500 */
70 if(ha
>= 0x7ff0000000000000LL
) { /* Inf or NaN */
71 w
= a
+b
; /* for sNaN */
72 if (issignaling (a
) || issignaling (b
))
74 if(ha
== 0x7ff0000000000000LL
)
76 if(hb
== 0x7ff0000000000000LL
)
80 /* scale a and b by 2**-600 */
86 else if(hb
< 0x23d0000000000000LL
) { /* b < 2**-450 */
87 if(hb
<= 0x000fffffffffffffLL
) { /* subnormal b or 0 */
93 } else { /* scale a and b by 2^600 */
100 /* medium size a and b */
103 ldbl_unpack (a
, &hi
, &lo
);
109 = a1*(a1+a2) + a2*a + b*b
110 = a1*a1 + a1*a2 + a2*a + b*b
111 = a1*a1 + a2*(a+a1) + b*b */
112 w
= sqrtl(a1
*a1
-(b
*(-b
)-a2
*(a
+a1
)));
115 ldbl_unpack (b
, &hi
, &lo
);
118 ldbl_unpack (a
, &hi
, &lo
);
122 = a*a + (a-b)*(a-b) - (a-b)*(a-b) + b*b
123 = a*a + w*w - (a*a - 2*a*b + b*b) + b*b
127 = w*w + a1*(b1+b2) + a2*b
128 = w*w + a1*b1 + a1*b2 + a2*b */
129 w
= sqrtl(a1
*b1
-(w
*(-w
)-(a1
*b2
+a2
*b
)));
134 math_check_force_underflow_nonneg (w
);
140 libm_alias_finite (__ieee754_hypotl
, __hypotl
)