[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

3490

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

1753

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.

1755

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.

1770

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

1807

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

1748

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

1780

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

1748

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)
+  // 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;
 }
 
 /* Normalized llvm frem (C fmod). */
 IEEEFloat::opStatus IEEEFloat::mod(const IEEEFloat &rhs) {
   opStatus fs;
   fs = modSpecials(rhs);
   unsigned int origSign = sign;

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(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<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