]>
git.ipfire.org Git - thirdparty/glibc.git/blob - sysdeps/ieee754/ldbl-96/e_hypotl.c
1 /* e_hypotl.c -- long double version of e_hypot.c.
5 * ====================================================
6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 * Developed at SunPro, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
12 * ====================================================
15 /* __ieee754_hypotl(x,y)
18 * If (assume round-to-nearest) z=x*x+y*y
19 * has error less than sqrt(2)/2 ulp, than
20 * sqrt(z) has error less than 1 ulp (exercise).
22 * So, compute sqrt(x*x+y*y) with some care as
23 * follows to get the error below 1 ulp:
26 * (if possible, set rounding to round-to-nearest)
28 * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
29 * where x1 = x with lower 32 bits cleared, x2 = x-x1; else
31 * t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
32 * where t1 = 2x with lower 32 bits cleared, t2 = 2x-t1,
33 * y1= y with lower 32 bits chopped, y2 = y-y1.
35 * NOTE: scaling may be necessary if some argument is too
39 * hypot(x,y) is INF if x or y is +INF or -INF; else
40 * hypot(x,y) is NAN if x or y is NAN.
43 * hypot(x,y) returns sqrt(x^2+y^2) with error less
44 * than 1 ulps (units in the last place)
48 #include <math_private.h>
49 #include <math-underflow.h>
50 #include <libm-alias-finite.h>
52 long double __ieee754_hypotl(long double x
, long double y
)
54 long double a
,b
,t1
,t2
,y1
,y2
,w
;
57 GET_LDOUBLE_EXP(ea
,x
);
59 GET_LDOUBLE_EXP(eb
,y
);
61 if(eb
> ea
) {a
=y
;b
=x
;j
=ea
; ea
=eb
;eb
=j
;} else {a
=x
;b
=y
;}
62 SET_LDOUBLE_EXP(a
,ea
); /* a <- |a| */
63 SET_LDOUBLE_EXP(b
,eb
); /* b <- |b| */
64 if((ea
-eb
)>0x46) {return a
+b
;} /* x/y > 2**70 */
66 if(__builtin_expect(ea
> 0x5f3f,0)) { /* a>2**8000 */
67 if(ea
== 0x7fff) { /* Inf or NaN */
68 uint32_t exp
__attribute__ ((unused
));
70 w
= a
+b
; /* for sNaN */
71 if (issignaling (a
) || issignaling (b
))
73 GET_LDOUBLE_WORDS(exp
,high
,low
,a
);
74 if(((high
&0x7fffffff)|low
)==0) w
= a
;
75 GET_LDOUBLE_WORDS(exp
,high
,low
,b
);
76 if(((eb
^0x7fff)|(high
&0x7fffffff)|low
)==0) w
= b
;
79 /* scale a and b by 2**-9600 */
80 ea
-= 0x2580; eb
-= 0x2580; k
+= 9600;
81 SET_LDOUBLE_EXP(a
,ea
);
82 SET_LDOUBLE_EXP(b
,eb
);
84 if(__builtin_expect(eb
< 0x20bf, 0)) { /* b < 2**-8000 */
85 if(eb
== 0) { /* subnormal b or 0 */
86 uint32_t exp
__attribute__ ((unused
));
88 GET_LDOUBLE_WORDS(exp
,high
,low
,b
);
89 if((high
|low
)==0) return a
;
90 SET_LDOUBLE_WORDS(t1
, 0x7ffd, 0x80000000, 0); /* t1=2^16382 */
94 GET_LDOUBLE_EXP (ea
, a
);
95 GET_LDOUBLE_EXP (eb
, b
);
105 } else { /* scale a and b by 2^9600 */
106 ea
+= 0x2580; /* a *= 2^9600 */
107 eb
+= 0x2580; /* b *= 2^9600 */
109 SET_LDOUBLE_EXP(a
,ea
);
110 SET_LDOUBLE_EXP(b
,eb
);
113 /* medium size a and b */
117 GET_LDOUBLE_MSW(high
,a
);
118 SET_LDOUBLE_WORDS(t1
,ea
,high
,0);
120 w
= sqrtl(t1
*t1
-(b
*(-b
)-t2
*(a
+t1
)));
123 GET_LDOUBLE_MSW(high
,b
);
125 SET_LDOUBLE_WORDS(y1
,eb
,high
,0);
127 GET_LDOUBLE_MSW(high
,a
);
128 SET_LDOUBLE_WORDS(t1
,ea
+1,high
,0);
130 w
= sqrtl(t1
*y1
-(w
*(-w
)-(t1
*y2
+t2
*b
)));
135 GET_LDOUBLE_EXP(exp
,t1
);
136 SET_LDOUBLE_EXP(t1
,exp
+k
);
138 math_check_force_underflow_nonneg (w
);
142 libm_alias_finite (__ieee754_hypotl
, __hypotl
)