Index: uspace/lib/softfloat/generic/sub.c =================================================================== --- uspace/lib/softfloat/generic/sub.c (revision 750636ac4a65360c44846b7681af6ae46854dd3a) +++ uspace/lib/softfloat/generic/sub.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 Substraction functions. */ @@ -36,6 +37,12 @@ #include #include - -/** Subtract two float32 numbers with same signs +#include + +/** + * Subtract two single-precision floats with the same signs. + * + * @param a First input operand. + * @param b Second input operand. + * @return Result of substraction. */ float32 subFloat32(float32 a, float32 b) @@ -52,7 +59,7 @@ /* TODO: fix SigNaN */ if (isFloat32SigNaN(b)) { - }; - return b; - }; + } + return b; + } if (b.parts.exp == FLOAT32_MAX_EXPONENT) { @@ -72,7 +79,7 @@ /* TODO: fix SigNaN */ if (isFloat32SigNaN(a) || isFloat32SigNaN(b)) { - }; - return a; - }; + } + return a; + } if (a.parts.exp == FLOAT32_MAX_EXPONENT) { @@ -82,5 +89,5 @@ result.binary = FLOAT32_NAN; return result; - }; + } return a; } @@ -92,16 +99,16 @@ frac2 = b.parts.fraction; exp2 = b.parts.exp; - }; + } if (exp1 == 0) { /* both are denormalized */ - result.parts.fraction = frac1-frac2; + result.parts.fraction = frac1 - frac2; if (result.parts.fraction > frac1) { /* TODO: underflow exception */ return result; - }; + } result.parts.exp = 0; return result; - }; + } /* add hidden bit */ @@ -114,5 +121,5 @@ /* normalized */ frac2 |= FLOAT32_HIDDEN_BIT_MASK; - }; + } /* create some space for rounding */ @@ -121,8 +128,9 @@ if (expdiff > FLOAT32_FRACTION_SIZE + 1) { - goto done; - }; + goto done; + } frac1 = frac1 - (frac2 >> expdiff); + done: /* TODO: find first nonzero digit and shift result and detect possibly underflow */ @@ -130,6 +138,6 @@ --exp1; frac1 <<= 1; - /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */ - }; + /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */ + } /* rounding - if first bit after fraction is set then round up */ @@ -139,7 +147,7 @@ ++exp1; frac1 >>= 1; - }; - - /*Clear hidden bit and shift */ + } + + /* Clear hidden bit and shift */ result.parts.fraction = ((frac1 >> 6) & (~FLOAT32_HIDDEN_BIT_MASK)); result.parts.exp = exp1; @@ -148,5 +156,10 @@ } -/** Subtract two float64 numbers with same signs +/** + * Subtract two double-precision floats with the same signs. + * + * @param a First input operand. + * @param b Second input operand. + * @return Result of substraction. */ float64 subFloat64(float64 a, float64 b) @@ -164,7 +177,7 @@ /* TODO: fix SigNaN */ if (isFloat64SigNaN(b)) { - }; - return b; - }; + } + return b; + } if (b.parts.exp == FLOAT64_MAX_EXPONENT) { @@ -184,7 +197,7 @@ /* TODO: fix SigNaN */ if (isFloat64SigNaN(a) || isFloat64SigNaN(b)) { - }; - return a; - }; + } + return a; + } if (a.parts.exp == FLOAT64_MAX_EXPONENT) { @@ -194,5 +207,5 @@ result.binary = FLOAT64_NAN; return result; - }; + } return a; } @@ -204,5 +217,5 @@ frac2 = b.parts.fraction; exp2 = b.parts.exp; - }; + } if (exp1 == 0) { @@ -212,8 +225,8 @@ /* TODO: underflow exception */ return result; - }; + } result.parts.exp = 0; return result; - }; + } /* add hidden bit */ @@ -226,5 +239,5 @@ /* normalized */ frac2 |= FLOAT64_HIDDEN_BIT_MASK; - }; + } /* create some space for rounding */ @@ -233,8 +246,9 @@ if (expdiff > FLOAT64_FRACTION_SIZE + 1) { - goto done; - }; + goto done; + } frac1 = frac1 - (frac2 >> expdiff); + done: /* TODO: find first nonzero digit and shift result and detect possibly underflow */ @@ -242,6 +256,6 @@ --exp1; frac1 <<= 1; - /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */ - }; + /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */ + } /* rounding - if first bit after fraction is set then round up */ @@ -251,7 +265,7 @@ ++exp1; frac1 >>= 1; - }; - - /*Clear hidden bit and shift */ + } + + /* Clear hidden bit and shift */ result.parts.fraction = ((frac1 >> 6) & (~FLOAT64_HIDDEN_BIT_MASK)); result.parts.exp = exp1; @@ -260,4 +274,152 @@ } +/** + * Subtract two quadruple-precision floats with the same signs. + * + * @param a First input operand. + * @param b Second input operand. + * @return Result of substraction. + */ +float128 subFloat128(float128 a, float128 b) +{ + int expdiff; + uint32_t exp1, exp2; + uint64_t frac1_hi, frac1_lo, frac2_hi, frac2_lo, tmp_hi, tmp_lo; + float128 result; + + result.binary.hi = 0; + result.binary.lo = 0; + + expdiff = a.parts.exp - b.parts.exp; + if ((expdiff < 0 ) || ((expdiff == 0) && + lt128(a.parts.frac_hi, a.parts.frac_lo, b.parts.frac_hi, b.parts.frac_lo))) { + if (isFloat128NaN(b)) { + /* TODO: fix SigNaN */ + if (isFloat128SigNaN(b)) { + } + return b; + } + + if (b.parts.exp == FLOAT128_MAX_EXPONENT) { + b.parts.sign = !b.parts.sign; /* num -(+-inf) = -+inf */ + return b; + } + + result.parts.sign = !a.parts.sign; + + frac1_hi = b.parts.frac_hi; + frac1_lo = b.parts.frac_lo; + exp1 = b.parts.exp; + frac2_hi = a.parts.frac_hi; + frac2_lo = a.parts.frac_lo; + exp2 = a.parts.exp; + expdiff *= -1; + } else { + if (isFloat128NaN(a)) { + /* TODO: fix SigNaN */ + if (isFloat128SigNaN(a) || isFloat128SigNaN(b)) { + } + return a; + } + + if (a.parts.exp == FLOAT128_MAX_EXPONENT) { + if (b.parts.exp == FLOAT128_MAX_EXPONENT) { + /* inf - inf => nan */ + /* TODO: fix exception */ + result.binary.hi = FLOAT128_NAN_HI; + result.binary.lo = FLOAT128_NAN_LO; + return result; + } + return a; + } + + result.parts.sign = a.parts.sign; + + frac1_hi = a.parts.frac_hi; + frac1_lo = a.parts.frac_lo; + exp1 = a.parts.exp; + frac2_hi = b.parts.frac_hi; + frac2_lo = b.parts.frac_lo; + exp2 = b.parts.exp; + } + + if (exp1 == 0) { + /* both are denormalized */ + sub128(frac1_hi, frac1_lo, frac2_hi, frac2_lo, &tmp_hi, &tmp_lo); + result.parts.frac_hi = tmp_hi; + result.parts.frac_lo = tmp_lo; + if (lt128(frac1_hi, frac1_lo, result.parts.frac_hi, result.parts.frac_lo)) { + /* TODO: underflow exception */ + return result; + } + result.parts.exp = 0; + return result; + } + + /* add hidden bit */ + or128(frac1_hi, frac1_lo, + FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, + &frac1_hi, &frac1_lo); + + if (exp2 == 0) { + /* denormalized */ + --expdiff; + } else { + /* normalized */ + or128(frac2_hi, frac2_lo, + FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, + &frac2_hi, &frac2_lo); + } + + /* create some space for rounding */ + lshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo); + lshift128(frac2_hi, frac2_lo, 6, &frac2_hi, &frac2_lo); + + if (expdiff > FLOAT128_FRACTION_SIZE + 1) { + goto done; + } + + rshift128(frac2_hi, frac2_lo, expdiff, &tmp_hi, &tmp_lo); + sub128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &frac1_hi, &frac1_lo); + +done: + /* TODO: find first nonzero digit and shift result and detect possibly underflow */ + lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 6, + &tmp_hi, &tmp_lo); + and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo); + while ((exp1 > 0) && (!lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo))) { + --exp1; + lshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo); + /* TODO: fix underflow - frac1 == 0 does not necessary means underflow... */ + + lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 6, + &tmp_hi, &tmp_lo); + and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo); + } + + /* rounding - if first bit after fraction is set then round up */ + add128(frac1_hi, frac1_lo, 0x0ll, 0x20ll, &frac1_hi, &frac1_lo); + + lshift128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, 7, + &tmp_hi, &tmp_lo); + and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo); + if (lt128(0x0ll, 0x0ll, tmp_hi, tmp_lo)) { + ++exp1; + rshift128(frac1_hi, frac1_lo, 1, &frac1_hi, &frac1_lo); + } + + /* Clear hidden bit and shift */ + rshift128(frac1_hi, frac1_lo, 6, &frac1_hi, &frac1_lo); + not128(FLOAT128_HIDDEN_BIT_MASK_HI, FLOAT128_HIDDEN_BIT_MASK_LO, + &tmp_hi, &tmp_lo); + and128(frac1_hi, frac1_lo, tmp_hi, tmp_lo, &tmp_hi, &tmp_lo); + result.parts.frac_hi = tmp_hi; + result.parts.frac_lo = tmp_lo; + + result.parts.exp = exp1; + + return result; +} + /** @} */