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