Fix APFloat mod

[simonbyrne]

Summary

The existing code for APFloat mod was incorrect in two ways:

• intermediate rounding gave incorrect results in some cases (e.g. frem double 0.3 0.01)
• the method of integer truncation could cause intermediate overflow

This caused subsequent errors in the constant propagation pass, see https://llvm.org/bugs/show_bug.cgi?id=3316

[simonbyrne]

Use fusedMultiplyAdd.

[scanon]

Numerically this is sound, and something that I've wanted to fix for a while but haven't had time to deal with. Thanks for doing it! Please add some test cases that would detect the old error (@gottesmm can likely point you in the right direction).

[simonbyrne]

Thanks. I'm not really that familiar with the layout of the repository: where should tests go?

[Meinersbur]

In D29346#664495, @simonbyrne wrote:

Thanks. I'm not really that familiar with the layout of the repository: where should tests go?

/llvm/unittests/ADT/APFloatTest.cpp

[simonbyrne]

Changed to use subtraction to avoid intermediate overflow, added tests.

[scanon]

Great, thanks Simon! This looks fine numerically. @gottesmm, can you double-check for style, etc?

[gottesmm]

Did you run this through git-clang-format?

[simonbyrne]

Updated to use git clang-format

[gottesmm]

Thanks for running it through git-clang-format. LGTM!

[skatkov]

This patch causes incorrect functional behavior, specifically, IEEE specification requires the sign of the remainder is the same as numerator in case remainder is zero.
It is incorrect with this patch. The following tests can be used to show the issue:

{
  APFloat f1(APFloat::IEEEdouble(), "-4.0");
  APFloat f2(APFloat::IEEEdouble(), "-2.0");
  APFloat expected(APFloat::IEEEdouble(), "-0.0");
  EXPECT_EQ(f1.mod(f2), APFloat::opOK);
  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(f1.mod(f2), APFloat::opOK);
  EXPECT_TRUE(f1.bitwiseIsEqual(expected));
}

I plan to revert this patch and add mentioned tests...

[skatkov]

Actually, I've upload a fix, please review it.

llvm/trunk/lib/Support/APFloat.cpp

IEEEFloat::opStatus IEEEFloat::remainder(const IEEEFloat &rhs) {
   assert(fs==opOK || fs==opInexact);   // likewise
 
   if (isZero())
     sign = origSign;    // IEEE754 requires this
   delete[] x;
   return fs;
 }
 
