[APFloat] Fix fusedMultiplyAdd when `this` equals to `Addend`

[ekatz]

Summary

The arguments to fusedMultiplyAdd are passed by reference. We must save the Addend value on the beginning of the function, before we modify this, as they may be the same reference.

Fix PR44051.

[hfinkel]

LGTM

[foad]

Typo in commit message: "this this"

llvm/lib/Support/APFloat.cpp

1003

Marginally more efficient to omit the , 0 here, to use the one-argument constructor.

Since we are unconditionally constructing the addend here anyway, would it make more sense to force all callers to pass in an addend (using zero if they're not interested in it)?

1039

Is it safe to skip this code when the addend is zero? If the result of the multiply is also zero, adding zero could affect the sign of the result, couldn't it?

[ekatz]

llvm/lib/Support/APFloat.cpp

1003

Correct. I'll change the constructor to the single parameter version.

Regarding the construction of the addend, that is a good idea.
I suggest making addend passed by value, with default to zero (keep in mind that this function is private, so there is no API change).
What do you think?

1039

Good observation, but adding/subtracting a negative zero, doesn't affect the sign (it is considered the same as positive zero).

[foad]

llvm/lib/Support/APFloat.cpp

1003

Sounds great if you can make it work. I wasn't sure exactly how to implement the "default to zero" part.

1039

A couple of cases where adding zero changes the sign, accoerding to the comment at the end of IEEEFloat::addOrSubtract:
-0 + +0 = +0 (except if rounding to -inf)
+0 + -0 = -0 (if rounding to -inf)
(I realise this code is pre-existing, so no need to address it as part of the current patch.)

[ekatz]

llvm/lib/Support/APFloat.cpp

1039

Good catch! Those corner cases never end...
I'll open a new bug for it.

[ekatz]

Pass added parameter by value (instead of - by ref), and have a default value of zero.

[ekatz]

Note that this is a safe overload, as multiplySignificand is private.

llvm/include/llvm/ADT/APFloat.h

private:
   /// @}
 
   /// \name Significand operations.
   /// @{
 
   integerPart addSignificand(const IEEEFloat &);
   integerPart subtractSignificand(const IEEEFloat &, integerPart);
   lostFraction addOrSubtractSignificand(const IEEEFloat &, bool subtract);
-  lostFraction multiplySignificand(const IEEEFloat &, const IEEEFloat *);
+  lostFraction multiplySignificand(const IEEEFloat &, IEEEFloat);
+  lostFraction multiplySignificand(const IEEEFloat&);
   lostFraction divideSignificand(const IEEEFloat &);
   void incrementSignificand();
   void initialize(const fltSemantics *);
   void shiftSignificandLeft(unsigned int);
   lostFraction shiftSignificandRight(unsigned int);
   unsigned int significandLSB() const;
   unsigned int significandMSB() const;
   void zeroSignificand();

llvm/lib/Support/APFloat.cpp

IEEEFloat::integerPart IEEEFloat::subtractSignificand(const IEEEFloat &rhs,
   return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
                            partCount());
 }
 
 /* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
    on to the full-precision result of the multiplication.  Returns the
    lost fraction.  */
 lostFraction IEEEFloat::multiplySignificand(const IEEEFloat &rhs,
-                                            const IEEEFloat *addend) {
+                                            IEEEFloat addend) {
   unsigned int omsb;        // One, not zero, based MSB.
   unsigned int partsCount, newPartsCount, precision;
   integerPart *lhsSignificand;
   integerPart scratch[4];
   integerPart *fullSignificand;
   lostFraction lost_fraction;
   bool ignored;
 
   assert(semantics == rhs.semantics);
 
   precision = semantics->precision;
 
   // Allocate space for twice as many bits as the original significand, plus one
   // extra bit for the addition to overflow into.
   newPartsCount = partCountForBits(precision * 2 + 1);
 
lostFraction IEEEFloat::multiplySignificand(const IEEEFloat &rhs,
   // the result of multiplication is:
   //   *this = c48 c47 c46 . c45 ... c0 * 2^(e1+e2)
   // Note that there are three significant bits at the left-hand side of the
   // radix point: two for the multiplication, and an overflow bit for the
   // addition (that will always be zero at this point). Move the radix point
   // toward left by two bits, and adjust exponent accordingly.
   exponent += 2;
 
-  if (addend && addend->isNonZero()) {
+  if (addend.isNonZero()) {
     // The intermediate result of the multiplication has "2 * precision"
     // signicant bit; adjust the addend to be consistent with mul result.
     //
     Significand savedSignificand = significand;
     const fltSemantics *savedSemantics = semantics;
     fltSemantics extendedSemantics;
     opStatus status;
     unsigned int extendedPrecision;
if (addend.isNonZero()) {
     extendedSemantics.precision = extendedPrecision;
 
     if (newPartsCount == 1)
       significand.part = fullSignificand[0];
     else
       significand.parts = fullSignificand;
     semantics = &extendedSemantics;
 
-    IEEEFloat extendedAddend(*addend);
-    status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored);
+    status = addend.convert(extendedSemantics, rmTowardZero, &ignored);
     assert(status == opOK);
     (void)status;
 
     // Shift the significand of the addend right by one bit. This guarantees
     // that the high bit of the significand is zero (same as fullSignificand),
     // so the addition will overflow (if it does overflow at all) into the top bit.
-    lost_fraction = extendedAddend.shiftSignificandRight(1);
+    lost_fraction = addend.shiftSignificandRight(1);
     assert(lost_fraction == lfExactlyZero &&
            "Lost precision while shifting addend for fused-multiply-add.");
 
-    lost_fraction = addOrSubtractSignificand(extendedAddend, false);
+    lost_fraction = addOrSubtractSignificand(addend, false);
 
     /* Restore our state.  */
     if (newPartsCount == 1)
       fullSignificand[0] = significand.part;
     significand = savedSignificand;
     semantics = savedSemantics;
 
     omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
lostFraction IEEEFloat::multiplySignificand(const IEEEFloat &rhs,
   APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
 
   if (newPartsCount > 4)
     delete [] fullSignificand;
 
   return lost_fraction;
 }
 
+lostFraction IEEEFloat::multiplySignificand(const IEEEFloat &rhs) {
+  return multiplySignificand(rhs, IEEEFloat(*semantics));
+}
+
 /* Multiply the significands of LHS and RHS to DST.  */
 lostFraction IEEEFloat::divideSignificand(const IEEEFloat &rhs) {
   unsigned int bit, i, partsCount;
   const integerPart *rhsSignificand;
   integerPart *lhsSignificand, *dividend, *divisor;
   integerPart scratch[4];
   lostFraction lost_fraction;
 

 IEEEFloat::opStatus IEEEFloat::multiply(const IEEEFloat &rhs,
                                         roundingMode rounding_mode) {
   opStatus fs;
 
   sign ^= rhs.sign;
   fs = multiplySpecials(rhs);
 
   if (isFiniteNonZero()) {
-    lostFraction lost_fraction = multiplySignificand(rhs, nullptr);
+    lostFraction lost_fraction = multiplySignificand(rhs);
     fs = normalize(rounding_mode, lost_fraction);
     if (lost_fraction != lfExactlyZero)
       fs = (opStatus) (fs | opInexact);
   }
 
   return fs;
 }
 
IEEEFloat::opStatus IEEEFloat::fusedMultiplyAdd(const IEEEFloat &multiplicand,
 
   /* If and only if all arguments are normal do we need to do an
      extended-precision calculation.  */
   if (isFiniteNonZero() &&
       multiplicand.isFiniteNonZero() &&
       addend.isFinite()) {
     lostFraction lost_fraction;
 
-    lost_fraction = multiplySignificand(multiplicand, &addend);
+    lost_fraction = multiplySignificand(multiplicand, addend);
     fs = normalize(rounding_mode, lost_fraction);
     if (lost_fraction != lfExactlyZero)
       fs = (opStatus) (fs | opInexact);
 
     /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
        positive zero unless rounding to minus infinity, except that
        adding two like-signed zeroes gives that zero.  */
     if (category == fcZero && !(fs & opUnderflow) && sign != addend.sign)
for (;; parts *= 2) {
     decSig.exponent += exp;
 
     lostFraction calcLostFraction;
     integerPart HUerr, HUdistance;
     unsigned int powHUerr;
 
     if (exp >= 0) {
       /* multiplySignificand leaves the precision-th bit set to 1.  */
-      calcLostFraction = decSig.multiplySignificand(pow5, nullptr);
+      calcLostFraction = decSig.multiplySignificand(pow5);
       powHUerr = powStatus != opOK;
     } else {
       calcLostFraction = decSig.divideSignificand(pow5);
       /* Denormal numbers have less precision.  */
       if (decSig.exponent < semantics->minExponent) {
         excessPrecision += (semantics->minExponent - decSig.exponent);
         truncatedBits = excessPrecision;
         if (excessPrecision > calcSemantics.precision)

llvm/unittests/ADT/APFloatTest.cpp

TEST(APFloatTest, FMA) {
   // Regression test that failed an assertion.
   {
     APFloat f1(-8.85242279E-41f);
     APFloat f2(2.0f);
     APFloat f3(8.85242279E-41f);
     f1.fusedMultiplyAdd(f2, f3, APFloat::rmNearestTiesToEven);
     EXPECT_EQ(-8.85242279E-41f, f1.convertToFloat());
   }
+
+  // Test using only a single instance of APFloat.
+  {
+    APFloat F(1.5);
+
+    F.fusedMultiplyAdd(F, F, APFloat::rmNearestTiesToEven);
+    EXPECT_EQ(3.75, F.convertToDouble());
+  }
 }
 
 TEST(APFloatTest, MinNum) {
   APFloat f1(1.0);
   APFloat f2(2.0);
   APFloat nan = APFloat::getNaN(APFloat::IEEEdouble());
 
   EXPECT_EQ(1.0, minnum(f1, f2).convertToDouble());