|
|
1
|
+
|
|
|
2
|
+/* @(#)e_hypot.c 1.3 95/01/18 */
|
|
|
3
|
+/*
|
|
|
4
|
+ * ====================================================
|
|
|
5
|
+ * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
|
|
|
6
|
+ *
|
|
|
7
|
+ * Developed at SunSoft, a Sun Microsystems, Inc. business.
|
|
|
8
|
+ * Permission to use, copy, modify, and distribute this
|
|
|
9
|
+ * software is freely granted, provided that this notice
|
|
|
10
|
+ * is preserved.
|
|
|
11
|
+ * ====================================================
|
|
|
12
|
+ */
|
|
|
13
|
+
|
|
|
14
|
+/* __ieee754_hypot(x,y)
|
|
|
15
|
+ *
|
|
|
16
|
+ * Method :
|
|
|
17
|
+ * If (assume round-to-nearest) z=x*x+y*y
|
|
|
18
|
+ * has error less than sqrt(2)/2 ulp, than
|
|
|
19
|
+ * sqrt(z) has error less than 1 ulp (exercise).
|
|
|
20
|
+ *
|
|
|
21
|
+ * So, compute sqrt(x*x+y*y) with some care as
|
|
|
22
|
+ * follows to get the error below 1 ulp:
|
|
|
23
|
+ *
|
|
|
24
|
+ * Assume x>y>0;
|
|
|
25
|
+ * (if possible, set rounding to round-to-nearest)
|
|
|
26
|
+ * 1. if x > 2y use
|
|
|
27
|
+ * x1*x1+(y*y+(x2*(x+x1))) for x*x+y*y
|
|
|
28
|
+ * where x1 = x with lower 32 bits cleared, x2 = x-x1; else
|
|
|
29
|
+ * 2. if x <= 2y use
|
|
|
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.
|
|
|
33
|
+ *
|
|
|
34
|
+ * NOTE: scaling may be necessary if some argument is too
|
|
|
35
|
+ * large or too tiny
|
|
|
36
|
+ *
|
|
|
37
|
+ * Special cases:
|
|
|
38
|
+ * hypot(x,y) is INF if x or y is +INF or -INF; else
|
|
|
39
|
+ * hypot(x,y) is NAN if x or y is NAN.
|
|
|
40
|
+ *
|
|
|
41
|
+ * Accuracy:
|
|
|
42
|
+ * hypot(x,y) returns sqrt(x^2+y^2) with error less
|
|
|
43
|
+ * than 1 ulps (units in the last place)
|
|
|
44
|
+ */
|
|
|
45
|
+
|
|
|
46
|
+#include "fdlibm.h"
|
|
|
47
|
+
|
|
|
48
|
+static inline void
|
|
|
49
|
+set_hiword(double* d, int hi)
|
|
|
50
|
+{
|
|
|
51
|
+ union { int i[2]; double d; } ux;
|
|
|
52
|
+
|
|
|
53
|
+ ux.d = *d;
|
|
|
54
|
+ ux.i[HIWORD] = hi;
|
|
|
55
|
+ *d = ux.d;
|
|
|
56
|
+}
|
|
|
57
|
+
|
|
|
58
|
+static inline int
|
|
|
59
|
+get_loword(double d)
|
|
|
60
|
+{
|
|
|
61
|
+ union { int i[2]; double d; } ux;
|
|
|
62
|
+ ux.d = d;
|
|
|
63
|
+ return ux.i[LOWORD];
|
|
|
64
|
+}
|
|
|
65
|
+
|
|
|
66
|
+static inline int
|
|
|
67
|
+get_hiword(double d)
|
|
|
68
|
+{
|
|
|
69
|
+ union { int i[2]; double d; } ux;
|
|
|
70
|
+ ux.d = d;
|
|
|
71
|
+
|
|
|
72
|
+ return ux.i[HIWORD];
|
|
|
73
|
+}
|
|
|
74
|
+
|
|
|
75
|
+
|
|
|
76
|
+#ifdef __STDC__
|
|
|
77
|
+ double __ieee754_hypot(double x, double y)
|
|
|
78
|
+#else
|
|
|
79
|
+ double __ieee754_hypot(x,y)
|
|
|
80
|
+ double x, y;
|
|
|
81
|
+#endif
|
|
|
82
|
+{
|
|
|
83
|
+ double a=x,b=y,t1,t2,y1,y2,w;
|
|
|
84
|
+ int j,k,ha,hb;
|
|
|
85
|
+
|
|
|
86
|
+ ha = get_hiword(x)&0x7fffffff; /* high word of x */
|
|
|
87
|
+ hb = get_hiword(y)&0x7fffffff; /* high word of y */
|
|
|
88
|
+ if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
|
|
|
89
|
+ set_hiword(&a,ha); /* a <- |a| */
|
|
|
90
|
+ set_hiword(&b, hb); /* b <- |b| */
|
|
|
91
|
+
|
|
|
92
|
+ if((ha-hb)>0x3c00000) {return a+b;} /* x/y > 2**60 */
|
|
|
93
|
+ k=0;
|
|
|
94
|
+ if(ha > 0x5f300000) { /* a>2**500 */
|
|
|
95
|
+ if(ha >= 0x7ff00000) { /* Inf or NaN */
|
|
|
96
|
+ w = a+b; /* for sNaN */
|
|
|
97
|
+ if(((ha&0xfffff)|get_loword(a))==0) w = a;
|
|
|
98
|
+ if(((hb^0x7ff00000)|get_loword(b))==0) w = b;
|
|
|
99
|
+ return w;
|
|
|
100
|
+ }
|
|
|
101
|
+ /* scale a and b by 2**-600 */
|
|
|
102
|
+ ha -= 0x25800000; hb -= 0x25800000; k += 600;
|
|
|
103
|
+ set_hiword(&a, ha);
|
|
|
104
|
+ set_hiword(&b, hb);
|
|
|
105
|
+ }
|
|
|
106
|
+ if(hb < 0x20b00000) { /* b < 2**-500 */
|
|
|
107
|
+ if(hb <= 0x000fffff) { /* subnormal b or 0 */
|
|
|
108
|
+ if((hb|(get_loword(b)))==0) return a;
|
|
|
109
|
+ t1=0;
|
|
|
110
|
+ set_hiword(&t1, 0x7fd00000); /* t1=2^1022 */
|
|
|
111
|
+ b *= t1;
|
|
|
112
|
+ a *= t1;
|
|
|
113
|
+ k -= 1022;
|
|
|
114
|
+ } else { /* scale a and b by 2^600 */
|
|
|
115
|
+ ha += 0x25800000; /* a *= 2^600 */
|
|
|
116
|
+ hb += 0x25800000; /* b *= 2^600 */
|
|
|
117
|
+ k -= 600;
|
|
|
118
|
+ set_hiword(&a, ha);
|
|
|
119
|
+ set_hiword(&b, hb);
|
|
|
120
|
+ }
|
|
|
121
|
+ }
|
|
|
122
|
+ /* medium size a and b */
|
|
|
123
|
+ w = a-b;
|
|
|
124
|
+ if (w>b) {
|
|
|
125
|
+ t1 = 0;
|
|
|
126
|
+ set_hiword(&t1, ha);
|
|
|
127
|
+ t2 = a-t1;
|
|
|
128
|
+ w = sqrt(t1*t1-(b*(-b)-t2*(a+t1)));
|
|
|
129
|
+ } else {
|
|
|
130
|
+ a = a+a;
|
|
|
131
|
+ y1 = 0;
|
|
|
132
|
+ set_hiword(&y1, hb);
|
|
|
133
|
+ y2 = b - y1;
|
|
|
134
|
+ t1 = 0;
|
|
|
135
|
+ set_hiword(&t1, ha+0x00100000);
|
|
|
136
|
+ t2 = a - t1;
|
|
|
137
|
+ w = sqrt(t1*y1-(w*(-w)-(t1*y2+t2*b)));
|
|
|
138
|
+ }
|
|
|
139
|
+ if(k!=0) {
|
|
|
140
|
+ t1 = 1.0;
|
|
|
141
|
+ set_hiword(&t1, get_hiword(t1) + (k<<20));
|
|
|
142
|
+ return t1*w;
|
|
|
143
|
+ } else return w;
|
|
|
144
|
+} |