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

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

Fixed some problems with 64 bit arithmetic but others still persisting.
  • Property mode set to 100644
File size: 7.5 KB
RevLine 
[b5440cf]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
[12c6f2d]29#include<sftypes.h>
30#include<add.h>
[e6a40ac]31#include<div.h>
[12c6f2d]32#include<comparison.h>
[e6a40ac]33#include<mul.h>
[e979fea]34#include<common.h>
35
[b5440cf]36
[12c6f2d]37float32 divFloat32(float32 a, float32 b)
38{
[1266543]39 float32 result;
40 __s32 aexp, bexp, cexp;
41 __u64 afrac, bfrac, cfrac;
[12c6f2d]42
[1266543]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
[e6a40ac]147 /* find first nonzero digit and shift result and detect possibly underflow */
[1266543]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;
[12c6f2d]189}
[b5440cf]190
[e6a40ac]191float64 divFloat64(float64 a, float64 b)
192{
193 float64 result;
194 __s32 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 if (isFloat64SigNaN(a)) {
202 /*FIXME: SigNaN*/
203 }
204 /*NaN*/
205 return a;
206 }
207
208 if (isFloat64NaN(b)) {
209 if (isFloat64SigNaN(b)) {
210 /*FIXME: SigNaN*/
211 }
212 /*NaN*/
213 return b;
214 }
215
216 if (isFloat64Infinity(a)) {
217 if (isFloat64Infinity(b)) {
218 /*FIXME: inf / inf */
219 result.binary = FLOAT64_NAN;
220 return result;
221 }
222 /* inf / num */
223 result.parts.exp = a.parts.exp;
224 result.parts.fraction = a.parts.fraction;
225 return result;
226 }
227
228 if (isFloat64Infinity(b)) {
229 if (isFloat64Zero(a)) {
230 /* FIXME 0 / inf */
231 result.parts.exp = 0;
232 result.parts.fraction = 0;
233 return result;
234 }
235 /* FIXME: num / inf*/
236 result.parts.exp = 0;
237 result.parts.fraction = 0;
238 return result;
239 }
240
241 if (isFloat64Zero(b)) {
242 if (isFloat64Zero(a)) {
243 /*FIXME: 0 / 0*/
244 result.binary = FLOAT64_NAN;
245 return result;
246 }
247 /* FIXME: division by zero */
248 result.parts.exp = 0;
249 result.parts.fraction = 0;
250 return result;
251 }
252
253
254 afrac = a.parts.fraction;
255 aexp = a.parts.exp;
256 bfrac = b.parts.fraction;
257 bexp = b.parts.exp;
258
259 /* denormalized numbers */
260 if (aexp == 0) {
261 if (afrac == 0) {
262 result.parts.exp = 0;
263 result.parts.fraction = 0;
264 return result;
265 }
266 /* normalize it*/
267
268 afrac <<= 1;
269 /* afrac is nonzero => it must stop */
270 while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) {
271 afrac <<= 1;
272 aexp--;
273 }
274 }
275
276 if (bexp == 0) {
277 bfrac <<= 1;
278 /* bfrac is nonzero => it must stop */
279 while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) {
280 bfrac <<= 1;
281 bexp--;
282 }
283 }
284
285 afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 );
286 bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1);
287
288 if ( bfrac <= (afrac << 1) ) {
289 afrac >>= 1;
290 aexp++;
291 }
292
293 cexp = aexp - bexp + FLOAT64_BIAS - 2;
294
295 cfrac = divFloat64estim(afrac, bfrac);
296
297 if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */
298 mul64integers( bfrac, cfrac, &remlo, &remhi);
299 /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/
300 remhi = afrac - remhi - ( remlo > 0);
301 remlo = - remlo;
302
303 while ((__s64) remhi < 0) {
304 cfrac--;
305 remlo += bfrac;
306 remhi += ( remlo < bfrac );
307 }
308 cfrac |= ( remlo != 0 );
309 }
310
[e979fea]311 /* round and shift */
312 result = finishFloat64(cexp, cfrac, result.parts.sign);
313 return result;
[e6a40ac]314
315}
316
317__u64 divFloat64estim(__u64 a, __u64 b)
318{
319 __u64 bhi;
320 __u64 remhi, remlo;
321 __u64 result;
322
323 if ( b <= a ) {
324 return 0xFFFFFFFFFFFFFFFFull;
325 }
326
327 bhi = b >> 32;
328 result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32;
329 mul64integers(b, result, &remlo, &remhi);
330
331 remhi = a - remhi - (remlo > 0);
332 remlo = - remlo;
333
334 b <<= 32;
335 while ( (__s64) remhi < 0 ) {
336 result -= 0x1ll << 32;
337 remlo += b;
338 remhi += bhi + ( remlo < b );
339 }
340 remhi = (remhi << 32) | (remlo >> 32);
341 if (( bhi << 32) <= remhi) {
342 result |= 0xFFFFFFFF;
343 } else {
344 result |= remhi / bhi;
345 }
346
347
348 return result;
349}
350
Note: See TracBrowser for help on using the repository browser.