1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
|
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunPro, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* From e_hypotl.c -- long double version of e_hypot.c.
* Conversion to long double by Jakub Jelinek, jakub@redhat.com.
* Conversion to __float128 by FX Coudert, fxcoudert@gcc.gnu.org.
*/
/* hypotq(x,y)
*
* Method :
* If (assume round-to-nearest) z=x*x+y*y
* has error less than sqrtl(2)/2 ulp, than
* sqrtl(z) has error less than 1 ulp (exercise).
*
* So, compute sqrtl(x*x+y*y) with some care as
* follows to get the error below 1 ulp:
*
* Assume x>y>0;
* (if possible, set rounding to round-to-nearest)
* 1. if x > 2y use
* x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
* where x1 = x with lower 64 bits cleared, x2 = x-x1; else
* 2. if x <= 2y use
* t1*y1+((x-y)*(x-y)+(t1*y2+t2*y))
* where t1 = 2x with lower 64 bits cleared, t2 = 2x-t1,
* y1= y with lower 64 bits chopped, y2 = y-y1.
*
* NOTE: scaling may be necessary if some argument is too
* large or too tiny
*
* Special cases:
* hypotq(x,y) is INF if x or y is +INF or -INF; else
* hypotq(x,y) is NAN if x or y is NAN.
*
* Accuracy:
* hypotq(x,y) returns sqrtl(x^2+y^2) with error less
* than 1 ulps (units in the last place)
*/
#include "quadmath-imp.h"
__float128
hypotq (__float128 x, __float128 y)
{
__float128 a, b, t1, t2, y1, y2, w;
int64_t j, k, ha, hb;
GET_FLT128_MSW64(ha,x);
ha &= 0x7fffffffffffffffLL;
GET_FLT128_MSW64(hb,y);
hb &= 0x7fffffffffffffffLL;
if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
SET_FLT128_MSW64(a,ha); /* a <- |a| */
SET_FLT128_MSW64(b,hb); /* b <- |b| */
if((ha-hb)>0x78000000000000LL) {return a+b;} /* x/y > 2**120 */
k=0;
if(ha > 0x5f3f000000000000LL) { /* a>2**8000 */
if(ha >= 0x7fff000000000000LL) { /* Inf or NaN */
uint64_t low;
w = a+b; /* for sNaN */
GET_FLT128_LSW64(low,a);
if(((ha&0xffffffffffffLL)|low)==0) w = a;
GET_FLT128_LSW64(low,b);
if(((hb^0x7fff000000000000LL)|low)==0) w = b;
return w;
}
/* scale a and b by 2**-9600 */
ha -= 0x2580000000000000LL;
hb -= 0x2580000000000000LL; k += 9600;
SET_FLT128_MSW64(a,ha);
SET_FLT128_MSW64(b,hb);
}
if(hb < 0x20bf000000000000LL) { /* b < 2**-8000 */
if(hb <= 0x0000ffffffffffffLL) { /* subnormal b or 0 */
uint64_t low;
GET_FLT128_LSW64(low,b);
if((hb|low)==0) return a;
t1=0;
SET_FLT128_MSW64(t1,0x7ffd000000000000LL); /* t1=2^16382 */
b *= t1;
a *= t1;
k -= 16382;
GET_FLT128_MSW64 (ha, a);
GET_FLT128_MSW64 (hb, b);
if (hb > ha)
{
t1 = a;
a = b;
b = t1;
j = ha;
ha = hb;
hb = j;
}
} else { /* scale a and b by 2^9600 */
ha += 0x2580000000000000LL; /* a *= 2^9600 */
hb += 0x2580000000000000LL; /* b *= 2^9600 */
k -= 9600;
SET_FLT128_MSW64(a,ha);
SET_FLT128_MSW64(b,hb);
}
}
/* medium size a and b */
w = a-b;
if (w>b) {
t1 = 0;
SET_FLT128_MSW64(t1,ha);
t2 = a-t1;
w = sqrtq(t1*t1-(b*(-b)-t2*(a+t1)));
} else {
a = a+a;
y1 = 0;
SET_FLT128_MSW64(y1,hb);
y2 = b - y1;
t1 = 0;
SET_FLT128_MSW64(t1,ha+0x0001000000000000LL);
t2 = a - t1;
w = sqrtq(t1*y1-(w*(-w)-(t1*y2+t2*b)));
}
if(k!=0) {
uint64_t high;
t1 = 1.0Q;
GET_FLT128_MSW64(high,t1);
SET_FLT128_MSW64(t1,high+(k<<48));
w *= t1;
math_check_force_underflow_nonneg (w);
return w;
} else return w;
}
|
