source: mainline/uspace/lib/softfloat/generic/div.c@ 6aae539d

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

softfloat redesign

  • avoid hardwired type sizes, actual sizes are determined at compile-time
  • add basic support for x87 extended-precision data types (stored as 96bit long double)
  • a lot of coding style changes (removal of CamelCase, etc.)
  • Property mode set to 100644
File size: 12.5 KB
RevLine 
[b5440cf]1/*
[df4ed85]2 * Copyright (c) 2005 Josef Cejka
[c67aff2]3 * Copyright (c) 2011 Petr Koupy
[b5440cf]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
[750636a]30/** @addtogroup softfloat
[846848a6]31 * @{
32 */
[c67aff2]33/** @file Division functions.
[846848a6]34 */
35
[750636a]36#include <sftypes.h>
37#include <add.h>
38#include <div.h>
39#include <comparison.h>
40#include <mul.h>
41#include <common.h>
[b5440cf]42
[88d5c1e]43/** Divide two single-precision floats.
44 *
[c67aff2]45 * @param a Nominator.
46 * @param b Denominator.
[88d5c1e]47 *
[c67aff2]48 * @return Result of division.
[88d5c1e]49 *
[c67aff2]50 */
[88d5c1e]51float32 div_float32(float32 a, float32 b)
[12c6f2d]52{
[1266543]53 float32 result;
[aa59fa0]54 int32_t aexp, bexp, cexp;
55 uint64_t afrac, bfrac, cfrac;
[12c6f2d]56
[1266543]57 result.parts.sign = a.parts.sign ^ b.parts.sign;
58
[88d5c1e]59 if (is_float32_nan(a)) {
60 if (is_float32_signan(a)) {
61 // FIXME: SigNaN
[1266543]62 }
[88d5c1e]63 /* NaN */
[1266543]64 return a;
65 }
66
[88d5c1e]67 if (is_float32_nan(b)) {
68 if (is_float32_signan(b)) {
69 // FIXME: SigNaN
[1266543]70 }
[88d5c1e]71 /* NaN */
[1266543]72 return b;
73 }
74
[88d5c1e]75 if (is_float32_infinity(a)) {
76 if (is_float32_infinity(b)) {
[1266543]77 /*FIXME: inf / inf */
[88d5c1e]78 result.bin = FLOAT32_NAN;
[1266543]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 }
[88d5c1e]86
87 if (is_float32_infinity(b)) {
88 if (is_float32_zero(a)) {
[1266543]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
[88d5c1e]100 if (is_float32_zero(b)) {
101 if (is_float32_zero(a)) {
[1266543]102 /*FIXME: 0 / 0*/
[88d5c1e]103 result.bin = FLOAT32_NAN;
[1266543]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) {
[c67aff2]120 result.parts.exp = 0;
121 result.parts.fraction = 0;
122 return result;
[1266543]123 }
[88d5c1e]124
[1266543]125 /* normalize it*/
126 afrac <<= 1;
[88d5c1e]127 /* afrac is nonzero => it must stop */
[c67aff2]128 while (!(afrac & FLOAT32_HIDDEN_BIT_MASK)) {
[1266543]129 afrac <<= 1;
130 aexp--;
131 }
132 }
[88d5c1e]133
[1266543]134 if (bexp == 0) {
135 bfrac <<= 1;
[88d5c1e]136 /* bfrac is nonzero => it must stop */
[c67aff2]137 while (!(bfrac & FLOAT32_HIDDEN_BIT_MASK)) {
[1266543]138 bfrac <<= 1;
139 bexp--;
140 }
141 }
[88d5c1e]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
[c67aff2]146 if (bfrac <= (afrac << 1)) {
[1266543]147 afrac >>= 1;
148 aexp++;
149 }
150
151 cexp = aexp - bexp + FLOAT32_BIAS - 2;
152
153 cfrac = (afrac << 32) / bfrac;
[c67aff2]154 if ((cfrac & 0x3F) == 0) {
155 cfrac |= (bfrac * cfrac != afrac << 32);
[1266543]156 }
157
158 /* pack and round */
159
[e6a40ac]160 /* find first nonzero digit and shift result and detect possibly underflow */
[c67aff2]161 while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)))) {
[1266543]162 cexp--;
163 cfrac <<= 1;
[c67aff2]164 /* TODO: fix underflow */
165 }
[1266543]166
167 cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/
168
169 if (cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)) {
170 ++cexp;
171 cfrac >>= 1;
[88d5c1e]172 }
173
[1266543]174 /* check overflow */
[c67aff2]175 if (cexp >= FLOAT32_MAX_EXPONENT) {
[1266543]176 /* FIXME: overflow, return infinity */
177 result.parts.exp = FLOAT32_MAX_EXPONENT;
178 result.parts.fraction = 0;
179 return result;
180 }
[88d5c1e]181
[1266543]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) {
[c67aff2]191 cexp++;
[1266543]192 cfrac >>= 1;
[88d5c1e]193 }
[1266543]194 } else {
[c67aff2]195 result.parts.exp = (uint32_t) cexp;
[1266543]196 }
197
198 result.parts.fraction = ((cfrac >> 6) & (~FLOAT32_HIDDEN_BIT_MASK));
199
[88d5c1e]200 return result;
[12c6f2d]201}
[b5440cf]202
[88d5c1e]203/** Divide two double-precision floats.
[c67aff2]204 *
205 * @param a Nominator.
206 * @param b Denominator.
[88d5c1e]207 *
[c67aff2]208 * @return Result of division.
[88d5c1e]209 *
[c67aff2]210 */
[88d5c1e]211float64 div_float64(float64 a, float64 b)
[e6a40ac]212{
213 float64 result;
[aa59fa0]214 int64_t aexp, bexp, cexp;
215 uint64_t afrac, bfrac, cfrac;
216 uint64_t remlo, remhi;
[c67aff2]217 uint64_t tmplo, tmphi;
[e6a40ac]218
219 result.parts.sign = a.parts.sign ^ b.parts.sign;
220
[88d5c1e]221 if (is_float64_nan(a)) {
222 if (is_float64_signan(b)) {
223 // FIXME: SigNaN
[f1f95f2]224 return b;
225 }
226
[88d5c1e]227 if (is_float64_signan(a)) {
228 // FIXME: SigNaN
[e6a40ac]229 }
[88d5c1e]230 /* NaN */
[e6a40ac]231 return a;
232 }
233
[88d5c1e]234 if (is_float64_nan(b)) {
235 if (is_float64_signan(b)) {
236 // FIXME: SigNaN
[e6a40ac]237 }
[88d5c1e]238 /* NaN */
[e6a40ac]239 return b;
240 }
241
[88d5c1e]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;
[e6a40ac]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 }
[88d5c1e]253
254 if (is_float64_infinity(b)) {
255 if (is_float64_zero(a)) {
[e6a40ac]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
[88d5c1e]267 if (is_float64_zero(b)) {
268 if (is_float64_zero(a)) {
[e6a40ac]269 /*FIXME: 0 / 0*/
[88d5c1e]270 result.bin = FLOAT64_NAN;
[e6a40ac]271 return result;
272 }
273 /* FIXME: division by zero */
274 result.parts.exp = 0;
275 result.parts.fraction = 0;
276 return result;
277 }
[88d5c1e]278
[e6a40ac]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) {
[f1f95f2]287 result.parts.exp = 0;
288 result.parts.fraction = 0;
289 return result;
[e6a40ac]290 }
[88d5c1e]291
[e6a40ac]292 /* normalize it*/
[f1f95f2]293 aexp++;
[c67aff2]294 /* afrac is nonzero => it must stop */
295 while (!(afrac & FLOAT64_HIDDEN_BIT_MASK)) {
[e6a40ac]296 afrac <<= 1;
297 aexp--;
298 }
299 }
[88d5c1e]300
[e6a40ac]301 if (bexp == 0) {
[f1f95f2]302 bexp++;
[c67aff2]303 /* bfrac is nonzero => it must stop */
304 while (!(bfrac & FLOAT64_HIDDEN_BIT_MASK)) {
[e6a40ac]305 bfrac <<= 1;
306 bexp--;
307 }
308 }
[88d5c1e]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
[c67aff2]313 if (bfrac <= (afrac << 1)) {
[e6a40ac]314 afrac >>= 1;
315 aexp++;
316 }
317
318 cexp = aexp - bexp + FLOAT64_BIAS - 2;
319
[c67aff2]320 cfrac = div128est(afrac, 0x0ll, bfrac);
[e6a40ac]321
[c67aff2]322 if ((cfrac & 0x1FF) <= 2) {
323 mul64(bfrac, cfrac, &tmphi, &tmplo);
324 sub128(afrac, 0x0ll, tmphi, tmplo, &remhi, &remlo);
[e6a40ac]325
[aa59fa0]326 while ((int64_t) remhi < 0) {
[e6a40ac]327 cfrac--;
[c67aff2]328 add128(remhi, remlo, 0x0ll, bfrac, &remhi, &remlo);
[e6a40ac]329 }
[c67aff2]330 cfrac |= (remlo != 0);
[e6a40ac]331 }
332
[e979fea]333 /* round and shift */
[88d5c1e]334 result = finish_float64(cexp, cfrac, result.parts.sign);
[e979fea]335 return result;
[e6a40ac]336}
337
[88d5c1e]338/** Divide two quadruple-precision floats.
[c67aff2]339 *
340 * @param a Nominator.
341 * @param b Denominator.
[88d5c1e]342 *
[c67aff2]343 * @return Result of division.
[88d5c1e]344 *
[c67aff2]345 */
[88d5c1e]346float128 div_float128(float128 a, float128 b)
[e6a40ac]347{
[c67aff2]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;
[88d5c1e]354
[c67aff2]355 result.parts.sign = a.parts.sign ^ b.parts.sign;
[88d5c1e]356
357 if (is_float128_nan(a)) {
358 if (is_float128_signan(b)) {
359 // FIXME: SigNaN
[c67aff2]360 return b;
361 }
[88d5c1e]362
363 if (is_float128_signan(a)) {
364 // FIXME: SigNaN
[c67aff2]365 }
[88d5c1e]366 /* NaN */
[c67aff2]367 return a;
[e6a40ac]368 }
[88d5c1e]369
370 if (is_float128_nan(b)) {
371 if (is_float128_signan(b)) {
372 // FIXME: SigNaN
[e6a40ac]373 }
[88d5c1e]374 /* NaN */
[c67aff2]375 return b;
[e6a40ac]376 }
[88d5c1e]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;
[c67aff2]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 }
[88d5c1e]391
392 if (is_float128_infinity(b)) {
393 if (is_float128_zero(a)) {
394 // FIXME 0 / inf
[c67aff2]395 result.parts.exp = 0;
396 result.parts.frac_hi = 0;
397 result.parts.frac_lo = 0;
398 return result;
399 }
[88d5c1e]400 // FIXME: num / inf
[c67aff2]401 result.parts.exp = 0;
402 result.parts.frac_hi = 0;
403 result.parts.frac_lo = 0;
404 return result;
405 }
[88d5c1e]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;
[c67aff2]412 return result;
413 }
[88d5c1e]414 // FIXME: division by zero
[c67aff2]415 result.parts.exp = 0;
416 result.parts.frac_hi = 0;
417 result.parts.frac_lo = 0;
418 return result;
419 }
[88d5c1e]420
[c67aff2]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;
[88d5c1e]427
[c67aff2]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 }
[88d5c1e]436
[c67aff2]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 }
[88d5c1e]448
[c67aff2]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 }
[88d5c1e]460
[c67aff2]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);
[88d5c1e]471
[c67aff2]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 }
[88d5c1e]476
[c67aff2]477 cexp = aexp - bexp + FLOAT128_BIAS - 2;
[88d5c1e]478
[c67aff2]479 cfrac_hi = div128est(afrac_hi, afrac_lo, bfrac_hi);
[88d5c1e]480
[c67aff2]481 mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_hi,
482 &tmp_lolo /* dummy */, &tmp_hihi, &tmp_hilo, &tmp_lohi);
[88d5c1e]483
484 /* sub192(afrac_hi, afrac_lo, 0,
[c67aff2]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;
[88d5c1e]492
[c67aff2]493 while ((int64_t) rem_hihi < 0) {
494 --cfrac_hi;
[88d5c1e]495 /* add192(rem_hihi, rem_hilo, rem_lohi,
[c67aff2]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 }
[88d5c1e]503
[c67aff2]504 cfrac_lo = div128est(rem_hilo, rem_lohi, bfrac_lo);
[88d5c1e]505
[c67aff2]506 if ((cfrac_lo & 0x3FFF) <= 4) {
507 mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_lo,
[88d5c1e]508 &tmp_hihi /* dummy */, &tmp_hilo, &tmp_lohi, &tmp_lolo);
509
[c67aff2]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;
[88d5c1e]518
[c67aff2]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 }
[88d5c1e]529
[c67aff2]530 cfrac_lo |= ((rem_hilo | rem_lohi | rem_lolo) != 0 );
531 }
[88d5c1e]532
[c67aff2]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);
[88d5c1e]536
537 result = finish_float128(cexp, cfrac_hi, cfrac_lo, result.parts.sign, shift_out);
[e6a40ac]538 return result;
539}
540
[231a60a]541/** @}
[846848a6]542 */
Note: See TracBrowser for help on using the repository browser.