source: mainline/uspace/lib/softfloat/div.c@ b3cf946

lfn serial ticket/834-toolchain-update topic/msim-upgrade topic/simplify-dev-export
Last change on this file since b3cf946 was 2416085, checked in by Martin Decky <martin@…>, 13 years ago

softfloat: move sources to a single directory
implement more ARM EABI bindings
  • Property mode set to 100644
File size: 12.5 KB
Line 
1/*
2 * Copyright (c) 2005 Josef Cejka
3 * Copyright (c) 2011 Petr Koupy
4 * All rights reserved.
5 *
6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions
8 * are met:
9 *
10 * - Redistributions of source code must retain the above copyright
11 * notice, this list of conditions and the following disclaimer.
12 * - Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the distribution.
15 * - The name of the author may not be used to endorse or promote products
16 * derived from this software without specific prior written permission.
17 *
18 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
19 * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
20 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
21 * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
22 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
23 * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
24 * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
25 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
26 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
27 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28 */
29
30/** @addtogroup softfloat
31 * @{
32 */
33/** @file Division functions.
34 */
35
36#include "sftypes.h"
37#include "add.h"
38#include "div.h"
39#include "comparison.h"
40#include "mul.h"
41#include "common.h"
42
43/** Divide two single-precision floats.
44 *
45 * @param a Nominator.
46 * @param b Denominator.
47 *
48 * @return Result of division.
49 *
50 */
51float32 div_float32(float32 a, float32 b)
52{
53 float32 result;
54 int32_t aexp, bexp, cexp;
55 uint64_t afrac, bfrac, cfrac;
56
57 result.parts.sign = a.parts.sign ^ b.parts.sign;
58
59 if (is_float32_nan(a)) {
60 if (is_float32_signan(a)) {
61 // FIXME: SigNaN
62 }
63 /* NaN */
64 return a;
65 }
66
67 if (is_float32_nan(b)) {
68 if (is_float32_signan(b)) {
69 // FIXME: SigNaN
70 }
71 /* NaN */
72 return b;
73 }
74
75 if (is_float32_infinity(a)) {
76 if (is_float32_infinity(b)) {
77 /*FIXME: inf / inf */
78 result.bin = FLOAT32_NAN;
79 return result;
80 }
81 /* inf / num */
82 result.parts.exp = a.parts.exp;
83 result.parts.fraction = a.parts.fraction;
84 return result;
85 }
86
87 if (is_float32_infinity(b)) {
88 if (is_float32_zero(a)) {
89 /* FIXME 0 / inf */
90 result.parts.exp = 0;
91 result.parts.fraction = 0;
92 return result;
93 }
94 /* FIXME: num / inf*/
95 result.parts.exp = 0;
96 result.parts.fraction = 0;
97 return result;
98 }
99
100 if (is_float32_zero(b)) {
101 if (is_float32_zero(a)) {
102 /*FIXME: 0 / 0*/
103 result.bin = FLOAT32_NAN;
104 return result;
105 }
106 /* FIXME: division by zero */
107 result.parts.exp = 0;
108 result.parts.fraction = 0;
109 return result;
110 }
111
112 afrac = a.parts.fraction;
113 aexp = a.parts.exp;
114 bfrac = b.parts.fraction;
115 bexp = b.parts.exp;
116
117 /* denormalized numbers */
118 if (aexp == 0) {
119 if (afrac == 0) {
120 result.parts.exp = 0;
121 result.parts.fraction = 0;
122 return result;
123 }
124
125 /* normalize it*/
126 afrac <<= 1;
127 /* afrac is nonzero => it must stop */
128 while (!(afrac & FLOAT32_HIDDEN_BIT_MASK)) {
129 afrac <<= 1;
130 aexp--;
131 }
132 }
133
134 if (bexp == 0) {
135 bfrac <<= 1;
136 /* bfrac is nonzero => it must stop */
137 while (!(bfrac & FLOAT32_HIDDEN_BIT_MASK)) {
138 bfrac <<= 1;
139 bexp--;
140 }
141 }
142
143 afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE - 1);
144 bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE);
145
146 if (bfrac <= (afrac << 1)) {
147 afrac >>= 1;
148 aexp++;
149 }
150
151 cexp = aexp - bexp + FLOAT32_BIAS - 2;
152
153 cfrac = (afrac << 32) / bfrac;
154 if ((cfrac & 0x3F) == 0) {
155 cfrac |= (bfrac * cfrac != afrac << 32);
156 }
157
158 /* pack and round */
159
160 /* find first nonzero digit and shift result and detect possibly underflow */
161 while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)))) {
162 cexp--;
163 cfrac <<= 1;
164 /* TODO: fix underflow */
165 }
166
167 cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
168
169 if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
170 ++cexp;
171 cfrac >>= 1;
172 }
173
174 /* check overflow */
175 if (cexp >= FLOAT32_MAX_EXPONENT) {
176 /* FIXME: overflow, return infinity */
177 result.parts.exp = FLOAT32_MAX_EXPONENT;
178 result.parts.fraction = 0;
179 return result;
180 }
181
182 if (cexp < 0) {
183 /* FIXME: underflow */
184 result.parts.exp = 0;
185 if ((cexp + FLOAT32_FRACTION_SIZE) < 0) {
186 result.parts.fraction = 0;
187 return result;
188 }
189 cfrac >>= 1;
190 while (cexp < 0) {
191 cexp++;
192 cfrac >>= 1;
193 }
194 } else {
195 result.parts.exp = (uint32_t) cexp;
196 }
197
198 result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
199
200 return result;
201}
202
203/** Divide two double-precision floats.
204 *
205 * @param a Nominator.
206 * @param b Denominator.
207 *
208 * @return Result of division.
209 *
210 */
211float64 div_float64(float64 a, float64 b)
212{
213 float64 result;
214 int64_t aexp, bexp, cexp;
215 uint64_t afrac, bfrac, cfrac;
216 uint64_t remlo, remhi;
217 uint64_t tmplo, tmphi;
218
219 result.parts.sign = a.parts.sign ^ b.parts.sign;
220
221 if (is_float64_nan(a)) {
222 if (is_float64_signan(b)) {
223 // FIXME: SigNaN
224 return b;
225 }
226
227 if (is_float64_signan(a)) {
228 // FIXME: SigNaN
229 }
230 /* NaN */
231 return a;
232 }
233
234 if (is_float64_nan(b)) {
235 if (is_float64_signan(b)) {
236 // FIXME: SigNaN
237 }
238 /* NaN */
239 return b;
240 }
241
242 if (is_float64_infinity(a)) {
243 if (is_float64_infinity(b) || is_float64_zero(b)) {
244 // FIXME: inf / inf
245 result.bin = FLOAT64_NAN;
246 return result;
247 }
248 /* inf / num */
249 result.parts.exp = a.parts.exp;
250 result.parts.fraction = a.parts.fraction;
251 return result;
252 }
253
254 if (is_float64_infinity(b)) {
255 if (is_float64_zero(a)) {
256 /* FIXME 0 / inf */
257 result.parts.exp = 0;
258 result.parts.fraction = 0;
259 return result;
260 }
261 /* FIXME: num / inf*/
262 result.parts.exp = 0;
263 result.parts.fraction = 0;
264 return result;
265 }
266
267 if (is_float64_zero(b)) {
268 if (is_float64_zero(a)) {
269 /*FIXME: 0 / 0*/
270 result.bin = FLOAT64_NAN;
271 return result;
272 }
273 /* FIXME: division by zero */
274 result.parts.exp = 0;
275 result.parts.fraction = 0;
276 return result;
277 }
278
279 afrac = a.parts.fraction;
280 aexp = a.parts.exp;
281 bfrac = b.parts.fraction;
282 bexp = b.parts.exp;
283
284 /* denormalized numbers */
285 if (aexp == 0) {
286 if (afrac == 0) {
287 result.parts.exp = 0;
288 result.parts.fraction = 0;
289 return result;
290 }
291
292 /* normalize it*/
293 aexp++;
294 /* afrac is nonzero => it must stop */
295 while (!(afrac & FLOAT64_HIDDEN_BIT_MASK)) {
296 afrac <<= 1;
297 aexp--;
298 }
299 }
300
301 if (bexp == 0) {
302 bexp++;
303 /* bfrac is nonzero => it must stop */
304 while (!(bfrac & FLOAT64_HIDDEN_BIT_MASK)) {
305 bfrac <<= 1;
306 bexp--;
307 }
308 }
309
310 afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 2);
311 bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 1);
312
313 if (bfrac <= (afrac << 1)) {
314 afrac >>= 1;
315 aexp++;
316 }
317
318 cexp = aexp - bexp + FLOAT64_BIAS - 2;
319
320 cfrac = div128est(afrac, 0x0ll, bfrac);
321
322 if ((cfrac & 0x1FF) <= 2) {
323 mul64(bfrac, cfrac, &tmphi, &tmplo);
324 sub128(afrac, 0x0ll, tmphi, tmplo, &remhi, &remlo);
325
326 while ((int64_t) remhi < 0) {
327 cfrac--;
328 add128(remhi, remlo, 0x0ll, bfrac, &remhi, &remlo);
329 }
330 cfrac |= (remlo != 0);
331 }
332
333 /* round and shift */
334 result = finish_float64(cexp, cfrac, result.parts.sign);
335 return result;
336}
337
338/** Divide two quadruple-precision floats.
339 *
340 * @param a Nominator.
341 * @param b Denominator.
342 *
343 * @return Result of division.
344 *
345 */
346float128 div_float128(float128 a, float128 b)
347{
348 float128 result;
349 int64_t aexp, bexp, cexp;
350 uint64_t afrac_hi, afrac_lo, bfrac_hi, bfrac_lo, cfrac_hi, cfrac_lo;
351 uint64_t shift_out;
352 uint64_t rem_hihi, rem_hilo, rem_lohi, rem_lolo;
353 uint64_t tmp_hihi, tmp_hilo, tmp_lohi, tmp_lolo;
354
355 result.parts.sign = a.parts.sign ^ b.parts.sign;
356
357 if (is_float128_nan(a)) {
358 if (is_float128_signan(b)) {
359 // FIXME: SigNaN
360 return b;
361 }
362
363 if (is_float128_signan(a)) {
364 // FIXME: SigNaN
365 }
366 /* NaN */
367 return a;
368 }
369
370 if (is_float128_nan(b)) {
371 if (is_float128_signan(b)) {
372 // FIXME: SigNaN
373 }
374 /* NaN */
375 return b;
376 }
377
378 if (is_float128_infinity(a)) {
379 if (is_float128_infinity(b) || is_float128_zero(b)) {
380 // FIXME: inf / inf
381 result.bin.hi = FLOAT128_NAN_HI;
382 result.bin.lo = FLOAT128_NAN_LO;
383 return result;
384 }
385 /* inf / num */
386 result.parts.exp = a.parts.exp;
387 result.parts.frac_hi = a.parts.frac_hi;
388 result.parts.frac_lo = a.parts.frac_lo;
389 return result;
390 }
391
392 if (is_float128_infinity(b)) {
393 if (is_float128_zero(a)) {
394 // FIXME 0 / inf
395 result.parts.exp = 0;
396 result.parts.frac_hi = 0;
397 result.parts.frac_lo = 0;
398 return result;
399 }
400 // FIXME: num / inf
401 result.parts.exp = 0;
402 result.parts.frac_hi = 0;
403 result.parts.frac_lo = 0;
404 return result;
405 }
406
407 if (is_float128_zero(b)) {
408 if (is_float128_zero(a)) {
409 // FIXME: 0 / 0
410 result.bin.hi = FLOAT128_NAN_HI;
411 result.bin.lo = FLOAT128_NAN_LO;
412 return result;
413 }
414 // FIXME: division by zero
415 result.parts.exp = 0;
416 result.parts.frac_hi = 0;
417 result.parts.frac_lo = 0;
418 return result;
419 }
420
421 afrac_hi = a.parts.frac_hi;
422 afrac_lo = a.parts.frac_lo;
423 aexp = a.parts.exp;
424 bfrac_hi = b.parts.frac_hi;
425 bfrac_lo = b.parts.frac_lo;
426 bexp = b.parts.exp;
427
428 /* denormalized numbers */
429 if (aexp == 0) {
430 if (eq128(afrac_hi, afrac_lo, 0x0ll, 0x0ll)) {
431 result.parts.exp = 0;
432 result.parts.frac_hi = 0;
433 result.parts.frac_lo = 0;
434 return result;
435 }
436
437 /* normalize it*/
438 aexp++;
439 /* afrac is nonzero => it must stop */
440 and128(afrac_hi, afrac_lo,
441 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
442 &tmp_hihi, &tmp_lolo);
443 while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) {
444 lshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo);
445 aexp--;
446 }
447 }
448
449 if (bexp == 0) {
450 bexp++;
451 /* bfrac is nonzero => it must stop */
452 and128(bfrac_hi, bfrac_lo,
453 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
454 &tmp_hihi, &tmp_lolo);
455 while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) {
456 lshift128(bfrac_hi, bfrac_lo, 1, &bfrac_hi, &bfrac_lo);
457 bexp--;
458 }
459 }
460
461 or128(afrac_hi, afrac_lo,
462 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
463 &afrac_hi, &afrac_lo);
464 lshift128(afrac_hi, afrac_lo,
465 (128 - FLOAT128_FRACTION_SIZE - 1), &afrac_hi, &afrac_lo);
466 or128(bfrac_hi, bfrac_lo,
467 FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO,
468 &bfrac_hi, &bfrac_lo);
469 lshift128(bfrac_hi, bfrac_lo,
470 (128 - FLOAT128_FRACTION_SIZE - 1), &bfrac_hi, &bfrac_lo);
471
472 if (le128(bfrac_hi, bfrac_lo, afrac_hi, afrac_lo)) {
473 rshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo);
474 aexp++;
475 }
476
477 cexp = aexp - bexp + FLOAT128_BIAS - 2;
478
479 cfrac_hi = div128est(afrac_hi, afrac_lo, bfrac_hi);
480
481 mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_hi,
482 &tmp_lolo /* dummy */, &tmp_hihi, &tmp_hilo, &tmp_lohi);
483
484 /* sub192(afrac_hi, afrac_lo, 0,
485 * tmp_hihi, tmp_hilo, tmp_lohi
486 * &rem_hihi, &rem_hilo, &rem_lohi); */
487 sub128(afrac_hi, afrac_lo, tmp_hihi, tmp_hilo, &rem_hihi, &rem_hilo);
488 if (tmp_lohi > 0) {
489 sub128(rem_hihi, rem_hilo, 0x0ll, 0x1ll, &rem_hihi, &rem_hilo);
490 }
491 rem_lohi = -tmp_lohi;
492
493 while ((int64_t) rem_hihi < 0) {
494 --cfrac_hi;
495 /* add192(rem_hihi, rem_hilo, rem_lohi,
496 * 0, bfrac_hi, bfrac_lo,
497 * &rem_hihi, &rem_hilo, &rem_lohi); */
498 add128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo, &rem_hilo, &rem_lohi);
499 if (lt128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo)) {
500 ++rem_hihi;
501 }
502 }
503
504 cfrac_lo = div128est(rem_hilo, rem_lohi, bfrac_lo);
505
506 if ((cfrac_lo & 0x3FFF) <= 4) {
507 mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_lo,
508 &tmp_hihi /* dummy */, &tmp_hilo, &tmp_lohi, &tmp_lolo);
509
510 /* sub192(rem_hilo, rem_lohi, 0,
511 * tmp_hilo, tmp_lohi, tmp_lolo,
512 * &rem_hilo, &rem_lohi, &rem_lolo); */
513 sub128(rem_hilo, rem_lohi, tmp_hilo, tmp_lohi, &rem_hilo, &rem_lohi);
514 if (tmp_lolo > 0) {
515 sub128(rem_hilo, rem_lohi, 0x0ll, 0x1ll, &rem_hilo, &rem_lohi);
516 }
517 rem_lolo = -tmp_lolo;
518
519 while ((int64_t) rem_hilo < 0) {
520 --cfrac_lo;
521 /* add192(rem_hilo, rem_lohi, rem_lolo,
522 * 0, bfrac_hi, bfrac_lo,
523 * &rem_hilo, &rem_lohi, &rem_lolo); */
524 add128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo, &rem_lohi, &rem_lolo);
525 if (lt128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo)) {
526 ++rem_hilo;
527 }
528 }
529
530 cfrac_lo |= ((rem_hilo | rem_lohi | rem_lolo) != 0 );
531 }
532
533 shift_out = cfrac_lo << (64 - (128 - FLOAT128_FRACTION_SIZE - 1));
534 rshift128(cfrac_hi, cfrac_lo, (128 - FLOAT128_FRACTION_SIZE - 1),
535 &cfrac_hi, &cfrac_lo);
536
537 result = finish_float128(cexp, cfrac_hi, cfrac_lo, result.parts.sign, shift_out);
538 return result;
539}
540
541/** @}
542 */
Note: See TracBrowser for help on using the repository browser.