Rev 1031 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 1031 | Rev 1657 | ||
---|---|---|---|
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 | /** @addtogroup softfloat |
|
- | 30 | * @{ |
|
- | 31 | */ |
|
- | 32 | /** @file |
|
- | 33 | */ |
|
- | 34 | ||
29 | #include<sftypes.h> |
35 | #include<sftypes.h> |
30 | #include<add.h> |
36 | #include<add.h> |
31 | #include<div.h> |
37 | #include<div.h> |
32 | #include<comparison.h> |
38 | #include<comparison.h> |
33 | #include<mul.h> |
39 | #include<mul.h> |
34 | #include<common.h> |
40 | #include<common.h> |
35 | 41 | ||
36 | 42 | ||
37 | float32 divFloat32(float32 a, float32 b) |
43 | float32 divFloat32(float32 a, float32 b) |
38 | { |
44 | { |
39 | float32 result; |
45 | float32 result; |
40 | int32_t aexp, bexp, cexp; |
46 | int32_t aexp, bexp, cexp; |
41 | uint64_t afrac, bfrac, cfrac; |
47 | uint64_t afrac, bfrac, cfrac; |
42 | 48 | ||
43 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
49 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
44 | 50 | ||
45 | if (isFloat32NaN(a)) { |
51 | if (isFloat32NaN(a)) { |
46 | if (isFloat32SigNaN(a)) { |
52 | if (isFloat32SigNaN(a)) { |
47 | /*FIXME: SigNaN*/ |
53 | /*FIXME: SigNaN*/ |
48 | } |
54 | } |
49 | /*NaN*/ |
55 | /*NaN*/ |
50 | return a; |
56 | return a; |
51 | } |
57 | } |
52 | 58 | ||
53 | if (isFloat32NaN(b)) { |
59 | if (isFloat32NaN(b)) { |
54 | if (isFloat32SigNaN(b)) { |
60 | if (isFloat32SigNaN(b)) { |
55 | /*FIXME: SigNaN*/ |
61 | /*FIXME: SigNaN*/ |
56 | } |
62 | } |
57 | /*NaN*/ |
63 | /*NaN*/ |
58 | return b; |
64 | return b; |
59 | } |
65 | } |
60 | 66 | ||
61 | if (isFloat32Infinity(a)) { |
67 | if (isFloat32Infinity(a)) { |
62 | if (isFloat32Infinity(b)) { |
68 | if (isFloat32Infinity(b)) { |
63 | /*FIXME: inf / inf */ |
69 | /*FIXME: inf / inf */ |
64 | result.binary = FLOAT32_NAN; |
70 | result.binary = FLOAT32_NAN; |
65 | return result; |
71 | return result; |
66 | } |
72 | } |
67 | /* inf / num */ |
73 | /* inf / num */ |
68 | result.parts.exp = a.parts.exp; |
74 | result.parts.exp = a.parts.exp; |
69 | result.parts.fraction = a.parts.fraction; |
75 | result.parts.fraction = a.parts.fraction; |
70 | return result; |
76 | return result; |
71 | } |
77 | } |
72 | 78 | ||
73 | if (isFloat32Infinity(b)) { |
79 | if (isFloat32Infinity(b)) { |
74 | if (isFloat32Zero(a)) { |
80 | if (isFloat32Zero(a)) { |
75 | /* FIXME 0 / inf */ |
81 | /* FIXME 0 / inf */ |
76 | result.parts.exp = 0; |
82 | result.parts.exp = 0; |
77 | result.parts.fraction = 0; |
83 | result.parts.fraction = 0; |
78 | return result; |
84 | return result; |
79 | } |
85 | } |
80 | /* FIXME: num / inf*/ |
86 | /* FIXME: num / inf*/ |
81 | result.parts.exp = 0; |
87 | result.parts.exp = 0; |
82 | result.parts.fraction = 0; |
88 | result.parts.fraction = 0; |
83 | return result; |
89 | return result; |
84 | } |
90 | } |
85 | 91 | ||
86 | if (isFloat32Zero(b)) { |
92 | if (isFloat32Zero(b)) { |
87 | if (isFloat32Zero(a)) { |
93 | if (isFloat32Zero(a)) { |
88 | /*FIXME: 0 / 0*/ |
94 | /*FIXME: 0 / 0*/ |
89 | result.binary = FLOAT32_NAN; |
95 | result.binary = FLOAT32_NAN; |
90 | return result; |
96 | return result; |
91 | } |
97 | } |
92 | /* FIXME: division by zero */ |
98 | /* FIXME: division by zero */ |
93 | result.parts.exp = 0; |
99 | result.parts.exp = 0; |
94 | result.parts.fraction = 0; |
100 | result.parts.fraction = 0; |
95 | return result; |
101 | return result; |
96 | } |
102 | } |
97 | 103 | ||
98 | 104 | ||
99 | afrac = a.parts.fraction; |
105 | afrac = a.parts.fraction; |
100 | aexp = a.parts.exp; |
106 | aexp = a.parts.exp; |
101 | bfrac = b.parts.fraction; |
107 | bfrac = b.parts.fraction; |
102 | bexp = b.parts.exp; |
108 | bexp = b.parts.exp; |
103 | 109 | ||
104 | /* denormalized numbers */ |
110 | /* denormalized numbers */ |
105 | if (aexp == 0) { |
111 | if (aexp == 0) { |
106 | if (afrac == 0) { |
112 | if (afrac == 0) { |
107 | result.parts.exp = 0; |
113 | result.parts.exp = 0; |
108 | result.parts.fraction = 0; |
114 | result.parts.fraction = 0; |
109 | return result; |
115 | return result; |
110 | } |
116 | } |
111 | /* normalize it*/ |
117 | /* normalize it*/ |
112 | 118 | ||
113 | afrac <<= 1; |
119 | afrac <<= 1; |
114 | /* afrac is nonzero => it must stop */ |
120 | /* afrac is nonzero => it must stop */ |
115 | while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) { |
121 | while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) { |
116 | afrac <<= 1; |
122 | afrac <<= 1; |
117 | aexp--; |
123 | aexp--; |
118 | } |
124 | } |
119 | } |
125 | } |
120 | 126 | ||
121 | if (bexp == 0) { |
127 | if (bexp == 0) { |
122 | bfrac <<= 1; |
128 | bfrac <<= 1; |
123 | /* bfrac is nonzero => it must stop */ |
129 | /* bfrac is nonzero => it must stop */ |
124 | while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) { |
130 | while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) { |
125 | bfrac <<= 1; |
131 | bfrac <<= 1; |
126 | bexp--; |
132 | bexp--; |
127 | } |
133 | } |
128 | } |
134 | } |
129 | 135 | ||
130 | afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 ); |
136 | afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 ); |
131 | bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE ); |
137 | bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE ); |
132 | 138 | ||
133 | if ( bfrac <= (afrac << 1) ) { |
139 | if ( bfrac <= (afrac << 1) ) { |
134 | afrac >>= 1; |
140 | afrac >>= 1; |
135 | aexp++; |
141 | aexp++; |
136 | } |
142 | } |
137 | 143 | ||
138 | cexp = aexp - bexp + FLOAT32_BIAS - 2; |
144 | cexp = aexp - bexp + FLOAT32_BIAS - 2; |
139 | 145 | ||
140 | cfrac = (afrac << 32) / bfrac; |
146 | cfrac = (afrac << 32) / bfrac; |
141 | if (( cfrac & 0x3F ) == 0) { |
147 | if (( cfrac & 0x3F ) == 0) { |
142 | cfrac |= ( bfrac * cfrac != afrac << 32 ); |
148 | cfrac |= ( bfrac * cfrac != afrac << 32 ); |
143 | } |
149 | } |
144 | 150 | ||
145 | /* pack and round */ |
151 | /* pack and round */ |
146 | 152 | ||
147 | /* find first nonzero digit and shift result and detect possibly underflow */ |
153 | /* find first nonzero digit and shift result and detect possibly underflow */ |
148 | while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) { |
154 | while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) { |
149 | cexp--; |
155 | cexp--; |
150 | cfrac <<= 1; |
156 | cfrac <<= 1; |
151 | /* TODO: fix underflow */ |
157 | /* TODO: fix underflow */ |
152 | }; |
158 | }; |
153 | 159 | ||
154 | cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/ |
160 | cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/ |
155 | 161 | ||
156 | if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) { |
162 | if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) { |
157 | ++cexp; |
163 | ++cexp; |
158 | cfrac >>= 1; |
164 | cfrac >>= 1; |
159 | } |
165 | } |
160 | 166 | ||
161 | /* check overflow */ |
167 | /* check overflow */ |
162 | if (cexp >= FLOAT32_MAX_EXPONENT ) { |
168 | if (cexp >= FLOAT32_MAX_EXPONENT ) { |
163 | /* FIXME: overflow, return infinity */ |
169 | /* FIXME: overflow, return infinity */ |
164 | result.parts.exp = FLOAT32_MAX_EXPONENT; |
170 | result.parts.exp = FLOAT32_MAX_EXPONENT; |
165 | result.parts.fraction = 0; |
171 | result.parts.fraction = 0; |
166 | return result; |
172 | return result; |
167 | } |
173 | } |
168 | 174 | ||
169 | if (cexp < 0) { |
175 | if (cexp < 0) { |
170 | /* FIXME: underflow */ |
176 | /* FIXME: underflow */ |
171 | result.parts.exp = 0; |
177 | result.parts.exp = 0; |
172 | if ((cexp + FLOAT32_FRACTION_SIZE) < 0) { |
178 | if ((cexp + FLOAT32_FRACTION_SIZE) < 0) { |
173 | result.parts.fraction = 0; |
179 | result.parts.fraction = 0; |
174 | return result; |
180 | return result; |
175 | } |
181 | } |
176 | cfrac >>= 1; |
182 | cfrac >>= 1; |
177 | while (cexp < 0) { |
183 | while (cexp < 0) { |
178 | cexp ++; |
184 | cexp ++; |
179 | cfrac >>= 1; |
185 | cfrac >>= 1; |
180 | } |
186 | } |
181 | 187 | ||
182 | } else { |
188 | } else { |
183 | result.parts.exp = (uint32_t)cexp; |
189 | result.parts.exp = (uint32_t)cexp; |
184 | } |
190 | } |
185 | 191 | ||
186 | result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); |
192 | result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); |
187 | 193 | ||
188 | return result; |
194 | return result; |
189 | } |
195 | } |
190 | 196 | ||
191 | float64 divFloat64(float64 a, float64 b) |
197 | float64 divFloat64(float64 a, float64 b) |
192 | { |
198 | { |
193 | float64 result; |
199 | float64 result; |
194 | int64_t aexp, bexp, cexp; |
200 | int64_t aexp, bexp, cexp; |
195 | uint64_t afrac, bfrac, cfrac; |
201 | uint64_t afrac, bfrac, cfrac; |
196 | uint64_t remlo, remhi; |
202 | uint64_t remlo, remhi; |
197 | 203 | ||
198 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
204 | result.parts.sign = a.parts.sign ^ b.parts.sign; |
199 | 205 | ||
200 | if (isFloat64NaN(a)) { |
206 | if (isFloat64NaN(a)) { |
201 | 207 | ||
202 | if (isFloat64SigNaN(b)) { |
208 | if (isFloat64SigNaN(b)) { |
203 | /*FIXME: SigNaN*/ |
209 | /*FIXME: SigNaN*/ |
204 | return b; |
210 | return b; |
205 | } |
211 | } |
206 | 212 | ||
207 | if (isFloat64SigNaN(a)) { |
213 | if (isFloat64SigNaN(a)) { |
208 | /*FIXME: SigNaN*/ |
214 | /*FIXME: SigNaN*/ |
209 | } |
215 | } |
210 | /*NaN*/ |
216 | /*NaN*/ |
211 | return a; |
217 | return a; |
212 | } |
218 | } |
213 | 219 | ||
214 | if (isFloat64NaN(b)) { |
220 | if (isFloat64NaN(b)) { |
215 | if (isFloat64SigNaN(b)) { |
221 | if (isFloat64SigNaN(b)) { |
216 | /*FIXME: SigNaN*/ |
222 | /*FIXME: SigNaN*/ |
217 | } |
223 | } |
218 | /*NaN*/ |
224 | /*NaN*/ |
219 | return b; |
225 | return b; |
220 | } |
226 | } |
221 | 227 | ||
222 | if (isFloat64Infinity(a)) { |
228 | if (isFloat64Infinity(a)) { |
223 | if (isFloat64Infinity(b) || isFloat64Zero(b)) { |
229 | if (isFloat64Infinity(b) || isFloat64Zero(b)) { |
224 | /*FIXME: inf / inf */ |
230 | /*FIXME: inf / inf */ |
225 | result.binary = FLOAT64_NAN; |
231 | result.binary = FLOAT64_NAN; |
226 | return result; |
232 | return result; |
227 | } |
233 | } |
228 | /* inf / num */ |
234 | /* inf / num */ |
229 | result.parts.exp = a.parts.exp; |
235 | result.parts.exp = a.parts.exp; |
230 | result.parts.fraction = a.parts.fraction; |
236 | result.parts.fraction = a.parts.fraction; |
231 | return result; |
237 | return result; |
232 | } |
238 | } |
233 | 239 | ||
234 | if (isFloat64Infinity(b)) { |
240 | if (isFloat64Infinity(b)) { |
235 | if (isFloat64Zero(a)) { |
241 | if (isFloat64Zero(a)) { |
236 | /* FIXME 0 / inf */ |
242 | /* FIXME 0 / inf */ |
237 | result.parts.exp = 0; |
243 | result.parts.exp = 0; |
238 | result.parts.fraction = 0; |
244 | result.parts.fraction = 0; |
239 | return result; |
245 | return result; |
240 | } |
246 | } |
241 | /* FIXME: num / inf*/ |
247 | /* FIXME: num / inf*/ |
242 | result.parts.exp = 0; |
248 | result.parts.exp = 0; |
243 | result.parts.fraction = 0; |
249 | result.parts.fraction = 0; |
244 | return result; |
250 | return result; |
245 | } |
251 | } |
246 | 252 | ||
247 | if (isFloat64Zero(b)) { |
253 | if (isFloat64Zero(b)) { |
248 | if (isFloat64Zero(a)) { |
254 | if (isFloat64Zero(a)) { |
249 | /*FIXME: 0 / 0*/ |
255 | /*FIXME: 0 / 0*/ |
250 | result.binary = FLOAT64_NAN; |
256 | result.binary = FLOAT64_NAN; |
251 | return result; |
257 | return result; |
252 | } |
258 | } |
253 | /* FIXME: division by zero */ |
259 | /* FIXME: division by zero */ |
254 | result.parts.exp = 0; |
260 | result.parts.exp = 0; |
255 | result.parts.fraction = 0; |
261 | result.parts.fraction = 0; |
256 | return result; |
262 | return result; |
257 | } |
263 | } |
258 | 264 | ||
259 | 265 | ||
260 | afrac = a.parts.fraction; |
266 | afrac = a.parts.fraction; |
261 | aexp = a.parts.exp; |
267 | aexp = a.parts.exp; |
262 | bfrac = b.parts.fraction; |
268 | bfrac = b.parts.fraction; |
263 | bexp = b.parts.exp; |
269 | bexp = b.parts.exp; |
264 | 270 | ||
265 | /* denormalized numbers */ |
271 | /* denormalized numbers */ |
266 | if (aexp == 0) { |
272 | if (aexp == 0) { |
267 | if (afrac == 0) { |
273 | if (afrac == 0) { |
268 | result.parts.exp = 0; |
274 | result.parts.exp = 0; |
269 | result.parts.fraction = 0; |
275 | result.parts.fraction = 0; |
270 | return result; |
276 | return result; |
271 | } |
277 | } |
272 | /* normalize it*/ |
278 | /* normalize it*/ |
273 | 279 | ||
274 | aexp++; |
280 | aexp++; |
275 | /* afrac is nonzero => it must stop */ |
281 | /* afrac is nonzero => it must stop */ |
276 | while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
282 | while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
277 | afrac <<= 1; |
283 | afrac <<= 1; |
278 | aexp--; |
284 | aexp--; |
279 | } |
285 | } |
280 | } |
286 | } |
281 | 287 | ||
282 | if (bexp == 0) { |
288 | if (bexp == 0) { |
283 | bexp++; |
289 | bexp++; |
284 | /* bfrac is nonzero => it must stop */ |
290 | /* bfrac is nonzero => it must stop */ |
285 | while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
291 | while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) { |
286 | bfrac <<= 1; |
292 | bfrac <<= 1; |
287 | bexp--; |
293 | bexp--; |
288 | } |
294 | } |
289 | } |
295 | } |
290 | 296 | ||
291 | afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 ); |
297 | afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 ); |
292 | bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1); |
298 | bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1); |
293 | 299 | ||
294 | if ( bfrac <= (afrac << 1) ) { |
300 | if ( bfrac <= (afrac << 1) ) { |
295 | afrac >>= 1; |
301 | afrac >>= 1; |
296 | aexp++; |
302 | aexp++; |
297 | } |
303 | } |
298 | 304 | ||
299 | cexp = aexp - bexp + FLOAT64_BIAS - 2; |
305 | cexp = aexp - bexp + FLOAT64_BIAS - 2; |
300 | 306 | ||
301 | cfrac = divFloat64estim(afrac, bfrac); |
307 | cfrac = divFloat64estim(afrac, bfrac); |
302 | 308 | ||
303 | if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */ |
309 | if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */ |
304 | mul64integers( bfrac, cfrac, &remlo, &remhi); |
310 | mul64integers( bfrac, cfrac, &remlo, &remhi); |
305 | /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ |
311 | /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ |
306 | remhi = afrac - remhi - ( remlo > 0); |
312 | remhi = afrac - remhi - ( remlo > 0); |
307 | remlo = - remlo; |
313 | remlo = - remlo; |
308 | 314 | ||
309 | while ((int64_t) remhi < 0) { |
315 | while ((int64_t) remhi < 0) { |
310 | cfrac--; |
316 | cfrac--; |
311 | remlo += bfrac; |
317 | remlo += bfrac; |
312 | remhi += ( remlo < bfrac ); |
318 | remhi += ( remlo < bfrac ); |
313 | } |
319 | } |
314 | cfrac |= ( remlo != 0 ); |
320 | cfrac |= ( remlo != 0 ); |
315 | } |
321 | } |
316 | 322 | ||
317 | /* round and shift */ |
323 | /* round and shift */ |
318 | result = finishFloat64(cexp, cfrac, result.parts.sign); |
324 | result = finishFloat64(cexp, cfrac, result.parts.sign); |
319 | return result; |
325 | return result; |
320 | 326 | ||
321 | } |
327 | } |
322 | 328 | ||
323 | uint64_t divFloat64estim(uint64_t a, uint64_t b) |
329 | uint64_t divFloat64estim(uint64_t a, uint64_t b) |
324 | { |
330 | { |
325 | uint64_t bhi; |
331 | uint64_t bhi; |
326 | uint64_t remhi, remlo; |
332 | uint64_t remhi, remlo; |
327 | uint64_t result; |
333 | uint64_t result; |
328 | 334 | ||
329 | if ( b <= a ) { |
335 | if ( b <= a ) { |
330 | return 0xFFFFFFFFFFFFFFFFull; |
336 | return 0xFFFFFFFFFFFFFFFFull; |
331 | } |
337 | } |
332 | 338 | ||
333 | bhi = b >> 32; |
339 | bhi = b >> 32; |
334 | result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32; |
340 | result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32; |
335 | mul64integers(b, result, &remlo, &remhi); |
341 | mul64integers(b, result, &remlo, &remhi); |
336 | 342 | ||
337 | remhi = a - remhi - (remlo > 0); |
343 | remhi = a - remhi - (remlo > 0); |
338 | remlo = - remlo; |
344 | remlo = - remlo; |
339 | 345 | ||
340 | b <<= 32; |
346 | b <<= 32; |
341 | while ( (int64_t) remhi < 0 ) { |
347 | while ( (int64_t) remhi < 0 ) { |
342 | result -= 0x1ll << 32; |
348 | result -= 0x1ll << 32; |
343 | remlo += b; |
349 | remlo += b; |
344 | remhi += bhi + ( remlo < b ); |
350 | remhi += bhi + ( remlo < b ); |
345 | } |
351 | } |
346 | remhi = (remhi << 32) | (remlo >> 32); |
352 | remhi = (remhi << 32) | (remlo >> 32); |
347 | if (( bhi << 32) <= remhi) { |
353 | if (( bhi << 32) <= remhi) { |
348 | result |= 0xFFFFFFFF; |
354 | result |= 0xFFFFFFFF; |
349 | } else { |
355 | } else { |
350 | result |= remhi / bhi; |
356 | result |= remhi / bhi; |
351 | } |
357 | } |
352 | 358 | ||
353 | 359 | ||
354 | return result; |
360 | return result; |
355 | } |
361 | } |
356 | 362 | ||
- | 363 | ||
- | 364 | /** @} |
|
- | 365 | */ |
|
- | 366 | ||
357 | 367 |