Index: llvm/lib/Support/APFloat.cpp =================================================================== --- llvm/lib/Support/APFloat.cpp +++ llvm/lib/Support/APFloat.cpp @@ -1738,40 +1738,110 @@ return fs; } -/* Normalized remainder. This is not currently correct in all cases. */ +/* Normalized remainder. */ IEEEFloat::opStatus IEEEFloat::remainder(const IEEEFloat &rhs) { opStatus fs; - IEEEFloat V = *this; unsigned int origSign = sign; - fs = V.divide(rhs, rmNearestTiesToEven); - if (fs == opDivByZero) + // First handle the special cases. + fs = modSpecials(rhs); + if (fs != opOK || !isFinite()) return fs; - int parts = partCount(); - integerPart *x = new integerPart[parts]; - bool ignored; - fs = V.convertToInteger(makeMutableArrayRef(x, parts), - parts * integerPartWidth, true, rmNearestTiesToEven, - &ignored); - if (fs == opInvalidOp) { - delete[] x; - return fs; + // If we can, try to use mod, instead of duplicating the code. + // We want to make sure the current value is less than twice the denom rhs. + if (rhs.isFinite()) { + IEEEFloat P2 = rhs; + if (P2.add(rhs, rmNearestTiesToEven) == opOK) { + fs = mod(P2); + assert(fs == opOK); + } } - fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true, - rmNearestTiesToEven); - assert(fs==opOK); // should always work + // Lets work with absolute numbers. + IEEEFloat P = rhs; + P.sign = false; + sign = false; - fs = V.multiply(rhs, rmNearestTiesToEven); - assert(fs==opOK || fs==opInexact); // should not overflow or underflow + // If the value is equal to the denom, then there is no remainder. + if (compare(P) == cmpEqual) { + makeZero(origSign); + return opOK; + } - fs = subtract(V, rmNearestTiesToEven); - assert(fs==opOK || fs==opInexact); // likewise + // + // To calculate the remainder we use the following scheme. + // + // The remainder is defained as follows: + // + // remainder = numer - rquot * denom = x - r * p + // + // Where r is the result of: x/p, rounded toward the nearest integral value + // (with halfway cases rounded toward the even number). + // + // Currently, (after x mod 2p): + // r is the number of 2p's present inside x, which is inherently, an even + // number of p's. + // + // We may split the remaining calculation into 4 options: + // - if x < 0.5p then we round to the nearest number with is 0, and are done. + // - if x == 0.5p then we round to the nearest even number which is 0, and we + // are done as well. + // - if 0.5p < x < p then we round to nearest number which is 1, and we have + // to subtract 1p at least once. + // - if x >= p then we must subtract p at least once, as x must be a + // remainder. + // + // By now, we were done, or we added 1 to r, which in turn, now an odd number. + // + // We can now split the remaining calculation to the following 3 options: + // - if x < 0.5p then we round to the nearest number with is 0, and are done. + // - if x == 0.5p then we round to the nearest even number. As r is odd, we + // must round up to the next even number. so we must subtract p once more. + // - if x > 0.5p (and inherently x < p) then we must round r up to the next + // integral, and subtract p once more. + // + + const IEEEFloat two(rhs.getSemantics(), 2); + IEEEFloat P2 = P; + if (P2.divide(two, rmNearestTiesToEven) == opOK) { + if (compare(P2) == cmpGreaterThan) { + fs = subtract(P, rmNearestTiesToEven); + assert(fs == opOK); + + cmpResult result = compare(P2); + if (result == cmpGreaterThan || result == cmpEqual) { + fs = subtract(P, rmNearestTiesToEven); + assert(fs == opOK); + } + } + } else { + // In case the division did not work, we will have to simulate it with twice + // the size of x, instead. + IEEEFloat V2 = *this; + fs = V2.add(*this, rmNearestTiesToEven); + assert(fs == opOK); + + if (V2.compare(P2) == cmpGreaterThan) { + fs = subtract(P, rmNearestTiesToEven); + assert(fs == opOK); + + V2 = *this; + fs = V2.add(*this, rmNearestTiesToEven); + assert(fs == opOK); + + cmpResult result = V2.compare(P2); + if (result == cmpGreaterThan || result == cmpEqual) { + fs = subtract(P, rmNearestTiesToEven); + assert(fs == opOK); + } + } + } if (isZero()) sign = origSign; // IEEE754 requires this - delete[] x; + else + sign ^= origSign; return fs; } Index: llvm/unittests/ADT/APFloatTest.cpp =================================================================== --- llvm/unittests/ADT/APFloatTest.cpp +++ llvm/unittests/ADT/APFloatTest.cpp @@ -3482,6 +3482,83 @@ } } +TEST(APFloatTest, remainder) { + { + APFloat f1(APFloat::IEEEdouble(), "1.5"); + APFloat f2(APFloat::IEEEdouble(), "1.0"); + APFloat expected(APFloat::IEEEdouble(), "-0.5"); + EXPECT_EQ(APFloat::opOK, f1.remainder(f2)); + EXPECT_TRUE(f1.bitwiseIsEqual(expected)); + } + { + APFloat f1(APFloat::IEEEdouble(), "0.5"); + APFloat f2(APFloat::IEEEdouble(), "1.0"); + APFloat expected(APFloat::IEEEdouble(), "0.5"); + EXPECT_EQ(APFloat::opOK, f1.remainder(f2)); + EXPECT_TRUE(f1.bitwiseIsEqual(expected)); + } + { + APFloat f1(APFloat::IEEEdouble(), "0x1.3333333333333p-2"); // 0.3 + APFloat f2(APFloat::IEEEdouble(), "0x1.47ae147ae147bp-7"); // 0.01 + APFloat expected(APFloat::IEEEdouble(), "-0x1.4p-56"); + EXPECT_EQ(APFloat::opOK, f1.remainder(f2)); + EXPECT_TRUE(f1.bitwiseIsEqual(expected)); + } + { + APFloat f1(APFloat::IEEEdouble(), "0x1p64"); // 1.8446744073709552e19 + APFloat f2(APFloat::IEEEdouble(), "1.5"); + APFloat expected(APFloat::IEEEdouble(), "-0.5"); + EXPECT_EQ(APFloat::opOK, f1.remainder(f2)); + EXPECT_TRUE(f1.bitwiseIsEqual(expected)); + } + { + APFloat f1(APFloat::IEEEdouble(), "0x1p1000"); + APFloat f2(APFloat::IEEEdouble(), "0x1p-1000"); + APFloat expected(APFloat::IEEEdouble(), "0.0"); + EXPECT_EQ(APFloat::opOK, f1.remainder(f2)); + EXPECT_TRUE(f1.bitwiseIsEqual(expected)); + } + { + APFloat f1(APFloat::IEEEdouble(), "0.0"); + APFloat f2(APFloat::IEEEdouble(), "1.0"); + APFloat expected(APFloat::IEEEdouble(), "0.0"); + EXPECT_EQ(APFloat::opOK, f1.remainder(f2)); + EXPECT_TRUE(f1.bitwiseIsEqual(expected)); + } + { + APFloat f1(APFloat::IEEEdouble(), "1.0"); + APFloat f2(APFloat::IEEEdouble(), "0.0"); + EXPECT_EQ(f1.remainder(f2), APFloat::opInvalidOp); + EXPECT_TRUE(f1.isNaN()); + } + { + APFloat f1(APFloat::IEEEdouble(), "0.0"); + APFloat f2(APFloat::IEEEdouble(), "0.0"); + EXPECT_EQ(f1.remainder(f2), APFloat::opInvalidOp); + EXPECT_TRUE(f1.isNaN()); + } + { + APFloat f1 = APFloat::getInf(APFloat::IEEEdouble(), false); + APFloat f2(APFloat::IEEEdouble(), "1.0"); + EXPECT_EQ(f1.remainder(f2), APFloat::opInvalidOp); + EXPECT_TRUE(f1.isNaN()); + } + { + APFloat f1(APFloat::IEEEdouble(), "-4.0"); + APFloat f2(APFloat::IEEEdouble(), "-2.0"); + APFloat expected(APFloat::IEEEdouble(), "-0.0"); + EXPECT_EQ(APFloat::opOK, f1.remainder(f2)); + EXPECT_TRUE(f1.bitwiseIsEqual(expected)); + } + { + APFloat f1(APFloat::IEEEdouble(), "-4.0"); + APFloat f2(APFloat::IEEEdouble(), "2.0"); + APFloat expected(APFloat::IEEEdouble(), "-0.0"); + EXPECT_EQ(APFloat::opOK, f1.remainder(f2)); + EXPECT_TRUE(f1.bitwiseIsEqual(expected)); + } +} + TEST(APFloatTest, PPCDoubleDoubleAddSpecial) { using DataType = std::tuple;