Commit | Line | Data |
---|---|---|
1da177e4 LT |
1 | /* IEEE754 floating point arithmetic |
2 | * double precision square root | |
3 | */ | |
4 | /* | |
5 | * MIPS floating point support | |
6 | * Copyright (C) 1994-2000 Algorithmics Ltd. | |
1da177e4 | 7 | * |
1da177e4 LT |
8 | * This program is free software; you can distribute it and/or modify it |
9 | * under the terms of the GNU General Public License (Version 2) as | |
10 | * published by the Free Software Foundation. | |
11 | * | |
12 | * This program is distributed in the hope it will be useful, but WITHOUT | |
13 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or | |
14 | * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License | |
15 | * for more details. | |
16 | * | |
17 | * You should have received a copy of the GNU General Public License along | |
18 | * with this program; if not, write to the Free Software Foundation, Inc., | |
3f7cac41 | 19 | * 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. |
1da177e4 LT |
20 | */ |
21 | ||
1da177e4 LT |
22 | #include "ieee754dp.h" |
23 | ||
24 | static const unsigned table[] = { | |
25 | 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, | |
26 | 29598, 36145, 43202, 50740, 58733, 67158, 75992, | |
27 | 85215, 83599, 71378, 60428, 50647, 41945, 34246, | |
28 | 27478, 21581, 16499, 12183, 8588, 5674, 3403, | |
29 | 1742, 661, 130 | |
30 | }; | |
31 | ||
2209bcb1 | 32 | union ieee754dp ieee754dp_sqrt(union ieee754dp x) |
1da177e4 | 33 | { |
cd21dfcf | 34 | struct _ieee754_csr oldcsr; |
2209bcb1 | 35 | union ieee754dp y, z, t; |
1da177e4 LT |
36 | unsigned scalx, yh; |
37 | COMPXDP; | |
38 | ||
39 | EXPLODEXDP; | |
9e8bad1f | 40 | ieee754_clearcx(); |
1da177e4 LT |
41 | FLUSHXDP; |
42 | ||
43 | /* x == INF or NAN? */ | |
44 | switch (xc) { | |
d5afa7e9 MR |
45 | case IEEE754_CLASS_SNAN: |
46 | return ieee754dp_nanxcpt(x); | |
47 | ||
1da177e4 LT |
48 | case IEEE754_CLASS_QNAN: |
49 | /* sqrt(Nan) = Nan */ | |
539bfb57 | 50 | return x; |
3f7cac41 | 51 | |
1da177e4 LT |
52 | case IEEE754_CLASS_ZERO: |
53 | /* sqrt(0) = 0 */ | |
54 | return x; | |
3f7cac41 | 55 | |
1da177e4 LT |
56 | case IEEE754_CLASS_INF: |
57 | if (xs) { | |
58 | /* sqrt(-Inf) = Nan */ | |
9e8bad1f | 59 | ieee754_setcx(IEEE754_INVALID_OPERATION); |
539bfb57 | 60 | return ieee754dp_indef(); |
1da177e4 LT |
61 | } |
62 | /* sqrt(+Inf) = Inf */ | |
63 | return x; | |
3f7cac41 | 64 | |
1da177e4 LT |
65 | case IEEE754_CLASS_DNORM: |
66 | DPDNORMX; | |
67 | /* fall through */ | |
3f7cac41 | 68 | |
1da177e4 LT |
69 | case IEEE754_CLASS_NORM: |
70 | if (xs) { | |
71 | /* sqrt(-x) = Nan */ | |
9e8bad1f | 72 | ieee754_setcx(IEEE754_INVALID_OPERATION); |
539bfb57 | 73 | return ieee754dp_indef(); |
1da177e4 LT |
74 | } |
75 | break; | |
76 | } | |
77 | ||
78 | /* save old csr; switch off INX enable & flag; set RN rounding */ | |
79 | oldcsr = ieee754_csr; | |
80 | ieee754_csr.mx &= ~IEEE754_INEXACT; | |
81 | ieee754_csr.sx &= ~IEEE754_INEXACT; | |
56a64733 | 82 | ieee754_csr.rm = FPU_CSR_RN; |
1da177e4 LT |
83 | |
84 | /* adjust exponent to prevent overflow */ | |
85 | scalx = 0; | |
86 | if (xe > 512) { /* x > 2**-512? */ | |
87 | xe -= 512; /* x = x / 2**512 */ | |
88 | scalx += 256; | |
70342287 | 89 | } else if (xe < -512) { /* x < 2**-512? */ |
1da177e4 LT |
90 | xe += 512; /* x = x * 2**512 */ |
91 | scalx -= 256; | |
92 | } | |
93 | ||
94 | y = x = builddp(0, xe + DP_EBIAS, xm & ~DP_HIDDEN_BIT); | |
95 | ||
96 | /* magic initial approximation to almost 8 sig. bits */ | |
97 | yh = y.bits >> 32; | |
98 | yh = (yh >> 1) + 0x1ff80000; | |
99 | yh = yh - table[(yh >> 15) & 31]; | |
100 | y.bits = ((u64) yh << 32) | (y.bits & 0xffffffff); | |
101 | ||
102 | /* Heron's rule once with correction to improve to ~18 sig. bits */ | |
103 | /* t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0; */ | |
104 | t = ieee754dp_div(x, y); | |
105 | y = ieee754dp_add(y, t); | |
106 | y.bits -= 0x0010000600000000LL; | |
107 | y.bits &= 0xffffffff00000000LL; | |
108 | ||
109 | /* triple to almost 56 sig. bits: y ~= sqrt(x) to within 1 ulp */ | |
70342287 | 110 | /* t=y*y; z=t; pt[n0]+=0x00100000; t+=z; z=(x-z)*y; */ |
1da177e4 | 111 | z = t = ieee754dp_mul(y, y); |
49548b09 | 112 | t.bexp += 0x001; |
1da177e4 LT |
113 | t = ieee754dp_add(t, z); |
114 | z = ieee754dp_mul(ieee754dp_sub(x, z), y); | |
115 | ||
70342287 | 116 | /* t=z/(t+x) ; pt[n0]+=0x00100000; y+=t; */ |
1da177e4 | 117 | t = ieee754dp_div(z, ieee754dp_add(t, x)); |
49548b09 | 118 | t.bexp += 0x001; |
1da177e4 LT |
119 | y = ieee754dp_add(y, t); |
120 | ||
121 | /* twiddle last bit to force y correctly rounded */ | |
122 | ||
123 | /* set RZ, clear INEX flag */ | |
56a64733 | 124 | ieee754_csr.rm = FPU_CSR_RZ; |
1da177e4 LT |
125 | ieee754_csr.sx &= ~IEEE754_INEXACT; |
126 | ||
127 | /* t=x/y; ...chopped quotient, possibly inexact */ | |
128 | t = ieee754dp_div(x, y); | |
129 | ||
130 | if (ieee754_csr.sx & IEEE754_INEXACT || t.bits != y.bits) { | |
131 | ||
132 | if (!(ieee754_csr.sx & IEEE754_INEXACT)) | |
133 | /* t = t-ulp */ | |
134 | t.bits -= 1; | |
135 | ||
136 | /* add inexact to result status */ | |
137 | oldcsr.cx |= IEEE754_INEXACT; | |
138 | oldcsr.sx |= IEEE754_INEXACT; | |
139 | ||
140 | switch (oldcsr.rm) { | |
56a64733 | 141 | case FPU_CSR_RU: |
1da177e4 LT |
142 | y.bits += 1; |
143 | /* drop through */ | |
56a64733 | 144 | case FPU_CSR_RN: |
1da177e4 LT |
145 | t.bits += 1; |
146 | break; | |
147 | } | |
148 | ||
149 | /* y=y+t; ...chopped sum */ | |
150 | y = ieee754dp_add(y, t); | |
151 | ||
152 | /* adjust scalx for correctly rounded sqrt(x) */ | |
153 | scalx -= 1; | |
154 | } | |
155 | ||
156 | /* py[n0]=py[n0]+scalx; ...scale back y */ | |
49548b09 | 157 | y.bexp += scalx; |
1da177e4 LT |
158 | |
159 | /* restore rounding mode, possibly set inexact */ | |
160 | ieee754_csr = oldcsr; | |
161 | ||
162 | return y; | |
163 | } |