[APFloat] Fix FP remainder operation

[ekatz]

Summary

Fix IEEEFloat::remainder() function.

Fix PR3359.

[ekatz]

Rebased to master.

[ekatz]

Ping

[ekatz]

Ping

[hfinkel]

Thanks for posting patches to fix problems with the APFloat implementation. FWIW, this code is hard to review, from my perspective, because it lacks comments (both your new code and the code it is replacing). As you clearly have gone through this and understand what it's doing, could you add a bunch of comments explaining what's going on? I think that would be very helpful for everyone else.

[arsenm]

llvm/unittests/ADT/APFloatTest.cpp

3379

These should be swapped due to a gtestism (same for the rest, it should be expected, actual)

[ekatz]

Added comments and fix tests to meet gtestism.

[ekatz]

Ping

[ekatz]

Ping

[efriedma]

I'd like to see more test coverage for denormal denominators, and denominators close to the max finite float.

llvm/lib/Support/APFloat.cpp

1755

Can we return early if rhs.isFinite() is false? The following code might work for rhs==infinity, but it makes everything more complicated to reason about.

1757

If the add is not opOK, that means the add returned opOverflow, so the value is already less than twice rhs? I guess that's right, but it would be nice to spell it out.

1774

Is the special case for compare(P) == cmpEqual necessary?

1812

If this check fails, that means P2 is denormal?

Can you make this work for denormal cases by just fiddling with the exponent of P2, instead of doing a proper divide? The resulting APFloat might not be properly normalized, but that doesn't matter here.

[ekatz]

Simplified logics according to the notes.
Also added more extensive tests.

[ekatz]

@efriedma , I updated the files according to your requests.

[efriedma]

llvm/lib/Support/APFloat.cpp

1751

The opDivByZero thing is sort of confusing. Maybe just something like if (!isFinite() || !rhs.isFinite()) return remainderSpecials(rhs);?

1785

Are all the subtractions here guaranteed to produce exact results? If they are, would it make sense to assert that?

[ekatz]

llvm/lib/Support/APFloat.cpp

1751

I thought about that, but I just used the same convention as in addOrSubtractSpecials.

If we were to add this if statement, then we are actually handling specials outside the function that supposed to do exactly that. And that alone makes even less sense than the strange opStatus (which is commented, and used only internally).

It would have made more sense to add a new opStatus just for that case, but it seems a bit too much, I guess (at least for those internal functions).

Which is why I decided to take this approach.

[ekatz]

Added asserts for the opStatus returned by the operations.

[efriedma]

LGTM

llvm/lib/Support/APFloat.cpp

if (isFiniteNonZero()) {
     fs = normalize(rounding_mode, lost_fraction);
     if (lost_fraction != lfExactlyZero)
       fs = (opStatus) (fs | opInexact);
   }
 
   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)
+  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 (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
+  IEEEFloat P = rhs;
+  P.sign = false;
+  sign = false;
 
-  fs = V.multiply(rhs, rmNearestTiesToEven);
-  assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
+  if (compare(P) == cmpEqual) {
+    makeZero(origSign);
+    return opOK;
+  }
 
-  fs = subtract(V, rmNearestTiesToEven);
-  assert(fs==opOK || fs==opInexact);   // likewise
+  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 {
+    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;
 }
 
 /* Normalized llvm frem (C fmod). */
 IEEEFloat::opStatus IEEEFloat::mod(const IEEEFloat &rhs) {
   opStatus fs;
   fs = modSpecials(rhs);
   unsigned int origSign = sign;
 
   while (isFiniteNonZero() && rhs.isFiniteNonZero() &&
          compareAbsoluteValue(rhs) != cmpLessThan) {
     IEEEFloat V = scalbn(rhs, ilogb(*this) - ilogb(rhs), rmNearestTiesToEven);
     if (compareAbsoluteValue(V) == cmpLessThan)

llvm/unittests/ADT/APFloatTest.cpp

TEST(APFloatTest, mod) {
     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));
   }
 }
 
+TEST(APFloatTest, remainder) {
+  {
+    APFloat f1(APFloat::IEEEdouble(), "1.5");
+    APFloat f2(APFloat::IEEEdouble(), "1.0");
+    APFloat expected(APFloat::IEEEdouble(), "-0.5");
+    EXPECT_EQ(f1.remainder(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.remainder(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.4p-56");
+    EXPECT_EQ(f1.remainder(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(), "-0.5");
+    EXPECT_EQ(f1.remainder(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.remainder(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.remainder(f2), APFloat::opOK);
+    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(f1.remainder(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.remainder(f2), APFloat::opOK);
+    EXPECT_TRUE(f1.bitwiseIsEqual(expected));
+  }
+}
+
 TEST(APFloatTest, PPCDoubleDoubleAddSpecial) {
   using DataType = std::tuple<uint64_t, uint64_t, uint64_t, uint64_t,
                               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