-/* Normalized llvm frem (C fmod).
-   This is not currently correct in all cases.  */
+/* Normalized llvm frem (C fmod). */
 IEEEFloat::opStatus IEEEFloat::mod(const IEEEFloat &rhs) {
   opStatus fs;
   fs = modSpecials(rhs);
 
-  if (isFiniteNonZero() && rhs.isFiniteNonZero()) {
-    IEEEFloat V = *this;
-    unsigned int origSign = sign;
-
-    fs = V.divide(rhs, rmNearestTiesToEven);
-    if (fs == opDivByZero)
-      return fs;
-
-    int parts = partCount();
-    integerPart *x = new integerPart[parts];
-    bool ignored;
-    fs = V.convertToInteger(makeMutableArrayRef(x, parts),
-                            parts * integerPartWidth, true, rmTowardZero,
-                            &ignored);
-    if (fs == opInvalidOp) {
-      delete[] x;
-      return fs;
-    }
-
-    fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
-                                          rmNearestTiesToEven);
-    assert(fs==opOK);   // should always work
-
-    fs = V.multiply(rhs, rmNearestTiesToEven);
-    assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
+  while (isFiniteNonZero() && rhs.isFiniteNonZero() &&
+         compareAbsoluteValue(rhs) != cmpLessThan) {
+    IEEEFloat V = scalbn(rhs, ilogb(*this) - ilogb(rhs), rmNearestTiesToEven);
+    if (compareAbsoluteValue(V) == cmpLessThan)
+      V = scalbn(V, -1, rmNearestTiesToEven);
+    V.sign = sign;
   
     fs = subtract(V, rmNearestTiesToEven);
-    assert(fs==opOK || fs==opInexact);   // likewise
-
-    if (isZero())
-      sign = origSign;    // IEEE754 requires this
-    delete[] x;
+    assert(fs==opOK);
   }
   return fs;
 }
 
 /* Normalized fused-multiply-add.  */
 IEEEFloat::opStatus IEEEFloat::fusedMultiplyAdd(const IEEEFloat &multiplicand,
                                                 const IEEEFloat &addend,
                                                 roundingMode rounding_mode) {

llvm/trunk/unittests/ADT/APFloatTest.cpp

TEST(APFloatTest, frexp) {
   EXPECT_EQ(-50, Exp);
   EXPECT_TRUE(APFloat(APFloat::IEEEdouble(), "0x1p-1").bitwiseIsEqual(Frac));
 
   Frac = frexp(APFloat(APFloat::IEEEdouble(), "0x1.c60f120d9f87cp+51"), Exp, RM);
   EXPECT_EQ(52, Exp);
   EXPECT_TRUE(APFloat(APFloat::IEEEdouble(), "0x1.c60f120d9f87cp-1").bitwiseIsEqual(Frac));
 }
 
+TEST(APFloatTest, mod) {
+  {
+    APFloat f1(APFloat::IEEEdouble(), "1.5");
+    APFloat f2(APFloat::IEEEdouble(), "1.0");
+    APFloat expected(APFloat::IEEEdouble(), "0.5");
+    EXPECT_EQ(f1.mod(f2), APFloat::opOK);
+    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(f1.mod(f2), APFloat::opOK);
+    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.47ae147ae1471p-7"); // 0.009999999999999983
+    EXPECT_EQ(f1.mod(f2), APFloat::opOK);
+    EXPECT_TRUE(f1.bitwiseIsEqual(expected));
+  }
+  {
+    APFloat f1(APFloat::IEEEdouble(), "0x1p64"); // 1.8446744073709552e19
+    APFloat f2(APFloat::IEEEdouble(), "1.5");
+    APFloat expected(APFloat::IEEEdouble(), "1.0");
+    EXPECT_EQ(f1.mod(f2), APFloat::opOK);
+    EXPECT_TRUE(f1.bitwiseIsEqual(expected));
+  }
+  {
+    APFloat f1(APFloat::IEEEdouble(), "0x1p1000");
+    APFloat f2(APFloat::IEEEdouble(), "0x1p-1000");
+    APFloat expected(APFloat::IEEEdouble(), "0.0");
+    EXPECT_EQ(f1.mod(f2), APFloat::opOK);
+    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(f1.mod(f2), APFloat::opOK);
+    EXPECT_TRUE(f1.bitwiseIsEqual(expected));
+  }
+  {
+    APFloat f1(APFloat::IEEEdouble(), "1.0");
+    APFloat f2(APFloat::IEEEdouble(), "0.0");
+    EXPECT_EQ(f1.mod(f2), APFloat::opInvalidOp);
+    EXPECT_TRUE(f1.isNaN());
+  }
+  {
+    APFloat f1(APFloat::IEEEdouble(), "0.0");
+    APFloat f2(APFloat::IEEEdouble(), "0.0");
+    EXPECT_EQ(f1.mod(f2), APFloat::opInvalidOp);
+    EXPECT_TRUE(f1.isNaN());
+  }
+  {
+    APFloat f1 = APFloat::getInf(APFloat::IEEEdouble(), false);
+    APFloat f2(APFloat::IEEEdouble(), "1.0");
+    EXPECT_EQ(f1.mod(f2), APFloat::opInvalidOp);
+    EXPECT_TRUE(f1.isNaN());
+  }
+}
+
 TEST(APFloatTest, PPCDoubleDoubleAddSpecial) {
   using DataType = std::tuple<uint64_t, uint64_t, uint64_t, uint64_t,
-                              APFloat::fltCategory, APFloat::roundingMode>;
-  DataType Data[] = {
+                              APFloat::fltCategory, APFloat::roundingMode>;  DataType Data[] = {
       // (1 + 0) + (-1 + 0) = fcZero
       std::make_tuple(0x3ff0000000000000ull, 0, 0xbff0000000000000ull, 0,
                       APFloat::fcZero, APFloat::rmNearestTiesToEven),
       // LDBL_MAX + (1.1 >> (1023 - 106) + 0)) = fcInfinity
       std::make_tuple(0x7fefffffffffffffull, 0x7c8ffffffffffffeull,
                       0x7948000000000000ull, 0ull, APFloat::fcInfinity,
                       APFloat::rmNearestTiesToEven),
       // TODO: change the 4th 0x75effffffffffffe to 0x75efffffffffffff when