Rev 828 | Rev 835 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 828 | Rev 829 | ||
---|---|---|---|
1 | /* |
1 | /* |
2 | * Copyright (C) 2005 Josef Cejka |
2 | * Copyright (C) 2005 Josef Cejka |
3 | * All rights reserved. |
3 | * All rights reserved. |
4 | * |
4 | * |
5 | * Redistribution and use in source and binary forms, with or without |
5 | * Redistribution and use in source and binary forms, with or without |
6 | * modification, are permitted provided that the following conditions |
6 | * modification, are permitted provided that the following conditions |
7 | * are met: |
7 | * are met: |
8 | * |
8 | * |
9 | * - Redistributions of source code must retain the above copyright |
9 | * - Redistributions of source code must retain the above copyright |
10 | * notice, this list of conditions and the following disclaimer. |
10 | * notice, this list of conditions and the following disclaimer. |
11 | * - Redistributions in binary form must reproduce the above copyright |
11 | * - Redistributions in binary form must reproduce the above copyright |
12 | * notice, this list of conditions and the following disclaimer in the |
12 | * notice, this list of conditions and the following disclaimer in the |
13 | * documentation and/or other materials provided with the distribution. |
13 | * documentation and/or other materials provided with the distribution. |
14 | * - The name of the author may not be used to endorse or promote products |
14 | * - The name of the author may not be used to endorse or promote products |
15 | * derived from this software without specific prior written permission. |
15 | * derived from this software without specific prior written permission. |
16 | * |
16 | * |
17 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR |
17 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR |
18 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES |
18 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES |
19 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. |
19 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. |
20 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, |
20 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, |
21 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT |
21 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT |
22 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |
22 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, |
23 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY |
23 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY |
24 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
24 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT |
25 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF |
25 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF |
26 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
26 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
27 | */ |
27 | */ |
28 | 28 | ||
29 | #include<sftypes.h> |
29 | #include<sftypes.h> |
30 | #include<add.h> |
30 | #include<add.h> |
31 | #include<div.h> |
31 | #include<div.h> |
32 | #include<comparison.h> |
32 | #include<comparison.h> |
33 | #include<mul.h> |
33 | #include<mul.h> |
- | 34 | #include<common.h> |
|
- | 35 | ||
34 | 36 | ||
35 | float32 divFloat32(float32 a, float32 b) |
37 | float32 divFloat32(float32 a, float32 b) |
36 | { |
38 | { |
37 | float32 result; |
39 | float32 result; |
38 | __s32 aexp, bexp, cexp; |
40 | __s32 aexp, bexp, cexp; |
39 | __u64 afrac, bfrac, cfrac; |
41 | __u64 afrac, bfrac, cfrac; |
40 | 42 | ||
41 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
43 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
42 | 44 | ||
43 | if (isFloat32NaN(a)) { |
45 | if (isFloat32NaN(a)) { |
44 | if (isFloat32SigNaN(a)) { |
46 | if (isFloat32SigNaN(a)) { |
45 | /*FIXME: SigNaN*/ |
47 | /*FIXME: SigNaN*/ |
46 | } |
48 | } |
47 | /*NaN*/ |
49 | /*NaN*/ |
48 | return a; |
50 | return a; |
49 | } |
51 | } |
50 | 52 | ||
51 | if (isFloat32NaN(b)) { |
53 | if (isFloat32NaN(b)) { |
52 | if (isFloat32SigNaN(b)) { |
54 | if (isFloat32SigNaN(b)) { |
53 | /*FIXME: SigNaN*/ |
55 | /*FIXME: SigNaN*/ |
54 | } |
56 | } |
55 | /*NaN*/ |
57 | /*NaN*/ |
56 | return b; |
58 | return b; |
57 | } |
59 | } |
58 | 60 | ||
59 | if (isFloat32Infinity(a)) { |
61 | if (isFloat32Infinity(a)) { |
60 | if (isFloat32Infinity(b)) { |
62 | if (isFloat32Infinity(b)) { |
61 | /*FIXME: inf / inf */ |
63 | /*FIXME: inf / inf */ |
62 | result.binary = FLOAT32_NAN; |
64 | result.binary = FLOAT32_NAN; |
63 | return result; |
65 | return result; |
64 | } |
66 | } |
65 | /* inf / num */ |
67 | /* inf / num */ |
66 | result.parts.exp = a.parts.exp; |
68 | result.parts.exp = a.parts.exp; |
67 | result.parts.fraction = a.parts.fraction; |
69 | result.parts.fraction = a.parts.fraction; |
68 | return result; |
70 | return result; |
69 | } |
71 | } |
70 | 72 | ||
71 | if (isFloat32Infinity(b)) { |
73 | if (isFloat32Infinity(b)) { |
72 | if (isFloat32Zero(a)) { |
74 | if (isFloat32Zero(a)) { |
73 | /* FIXME 0 / inf */ |
75 | /* FIXME 0 / inf */ |
74 | result.parts.exp = 0; |
76 | result.parts.exp = 0; |
75 | result.parts.fraction = 0; |
77 | result.parts.fraction = 0; |
76 | return result; |
78 | return result; |
77 | } |
79 | } |
78 | /* FIXME: num / inf*/ |
80 | /* FIXME: num / inf*/ |
79 | result.parts.exp = 0; |
81 | result.parts.exp = 0; |
80 | result.parts.fraction = 0; |
82 | result.parts.fraction = 0; |
81 | return result; |
83 | return result; |
82 | } |
84 | } |
83 | 85 | ||
84 | if (isFloat32Zero(b)) { |
86 | if (isFloat32Zero(b)) { |
85 | if (isFloat32Zero(a)) { |
87 | if (isFloat32Zero(a)) { |
86 | /*FIXME: 0 / 0*/ |
88 | /*FIXME: 0 / 0*/ |
87 | result.binary = FLOAT32_NAN; |
89 | result.binary = FLOAT32_NAN; |
88 | return result; |
90 | return result; |
89 | } |
91 | } |
90 | /* FIXME: division by zero */ |
92 | /* FIXME: division by zero */ |
91 | result.parts.exp = 0; |
93 | result.parts.exp = 0; |
92 | result.parts.fraction = 0; |
94 | result.parts.fraction = 0; |
93 | return result; |
95 | return result; |
94 | } |
96 | } |
95 | 97 | ||
96 | 98 | ||
97 | afrac = a.parts.fraction; |
99 | afrac = a.parts.fraction; |
98 | aexp = a.parts.exp; |
100 | aexp = a.parts.exp; |
99 | bfrac = b.parts.fraction; |
101 | bfrac = b.parts.fraction; |
100 | bexp = b.parts.exp; |
102 | bexp = b.parts.exp; |
101 | 103 | ||
102 | /* denormalized numbers */ |
104 | /* denormalized numbers */ |
103 | if (aexp == 0) { |
105 | if (aexp == 0) { |
104 | if (afrac == 0) { |
106 | if (afrac == 0) { |
105 | result.parts.exp = 0; |
107 | result.parts.exp = 0; |
106 | result.parts.fraction = 0; |
108 | result.parts.fraction = 0; |
107 | return result; |
109 | return result; |
108 | } |
110 | } |
109 | /* normalize it*/ |
111 | /* normalize it*/ |
110 | 112 | ||
111 | afrac <<= 1; |
113 | afrac <<= 1; |
112 | /* afrac is nonzero => it must stop */ |
114 | /* afrac is nonzero => it must stop */ |
113 | while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) { |
115 | while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) { |
114 | afrac <<= 1; |
116 | afrac <<= 1; |
115 | aexp--; |
117 | aexp--; |
116 | } |
118 | } |
117 | } |
119 | } |
118 | 120 | ||
119 | if (bexp == 0) { |
121 | if (bexp == 0) { |
120 | bfrac <<= 1; |
122 | bfrac <<= 1; |
121 | /* bfrac is nonzero => it must stop */ |
123 | /* bfrac is nonzero => it must stop */ |
122 | while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) { |
124 | while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) { |
123 | bfrac <<= 1; |
125 | bfrac <<= 1; |
124 | bexp--; |
126 | bexp--; |
125 | } |
127 | } |
126 | } |
128 | } |
127 | 129 | ||
128 | afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 ); |
130 | afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 ); |
129 | bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE ); |
131 | bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE ); |
130 | 132 | ||
131 | if ( bfrac <= (afrac << 1) ) { |
133 | if ( bfrac <= (afrac << 1) ) { |
132 | afrac >>= 1; |
134 | afrac >>= 1; |
133 | aexp++; |
135 | aexp++; |
134 | } |
136 | } |
135 | 137 | ||
136 | cexp = aexp - bexp + FLOAT32_BIAS - 2; |
138 | cexp = aexp - bexp + FLOAT32_BIAS - 2; |
137 | 139 | ||
138 | cfrac = (afrac << 32) / bfrac; |
140 | cfrac = (afrac << 32) / bfrac; |
139 | if (( cfrac & 0x3F ) == 0) { |
141 | if (( cfrac & 0x3F ) == 0) { |
140 | cfrac |= ( bfrac * cfrac != afrac << 32 ); |
142 | cfrac |= ( bfrac * cfrac != afrac << 32 ); |
141 | } |
143 | } |
142 | 144 | ||
143 | /* pack and round */ |
145 | /* pack and round */ |
144 | 146 | ||
145 | /* find first nonzero digit and shift result and detect possibly underflow */ |
147 | /* find first nonzero digit and shift result and detect possibly underflow */ |
146 | while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) { |
148 | while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) { |
147 | cexp--; |
149 | cexp--; |
148 | cfrac <<= 1; |
150 | cfrac <<= 1; |
149 | /* TODO: fix underflow */ |
151 | /* TODO: fix underflow */ |
150 | }; |
152 | }; |
151 | 153 | ||
152 | cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/ |
154 | cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/ |
153 | 155 | ||
154 | if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) { |
156 | if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) { |
155 | ++cexp; |
157 | ++cexp; |
156 | cfrac >>= 1; |
158 | cfrac >>= 1; |
157 | } |
159 | } |
158 | 160 | ||
159 | /* check overflow */ |
161 | /* check overflow */ |
160 | if (cexp >= FLOAT32_MAX_EXPONENT ) { |
162 | if (cexp >= FLOAT32_MAX_EXPONENT ) { |
161 | /* FIXME: overflow, return infinity */ |
163 | /* FIXME: overflow, return infinity */ |
162 | result.parts.exp = FLOAT32_MAX_EXPONENT; |
164 | result.parts.exp = FLOAT32_MAX_EXPONENT; |
163 | result.parts.fraction = 0; |
165 | result.parts.fraction = 0; |
164 | return result; |
166 | return result; |
165 | } |
167 | } |
166 | 168 | ||
167 | if (cexp < 0) { |
169 | if (cexp < 0) { |
168 | /* FIXME: underflow */ |
170 | /* FIXME: underflow */ |
169 | result.parts.exp = 0; |
171 | result.parts.exp = 0; |
170 | if ((cexp + FLOAT32_FRACTION_SIZE) < 0) { |
172 | if ((cexp + FLOAT32_FRACTION_SIZE) < 0) { |
171 | result.parts.fraction = 0; |
173 | result.parts.fraction = 0; |
172 | return result; |
174 | return result; |
173 | } |
175 | } |
174 | cfrac >>= 1; |
176 | cfrac >>= 1; |
175 | while (cexp < 0) { |
177 | while (cexp < 0) { |
176 | cexp ++; |
178 | cexp ++; |
177 | cfrac >>= 1; |
179 | cfrac >>= 1; |
178 | } |
180 | } |
179 | 181 | ||
180 | } else { |
182 | } else { |
181 | result.parts.exp = (__u32)cexp; |
183 | result.parts.exp = (__u32)cexp; |
182 | } |
184 | } |
183 | 185 | ||
184 | result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); |
186 | result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); |
185 | 187 | ||
186 | return result; |
188 | return result; |
187 | } |
189 | } |
188 | 190 | ||
189 | float64 divFloat64(float64 a, float64 b) |
191 | float64 divFloat64(float64 a, float64 b) |
190 | { |
192 | { |
191 | float64 result; |
193 | float64 result; |
192 | __s32 aexp, bexp, cexp; |
194 | __s32 aexp, bexp, cexp; |
193 | __u64 afrac, bfrac, cfrac; |
195 | __u64 afrac, bfrac, cfrac; |
194 | __u64 remlo, remhi; |
196 | __u64 remlo, remhi; |
195 | 197 | ||
196 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
198 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
197 | 199 | ||
198 | if (isFloat64NaN(a)) { |
200 | if (isFloat64NaN(a)) { |
199 | if (isFloat64SigNaN(a)) { |
201 | if (isFloat64SigNaN(a)) { |
200 | /*FIXME: SigNaN*/ |
202 | /*FIXME: SigNaN*/ |
201 | } |
203 | } |
202 | /*NaN*/ |
204 | /*NaN*/ |
203 | return a; |
205 | return a; |
204 | } |
206 | } |
205 | 207 | ||
206 | if (isFloat64NaN(b)) { |
208 | if (isFloat64NaN(b)) { |
207 | if (isFloat64SigNaN(b)) { |
209 | if (isFloat64SigNaN(b)) { |
208 | /*FIXME: SigNaN*/ |
210 | /*FIXME: SigNaN*/ |
209 | } |
211 | } |
210 | /*NaN*/ |
212 | /*NaN*/ |
211 | return b; |
213 | return b; |
212 | } |
214 | } |
213 | 215 | ||
214 | if (isFloat64Infinity(a)) { |
216 | if (isFloat64Infinity(a)) { |
215 | if (isFloat64Infinity(b)) { |
217 | if (isFloat64Infinity(b)) { |
216 | /*FIXME: inf / inf */ |
218 | /*FIXME: inf / inf */ |
217 | result.binary = FLOAT64_NAN; |
219 | result.binary = FLOAT64_NAN; |
218 | return result; |
220 | return result; |
219 | } |
221 | } |
220 | /* inf / num */ |
222 | /* inf / num */ |
221 | result.parts.exp = a.parts.exp; |
223 | result.parts.exp = a.parts.exp; |
222 | result.parts.fraction = a.parts.fraction; |
224 | result.parts.fraction = a.parts.fraction; |
223 | return result; |
225 | return result; |
224 | } |
226 | } |
225 | 227 | ||
226 | if (isFloat64Infinity(b)) { |
228 | if (isFloat64Infinity(b)) { |
227 | if (isFloat64Zero(a)) { |
229 | if (isFloat64Zero(a)) { |
228 | /* FIXME 0 / inf */ |
230 | /* FIXME 0 / inf */ |
229 | result.parts.exp = 0; |
231 | result.parts.exp = 0; |
230 | result.parts.fraction = 0; |
232 | result.parts.fraction = 0; |
231 | return result; |
233 | return result; |
232 | } |
234 | } |
233 | /* FIXME: num / inf*/ |
235 | /* FIXME: num / inf*/ |
234 | result.parts.exp = 0; |
236 | result.parts.exp = 0; |
235 | result.parts.fraction = 0; |
237 | result.parts.fraction = 0; |
236 | return result; |
238 | return result; |
237 | } |
239 | } |
238 | 240 | ||
239 | if (isFloat64Zero(b)) { |
241 | if (isFloat64Zero(b)) { |
240 | if (isFloat64Zero(a)) { |
242 | if (isFloat64Zero(a)) { |
241 | /*FIXME: 0 / 0*/ |
243 | /*FIXME: 0 / 0*/ |
242 | result.binary = FLOAT64_NAN; |
244 | result.binary = FLOAT64_NAN; |
243 | return result; |
245 | return result; |
244 | } |
246 | } |
245 | /* FIXME: division by zero */ |
247 | /* FIXME: division by zero */ |
246 | result.parts.exp = 0; |
248 | result.parts.exp = 0; |
247 | result.parts.fraction = 0; |
249 | result.parts.fraction = 0; |
248 | return result; |
250 | return result; |
249 | } |
251 | } |
250 | 252 | ||
251 | 253 | ||
252 | afrac = a.parts.fraction; |
254 | afrac = a.parts.fraction; |
253 | aexp = a.parts.exp; |
255 | aexp = a.parts.exp; |
254 | bfrac = b.parts.fraction; |
256 | bfrac = b.parts.fraction; |
255 | bexp = b.parts.exp; |
257 | bexp = b.parts.exp; |
256 | 258 | ||
257 | /* denormalized numbers */ |
259 | /* denormalized numbers */ |
258 | if (aexp == 0) { |
260 | if (aexp == 0) { |
259 | if (afrac == 0) { |
261 | if (afrac == 0) { |
260 | result.parts.exp = 0; |
262 | result.parts.exp = 0; |
261 | result.parts.fraction = 0; |
263 | result.parts.fraction = 0; |
262 | return result; |
264 | return result; |
263 | } |
265 | } |
264 | /* normalize it*/ |
266 | /* normalize it*/ |
265 | 267 | ||
266 | afrac <<= 1; |
268 | afrac <<= 1; |
267 | /* afrac is nonzero => it must stop */ |
269 | /* afrac is nonzero => it must stop */ |
268 | while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
270 | while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
269 | afrac <<= 1; |
271 | afrac <<= 1; |
270 | aexp--; |
272 | aexp--; |
271 | } |
273 | } |
272 | } |
274 | } |
273 | 275 | ||
274 | if (bexp == 0) { |
276 | if (bexp == 0) { |
275 | bfrac <<= 1; |
277 | bfrac <<= 1; |
276 | /* bfrac is nonzero => it must stop */ |
278 | /* bfrac is nonzero => it must stop */ |
277 | while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
279 | while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
278 | bfrac <<= 1; |
280 | bfrac <<= 1; |
279 | bexp--; |
281 | bexp--; |
280 | } |
282 | } |
281 | } |
283 | } |
282 | 284 | ||
283 | afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 ); |
285 | afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 ); |
284 | bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1); |
286 | bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1); |
285 | 287 | ||
286 | if ( bfrac <= (afrac << 1) ) { |
288 | if ( bfrac <= (afrac << 1) ) { |
287 | afrac >>= 1; |
289 | afrac >>= 1; |
288 | aexp++; |
290 | aexp++; |
289 | } |
291 | } |
290 | 292 | ||
291 | cexp = aexp - bexp + FLOAT64_BIAS - 2; |
293 | cexp = aexp - bexp + FLOAT64_BIAS - 2; |
292 | 294 | ||
293 | cfrac = divFloat64estim(afrac, bfrac); |
295 | cfrac = divFloat64estim(afrac, bfrac); |
294 | 296 | ||
295 | if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */ |
297 | if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */ |
296 | mul64integers( bfrac, cfrac, &remlo, &remhi); |
298 | mul64integers( bfrac, cfrac, &remlo, &remhi); |
297 | /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ |
299 | /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ |
298 | remhi = afrac - remhi - ( remlo > 0); |
300 | remhi = afrac - remhi - ( remlo > 0); |
299 | remlo = - remlo; |
301 | remlo = - remlo; |
300 | 302 | ||
301 | while ((__s64) remhi < 0) { |
303 | while ((__s64) remhi < 0) { |
302 | cfrac--; |
304 | cfrac--; |
303 | remlo += bfrac; |
305 | remlo += bfrac; |
304 | remhi += ( remlo < bfrac ); |
306 | remhi += ( remlo < bfrac ); |
305 | } |
307 | } |
306 | cfrac |= ( remlo != 0 ); |
308 | cfrac |= ( remlo != 0 ); |
307 | } |
309 | } |
308 | 310 | ||
309 | /* pack and round */ |
311 | /* round and shift */ |
310 | - | ||
311 | /* find first nonzero digit and shift result and detect possibly underflow */ |
- | |
312 | while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ) )))) { |
- | |
313 | cexp--; |
- | |
314 | cfrac <<= 1; |
- | |
315 | /* TODO: fix underflow */ |
- | |
316 | }; |
- | |
317 | - | ||
318 | - | ||
319 | cfrac >>= 1; |
- | |
320 | ++cexp; |
- | |
321 | cfrac += (0x1 << (64 - FLOAT64_FRACTION_SIZE - 3)); |
- | |
322 | - | ||
323 | if (cfrac & (FLOAT64_HIDDEN_BIT_MASK << (64 - FLOAT64_FRACTION_SIZE - 1 ))) { |
- | |
324 | ++cexp; |
- | |
325 | cfrac >>= 1; |
- | |
326 | } |
- | |
327 | - | ||
328 | /* check overflow */ |
- | |
329 | if (cexp >= FLOAT64_MAX_EXPONENT ) { |
- | |
330 | /* FIXME: overflow, return infinity */ |
- | |
331 | result.parts.exp = FLOAT64_MAX_EXPONENT; |
312 | result = finishFloat64(cexp, cfrac, result.parts.sign); |
332 | result.parts.fraction = 0; |
- | |
333 | return result; |
313 | return result; |
334 | } |
- | |
335 | 314 | ||
336 | if (cexp < 0) { |
- | |
337 | /* FIXME: underflow */ |
- | |
338 | result.parts.exp = 0; |
- | |
339 | if ((cexp + FLOAT64_FRACTION_SIZE) < 0) { |
- | |
340 | result.parts.fraction = 0; |
- | |
341 | return result; |
- | |
342 | } |
- | |
343 | cfrac >>= 1; |
- | |
344 | while (cexp < 0) { |
- | |
345 | cexp ++; |
- | |
346 | cfrac >>= 1; |
- | |
347 | } |
- | |
348 | return result; |
- | |
349 | - | ||
350 | } else { |
- | |
351 | cexp ++; /*normalized*/ |
- | |
352 | result.parts.exp = (__u32)cexp; |
- | |
353 | } |
- | |
354 | - | ||
355 | result.parts.fraction = ((cfrac >>(64 - FLOAT64_FRACTION_SIZE - 2 ) ) & (~FLOAT64_HIDDEN_BIT_MASK)); |
- | |
356 | - | ||
357 | return result; |
- | |
358 | } |
315 | } |
359 | 316 | ||
360 | __u64 divFloat64estim(__u64 a, __u64 b) |
317 | __u64 divFloat64estim(__u64 a, __u64 b) |
361 | { |
318 | { |
362 | __u64 bhi; |
319 | __u64 bhi; |
363 | __u64 remhi, remlo; |
320 | __u64 remhi, remlo; |
364 | __u64 result; |
321 | __u64 result; |
365 | 322 | ||
366 | if ( b <= a ) { |
323 | if ( b <= a ) { |
367 | return 0xFFFFFFFFFFFFFFFFull; |
324 | return 0xFFFFFFFFFFFFFFFFull; |
368 | } |
325 | } |
369 | 326 | ||
370 | bhi = b >> 32; |
327 | bhi = b >> 32; |
371 | result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32; |
328 | result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32; |
372 | mul64integers(b, result, &remlo, &remhi); |
329 | mul64integers(b, result, &remlo, &remhi); |
373 | 330 | ||
374 | remhi = a - remhi - (remlo > 0); |
331 | remhi = a - remhi - (remlo > 0); |
375 | remlo = - remlo; |
332 | remlo = - remlo; |
376 | 333 | ||
377 | b <<= 32; |
334 | b <<= 32; |
378 | while ( (__s64) remhi < 0 ) { |
335 | while ( (__s64) remhi < 0 ) { |
379 | result -= 0x1ll << 32; |
336 | result -= 0x1ll << 32; |
380 | remlo += b; |
337 | remlo += b; |
381 | remhi += bhi + ( remlo < b ); |
338 | remhi += bhi + ( remlo < b ); |
382 | } |
339 | } |
383 | remhi = (remhi << 32) | (remlo >> 32); |
340 | remhi = (remhi << 32) | (remlo >> 32); |
384 | if (( bhi << 32) <= remhi) { |
341 | if (( bhi << 32) <= remhi) { |
385 | result |= 0xFFFFFFFF; |
342 | result |= 0xFFFFFFFF; |
386 | } else { |
343 | } else { |
387 | result |= remhi / bhi; |
344 | result |= remhi / bhi; |
388 | } |
345 | } |
389 | 346 | ||
390 | 347 | ||
391 | return result; |
348 | return result; |
392 | } |
349 | } |
393 | 350 | ||
394 | 351 |