source: mainline/softfloat/generic/div.c@ d9f51ccc

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since d9f51ccc was f1f95f2, checked in by Josef Cejka <malyzelenyhnus@…>, 19 years ago

Fixed bugs in 64bit float division.
  • Property mode set to 100644
File size: 7.6 KB
Line 
1/*
2 * Copyright (C) 2005 Josef Cejka
3 * All rights reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
7 * are met:
8 *
9 * - Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 * - Redistributions in binary form must reproduce the above copyright
12 * notice, this list of conditions and the following disclaimer in the
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
15 * derived from this software without specific prior written permission.
16 *
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
19 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
20 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
21 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
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
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
26 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
27 */
28
29#include<sftypes.h>
30#include<add.h>
31#include<div.h>
32#include<comparison.h>
33#include<mul.h>
34#include<common.h>
35
36
37float32 divFloat32(float32 a, float32 b)
38{
39 float32 result;
40 __s32 aexp, bexp, cexp;
41 __u64 afrac, bfrac, cfrac;
42
43 result.parts.sign = a.parts.sign ^ b.parts.sign;
44
45 if (isFloat32NaN(a)) {
46 if (isFloat32SigNaN(a)) {
47 /*FIXME: SigNaN*/
48 }
49 /*NaN*/
50 return a;
51 }
52
53 if (isFloat32NaN(b)) {
54 if (isFloat32SigNaN(b)) {
55 /*FIXME: SigNaN*/
56 }
57 /*NaN*/
58 return b;
59 }
60
61 if (isFloat32Infinity(a)) {
62 if (isFloat32Infinity(b)) {
63 /*FIXME: inf / inf */
64 result.binary = FLOAT32_NAN;
65 return result;
66 }
67 /* inf / num */
68 result.parts.exp = a.parts.exp;
69 result.parts.fraction = a.parts.fraction;
70 return result;
71 }
72
73 if (isFloat32Infinity(b)) {
74 if (isFloat32Zero(a)) {
75 /* FIXME 0 / inf */
76 result.parts.exp = 0;
77 result.parts.fraction = 0;
78 return result;
79 }
80 /* FIXME: num / inf*/
81 result.parts.exp = 0;
82 result.parts.fraction = 0;
83 return result;
84 }
85
86 if (isFloat32Zero(b)) {
87 if (isFloat32Zero(a)) {
88 /*FIXME: 0 / 0*/
89 result.binary = FLOAT32_NAN;
90 return result;
91 }
92 /* FIXME: division by zero */
93 result.parts.exp = 0;
94 result.parts.fraction = 0;
95 return result;
96 }
97
98
99 afrac = a.parts.fraction;
100 aexp = a.parts.exp;
101 bfrac = b.parts.fraction;
102 bexp = b.parts.exp;
103
104 /* denormalized numbers */
105 if (aexp == 0) {
106 if (afrac == 0) {
107 result.parts.exp = 0;
108 result.parts.fraction = 0;
109 return result;
110 }
111 /* normalize it*/
112
113 afrac <<= 1;
114 /* afrac is nonzero => it must stop */
115 while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) {
116 afrac <<= 1;
117 aexp--;
118 }
119 }
120
121 if (bexp == 0) {
122 bfrac <<= 1;
123 /* bfrac is nonzero => it must stop */
124 while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) {
125 bfrac <<= 1;
126 bexp--;
127 }
128 }
129
130 afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 );
131 bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE );
132
133 if ( bfrac <= (afrac << 1) ) {
134 afrac >>= 1;
135 aexp++;
136 }
137
138 cexp = aexp - bexp + FLOAT32_BIAS - 2;
139
140 cfrac = (afrac << 32) / bfrac;
141 if (( cfrac & 0x3F ) == 0) {
142 cfrac |= ( bfrac * cfrac != afrac << 32 );
143 }
144
145 /* pack and round */
146
147 /* find first nonzero digit and shift result and detect possibly underflow */
148 while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) {
149 cexp--;
150 cfrac <<= 1;
151 /* TODO: fix underflow */
152 };
153
154 cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
155
156 if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
157 ++cexp;
158 cfrac >>= 1;
159 }
160
161 /* check overflow */
162 if (cexp >= FLOAT32_MAX_EXPONENT ) {
163 /* FIXME: overflow, return infinity */
164 result.parts.exp = FLOAT32_MAX_EXPONENT;
165 result.parts.fraction = 0;
166 return result;
167 }
168
169 if (cexp < 0) {
170 /* FIXME: underflow */
171 result.parts.exp = 0;
172 if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
173 result.parts.fraction = 0;
174 return result;
175 }
176 cfrac >>= 1;
177 while (cexp < 0) {
178 cexp ++;
179 cfrac >>= 1;
180 }
181
182 } else {
183 result.parts.exp = (__u32)cexp;
184 }
185
186 result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
187
188 return result;
189}
190
191float64 divFloat64(float64 a, float64 b)
192{
193 float64 result;
194 __s64 aexp, bexp, cexp;
195 __u64 afrac, bfrac, cfrac;
196 __u64 remlo, remhi;
197
198 result.parts.sign = a.parts.sign ^ b.parts.sign;
199
200 if (isFloat64NaN(a)) {
201
202 if (isFloat64SigNaN(b)) {
203 /*FIXME: SigNaN*/
204 return b;
205 }
206
207 if (isFloat64SigNaN(a)) {
208 /*FIXME: SigNaN*/
209 }
210 /*NaN*/
211 return a;
212 }
213
214 if (isFloat64NaN(b)) {
215 if (isFloat64SigNaN(b)) {
216 /*FIXME: SigNaN*/
217 }
218 /*NaN*/
219 return b;
220 }
221
222 if (isFloat64Infinity(a)) {
223 if (isFloat64Infinity(b) || isFloat64Zero(b)) {
224 /*FIXME: inf / inf */
225 result.binary = FLOAT64_NAN;
226 return result;
227 }
228 /* inf / num */
229 result.parts.exp = a.parts.exp;
230 result.parts.fraction = a.parts.fraction;
231 return result;
232 }
233
234 if (isFloat64Infinity(b)) {
235 if (isFloat64Zero(a)) {
236 /* FIXME 0 / inf */
237 result.parts.exp = 0;
238 result.parts.fraction = 0;
239 return result;
240 }
241 /* FIXME: num / inf*/
242 result.parts.exp = 0;
243 result.parts.fraction = 0;
244 return result;
245 }
246
247 if (isFloat64Zero(b)) {
248 if (isFloat64Zero(a)) {
249 /*FIXME: 0 / 0*/
250 result.binary = FLOAT64_NAN;
251 return result;
252 }
253 /* FIXME: division by zero */
254 result.parts.exp = 0;
255 result.parts.fraction = 0;
256 return result;
257 }
258
259
260 afrac = a.parts.fraction;
261 aexp = a.parts.exp;
262 bfrac = b.parts.fraction;
263 bexp = b.parts.exp;
264
265 /* denormalized numbers */
266 if (aexp == 0) {
267 if (afrac == 0) {
268 result.parts.exp = 0;
269 result.parts.fraction = 0;
270 return result;
271 }
272 /* normalize it*/
273
274 aexp++;
275 /* afrac is nonzero => it must stop */
276 while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
277 afrac <<= 1;
278 aexp--;
279 }
280 }
281
282 if (bexp == 0) {
283 bexp++;
284 /* bfrac is nonzero => it must stop */
285 while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
286 bfrac <<= 1;
287 bexp--;
288 }
289 }
290
291 afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
292 bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
293
294 if ( bfrac <= (afrac << 1) ) {
295 afrac >>= 1;
296 aexp++;
297 }
298
299 cexp = aexp - bexp + FLOAT64_BIAS - 2;
300
301 cfrac = divFloat64estim(afrac, bfrac);
302
303 if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
304 mul64integers( bfrac, cfrac, &remlo, &remhi);
305 /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/
306 remhi = afrac - remhi - ( remlo > 0);
307 remlo = - remlo;
308
309 while ((__s64) remhi < 0) {
310 cfrac--;
311 remlo += bfrac;
312 remhi += ( remlo < bfrac );
313 }
314 cfrac |= ( remlo != 0 );
315 }
316
317 /* round and shift */
318 result = finishFloat64(cexp, cfrac, result.parts.sign);
319 return result;
320
321}
322
323__u64 divFloat64estim(__u64 a, __u64 b)
324{
325 __u64 bhi;
326 __u64 remhi, remlo;
327 __u64 result;
328
329 if ( b <= a ) {
330 return 0xFFFFFFFFFFFFFFFFull;
331 }
332
333 bhi = b >> 32;
334 result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
335 mul64integers(b, result, &remlo, &remhi);
336
337 remhi = a - remhi - (remlo > 0);
338 remlo = - remlo;
339
340 b <<= 32;
341 while ( (__s64) remhi < 0 ) {
342 result -= 0x1ll << 32;
343 remlo += b;
344 remhi += bhi + ( remlo < b );
345 }
346 remhi = (remhi << 32) | (remlo >> 32);
347 if (( bhi << 32) <= remhi) {
348 result |= 0xFFFFFFFF;
349 } else {
350 result |= remhi / bhi;
351 }
352
353
354 return result;
355}
356
Note: See TracBrowser for help on using the repository browser.