Index: uspace/lib/softfloat/generic/div.c =================================================================== --- uspace/lib/softfloat/generic/div.c (revision 750636ac4a65360c44846b7681af6ae46854dd3a) +++ uspace/lib/softfloat/generic/div.c (revision 4863bb52d2de7a13a21dab7be2475975a848d49f) @@ -1,4 +1,5 @@ /* * Copyright (c) 2005 Josef Cejka + * Copyright (c) 2011 Petr Koupy * All rights reserved. * @@ -30,5 +31,5 @@ * @{ */ -/** @file +/** @file Division functions. */ @@ -40,4 +41,11 @@ #include +/** + * Divide two single-precision floats. + * + * @param a Nominator. + * @param b Denominator. + * @return Result of division. + */ float32 divFloat32(float32 a, float32 b) { @@ -100,5 +108,4 @@ return result; } - afrac = a.parts.fraction; @@ -110,13 +117,13 @@ if (aexp == 0) { if (afrac == 0) { - result.parts.exp = 0; - result.parts.fraction = 0; - return result; - } + result.parts.exp = 0; + result.parts.fraction = 0; + return result; + } + /* normalize it*/ - afrac <<= 1; - /* afrac is nonzero => it must stop */ - while (! (afrac & FLOAT32_HIDDEN_BIT_MASK) ) { + /* afrac is nonzero => it must stop */ + while (!(afrac & FLOAT32_HIDDEN_BIT_MASK)) { afrac <<= 1; aexp--; @@ -126,6 +133,6 @@ if (bexp == 0) { bfrac <<= 1; - /* bfrac is nonzero => it must stop */ - while (! (bfrac & FLOAT32_HIDDEN_BIT_MASK) ) { + /* bfrac is nonzero => it must stop */ + while (!(bfrac & FLOAT32_HIDDEN_BIT_MASK)) { bfrac <<= 1; bexp--; @@ -133,8 +140,8 @@ } - afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE - 1 ); - bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK ) << (32 - FLOAT32_FRACTION_SIZE ); - - if ( bfrac <= (afrac << 1) ) { + afrac = (afrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE - 1); + bfrac = (bfrac | FLOAT32_HIDDEN_BIT_MASK) << (32 - FLOAT32_FRACTION_SIZE); + + if (bfrac <= (afrac << 1)) { afrac >>= 1; aexp++; @@ -144,6 +151,6 @@ cfrac = (afrac << 32) / bfrac; - if (( cfrac & 0x3F ) == 0) { - cfrac |= ( bfrac * cfrac != afrac << 32 ); + if ((cfrac & 0x3F) == 0) { + cfrac |= (bfrac * cfrac != afrac << 32); } @@ -151,9 +158,9 @@ /* find first nonzero digit and shift result and detect possibly underflow */ - while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7 )))) { + while ((cexp > 0) && (cfrac) && (!(cfrac & (FLOAT32_HIDDEN_BIT_MASK << 7)))) { cexp--; cfrac <<= 1; - /* TODO: fix underflow */ - }; + /* TODO: fix underflow */ + } cfrac += (0x1 << 6); /* FIXME: 7 is not sure*/ @@ -162,8 +169,8 @@ ++cexp; cfrac >>= 1; - } + } /* check overflow */ - if (cexp >= FLOAT32_MAX_EXPONENT ) { + if (cexp >= FLOAT32_MAX_EXPONENT) { /* FIXME: overflow, return infinity */ result.parts.exp = FLOAT32_MAX_EXPONENT; @@ -181,10 +188,9 @@ cfrac >>= 1; while (cexp < 0) { - cexp ++; + cexp++; cfrac >>= 1; - } - + } } else { - result.parts.exp = (uint32_t)cexp; + result.parts.exp = (uint32_t) cexp; } @@ -194,4 +200,11 @@ } +/** + * Divide two double-precision floats. + * + * @param a Nominator. + * @param b Denominator. + * @return Result of division. + */ float64 divFloat64(float64 a, float64 b) { @@ -200,9 +213,9 @@ uint64_t afrac, bfrac, cfrac; uint64_t remlo, remhi; + uint64_t tmplo, tmphi; result.parts.sign = a.parts.sign ^ b.parts.sign; if (isFloat64NaN(a)) { - if (isFloat64SigNaN(b)) { /*FIXME: SigNaN*/ @@ -262,5 +275,4 @@ } - afrac = a.parts.fraction; aexp = a.parts.exp; @@ -275,9 +287,9 @@ return result; } + /* normalize it*/ - aexp++; - /* afrac is nonzero => it must stop */ - while (! (afrac & FLOAT64_HIDDEN_BIT_MASK) ) { + /* afrac is nonzero => it must stop */ + while (!(afrac & FLOAT64_HIDDEN_BIT_MASK)) { afrac <<= 1; aexp--; @@ -287,6 +299,6 @@ if (bexp == 0) { bexp++; - /* bfrac is nonzero => it must stop */ - while (! (bfrac & FLOAT64_HIDDEN_BIT_MASK) ) { + /* bfrac is nonzero => it must stop */ + while (!(bfrac & FLOAT64_HIDDEN_BIT_MASK)) { bfrac <<= 1; bexp--; @@ -294,8 +306,8 @@ } - afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 2 ); - bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK ) << (64 - FLOAT64_FRACTION_SIZE - 1); - - if ( bfrac <= (afrac << 1) ) { + afrac = (afrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 2); + bfrac = (bfrac | FLOAT64_HIDDEN_BIT_MASK) << (64 - FLOAT64_FRACTION_SIZE - 1); + + if (bfrac <= (afrac << 1)) { afrac >>= 1; aexp++; @@ -304,18 +316,15 @@ cexp = aexp - bexp + FLOAT64_BIAS - 2; - cfrac = divFloat64estim(afrac, bfrac); - - if (( cfrac & 0x1FF ) <= 2) { /*FIXME:?? */ - mul64integers( bfrac, cfrac, &remlo, &remhi); - /* (__u128)afrac << 64 - ( ((__u128)remhi<<64) + (__u128)remlo )*/ - remhi = afrac - remhi - ( remlo > 0); - remlo = - remlo; + cfrac = div128est(afrac, 0x0ll, bfrac); + + if ((cfrac & 0x1FF) <= 2) { + mul64(bfrac, cfrac, &tmphi, &tmplo); + sub128(afrac, 0x0ll, tmphi, tmplo, &remhi, &remlo); while ((int64_t) remhi < 0) { cfrac--; - remlo += bfrac; - remhi += ( remlo < bfrac ); - } - cfrac |= ( remlo != 0 ); + add128(remhi, remlo, 0x0ll, bfrac, &remhi, &remlo); + } + cfrac |= (remlo != 0); } @@ -323,38 +332,205 @@ result = finishFloat64(cexp, cfrac, result.parts.sign); return result; - } -uint64_t divFloat64estim(uint64_t a, uint64_t b) +/** + * Divide two quadruple-precision floats. + * + * @param a Nominator. + * @param b Denominator. + * @return Result of division. + */ +float128 divFloat128(float128 a, float128 b) { - uint64_t bhi; - uint64_t remhi, remlo; - uint64_t result; - - if ( b <= a ) { - return 0xFFFFFFFFFFFFFFFFull; - } - - bhi = b >> 32; - result = ((bhi << 32) <= a) ?( 0xFFFFFFFFull << 32) : ( a / bhi) << 32; - mul64integers(b, result, &remlo, &remhi); - - remhi = a - remhi - (remlo > 0); - remlo = - remlo; - - b <<= 32; - while ( (int64_t) remhi < 0 ) { - result -= 0x1ll << 32; - remlo += b; - remhi += bhi + ( remlo < b ); - } - remhi = (remhi << 32) | (remlo >> 32); - if (( bhi << 32) <= remhi) { - result |= 0xFFFFFFFF; - } else { - result |= remhi / bhi; - } - - + float128 result; + int64_t aexp, bexp, cexp; + uint64_t afrac_hi, afrac_lo, bfrac_hi, bfrac_lo, cfrac_hi, cfrac_lo; + uint64_t shift_out; + uint64_t rem_hihi, rem_hilo, rem_lohi, rem_lolo; + uint64_t tmp_hihi, tmp_hilo, tmp_lohi, tmp_lolo; + + result.parts.sign = a.parts.sign ^ b.parts.sign; + + if (isFloat128NaN(a)) { + if (isFloat128SigNaN(b)) { + /*FIXME: SigNaN*/ + return b; + } + + if (isFloat128SigNaN(a)) { + /*FIXME: SigNaN*/ + } + /*NaN*/ + return a; + } + + if (isFloat128NaN(b)) { + if (isFloat128SigNaN(b)) { + /*FIXME: SigNaN*/ + } + /*NaN*/ + return b; + } + + if (isFloat128Infinity(a)) { + if (isFloat128Infinity(b) || isFloat128Zero(b)) { + /*FIXME: inf / inf */ + result.binary.hi = FLOAT128_NAN_HI; + result.binary.lo = FLOAT128_NAN_LO; + return result; + } + /* inf / num */ + result.parts.exp = a.parts.exp; + result.parts.frac_hi = a.parts.frac_hi; + result.parts.frac_lo = a.parts.frac_lo; + return result; + } + + if (isFloat128Infinity(b)) { + if (isFloat128Zero(a)) { + /* FIXME 0 / inf */ + result.parts.exp = 0; + result.parts.frac_hi = 0; + result.parts.frac_lo = 0; + return result; + } + /* FIXME: num / inf*/ + result.parts.exp = 0; + result.parts.frac_hi = 0; + result.parts.frac_lo = 0; + return result; + } + + if (isFloat128Zero(b)) { + if (isFloat128Zero(a)) { + /*FIXME: 0 / 0*/ + result.binary.hi = FLOAT128_NAN_HI; + result.binary.lo = FLOAT128_NAN_LO; + return result; + } + /* FIXME: division by zero */ + result.parts.exp = 0; + result.parts.frac_hi = 0; + result.parts.frac_lo = 0; + return result; + } + + afrac_hi = a.parts.frac_hi; + afrac_lo = a.parts.frac_lo; + aexp = a.parts.exp; + bfrac_hi = b.parts.frac_hi; + bfrac_lo = b.parts.frac_lo; + bexp = b.parts.exp; + + /* denormalized numbers */ + if (aexp == 0) { + if (eq128(afrac_hi, afrac_lo, 0x0ll, 0x0ll)) { + result.parts.exp = 0; + result.parts.frac_hi = 0; + result.parts.frac_lo = 0; + return result; + } + + /* normalize it*/ + aexp++; + /* afrac is nonzero => it must stop */ + and128(afrac_hi, afrac_lo, + FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, + &tmp_hihi, &tmp_lolo); + while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) { + lshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo); + aexp--; + } + } + + if (bexp == 0) { + bexp++; + /* bfrac is nonzero => it must stop */ + and128(bfrac_hi, bfrac_lo, + FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, + &tmp_hihi, &tmp_lolo); + while (!lt128(0x0ll, 0x0ll, tmp_hihi, tmp_lolo)) { + lshift128(bfrac_hi, bfrac_lo, 1, &bfrac_hi, &bfrac_lo); + bexp--; + } + } + + or128(afrac_hi, afrac_lo, + FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, + &afrac_hi, &afrac_lo); + lshift128(afrac_hi, afrac_lo, + (128 - FLOAT128_FRACTION_SIZE - 1), &afrac_hi, &afrac_lo); + or128(bfrac_hi, bfrac_lo, + FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, + &bfrac_hi, &bfrac_lo); + lshift128(bfrac_hi, bfrac_lo, + (128 - FLOAT128_FRACTION_SIZE - 1), &bfrac_hi, &bfrac_lo); + + if (le128(bfrac_hi, bfrac_lo, afrac_hi, afrac_lo)) { + rshift128(afrac_hi, afrac_lo, 1, &afrac_hi, &afrac_lo); + aexp++; + } + + cexp = aexp - bexp + FLOAT128_BIAS - 2; + + cfrac_hi = div128est(afrac_hi, afrac_lo, bfrac_hi); + + mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_hi, + &tmp_lolo /* dummy */, &tmp_hihi, &tmp_hilo, &tmp_lohi); + + /* sub192(afrac_hi, afrac_lo, 0, + * tmp_hihi, tmp_hilo, tmp_lohi + * &rem_hihi, &rem_hilo, &rem_lohi); */ + sub128(afrac_hi, afrac_lo, tmp_hihi, tmp_hilo, &rem_hihi, &rem_hilo); + if (tmp_lohi > 0) { + sub128(rem_hihi, rem_hilo, 0x0ll, 0x1ll, &rem_hihi, &rem_hilo); + } + rem_lohi = -tmp_lohi; + + while ((int64_t) rem_hihi < 0) { + --cfrac_hi; + /* add192(rem_hihi, rem_hilo, rem_lohi, + * 0, bfrac_hi, bfrac_lo, + * &rem_hihi, &rem_hilo, &rem_lohi); */ + add128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo, &rem_hilo, &rem_lohi); + if (lt128(rem_hilo, rem_lohi, bfrac_hi, bfrac_lo)) { + ++rem_hihi; + } + } + + cfrac_lo = div128est(rem_hilo, rem_lohi, bfrac_lo); + + if ((cfrac_lo & 0x3FFF) <= 4) { + mul128(bfrac_hi, bfrac_lo, 0x0ll, cfrac_lo, + &tmp_hihi /* dummy */, &tmp_hilo, &tmp_lohi, &tmp_lolo); + + /* sub192(rem_hilo, rem_lohi, 0, + * tmp_hilo, tmp_lohi, tmp_lolo, + * &rem_hilo, &rem_lohi, &rem_lolo); */ + sub128(rem_hilo, rem_lohi, tmp_hilo, tmp_lohi, &rem_hilo, &rem_lohi); + if (tmp_lolo > 0) { + sub128(rem_hilo, rem_lohi, 0x0ll, 0x1ll, &rem_hilo, &rem_lohi); + } + rem_lolo = -tmp_lolo; + + while ((int64_t) rem_hilo < 0) { + --cfrac_lo; + /* add192(rem_hilo, rem_lohi, rem_lolo, + * 0, bfrac_hi, bfrac_lo, + * &rem_hilo, &rem_lohi, &rem_lolo); */ + add128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo, &rem_lohi, &rem_lolo); + if (lt128(rem_lohi, rem_lolo, bfrac_hi, bfrac_lo)) { + ++rem_hilo; + } + } + + cfrac_lo |= ((rem_hilo | rem_lohi | rem_lolo) != 0 ); + } + + shift_out = cfrac_lo << (64 - (128 - FLOAT128_FRACTION_SIZE - 1)); + rshift128(cfrac_hi, cfrac_lo, (128 - FLOAT128_FRACTION_SIZE - 1), + &cfrac_hi, &cfrac_lo); + + result = finishFloat128(cexp, cfrac_hi, cfrac_lo, result.parts.sign, shift_out); return result; }