[APFloat] Add PPCDoubleDouble multiplication

[timshen]

Rebase and add more comments.

Now that D27872 is checked in (and bots don't seem to freak out), this patch should be pretty straightforward.

[hfinkel]

LGTM

llvm/lib/Support/APFloat.cpp

if (Semantics == RHS.Semantics && RHS.Floats) {
     Floats[1] = RHS.Floats[1];
   } else if (this != &RHS) {
     this->~DoubleAPFloat();
     new (this) DoubleAPFloat(RHS);
   }
   return *this;
 }
 
+// Implement addition, subtraction, multiplication and division based on:
 // "Software for Doubled-Precision Floating-Point Computations",
 // by Seppo Linnainmaa, ACM TOMS vol 7 no 3, September 1981, pages 272-283.
 APFloat::opStatus DoubleAPFloat::addImpl(const APFloat &a, const APFloat &aa,
                                          const APFloat &c, const APFloat &cc,
                                          roundingMode RM) {
   int Status = opOK;
   APFloat z = a;
   Status |= z.add(c, RM);
APFloat::opStatus DoubleAPFloat::subtract(const DoubleAPFloat &RHS,
   changeSign();
   auto Ret = add(RHS, RM);
   changeSign();
   return Ret;
 }
 
 APFloat::opStatus DoubleAPFloat::multiply(const DoubleAPFloat &RHS,
                                           APFloat::roundingMode RM) {
-  assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
-  APFloat Tmp(semPPCDoubleDoubleLegacy, bitcastToAPInt());
-  auto Ret =
-      Tmp.multiply(APFloat(semPPCDoubleDoubleLegacy, RHS.bitcastToAPInt()), RM);
-  *this = DoubleAPFloat(semPPCDoubleDouble, Tmp.bitcastToAPInt());
-  return Ret;
+  const auto &LHS = *this;
+  auto &Out = *this;
+  /* Interesting observation: For special categories, finding the lowest
+     common ancestor of the following layered graph gives the correct
+     return category:
+
+        NaN
+       /   \
+     Zero  Inf
+       \   /
+       Normal
+
+     e.g. NaN * NaN = NaN
+          Zero * Inf = NaN
+          Normal * Zero = Zero
+          Normal * Inf = Inf
+  */
+  if (LHS.getCategory() == fcNaN) {
+    Out = LHS;
+    return opOK;
+  }
+  if (RHS.getCategory() == fcNaN) {
+    Out = RHS;
+    return opOK;
+  }
+  if ((LHS.getCategory() == fcZero && RHS.getCategory() == fcInfinity) ||
+      (LHS.getCategory() == fcInfinity && RHS.getCategory() == fcZero)) {
+    Out.makeNaN(false, false, nullptr);
+    return opOK;
+  }
+  if (LHS.getCategory() == fcZero || LHS.getCategory() == fcInfinity) {
+    Out = LHS;
+    return opOK;
+  }
+  if (RHS.getCategory() == fcZero || RHS.getCategory() == fcInfinity) {
+    Out = RHS;
+    return opOK;
+  }
+  assert(LHS.getCategory() == fcNormal && RHS.getCategory() == fcNormal &&
+         "Special cases not handled exhaustively");
+
+  int Status = opOK;
+  APFloat A = Floats[0], B = Floats[1], C = RHS.Floats[0], D = RHS.Floats[1];
+  // t = a * c
+  APFloat T = A;
+  Status |= T.multiply(C, RM);
+  if (!T.isFiniteNonZero()) {
+    Floats[0] = T;
+    Floats[1].makeZero(false);
+    return (opStatus)Status;
+  }
+
+  // tau = fmsub(a, c, t), that is -fmadd(-a, c, t).
+  APFloat Tau = A;
+  T.changeSign();
+  Status |= Tau.fusedMultiplyAdd(C, T, RM);
+  T.changeSign();
+  {
+    // v = a * d
+    APFloat V = A;
+    Status |= V.multiply(D, RM);
+    // w = b * c
+    APFloat W = B;
+    Status |= W.multiply(C, RM);
+    Status |= V.add(W, RM);
+    // tau += v + w
+    Status |= Tau.add(V, RM);
+  }
+  // u = t + tau
+  APFloat U = T;
+  Status |= U.add(Tau, RM);
+
+  Floats[0] = U;
+  if (!U.isFinite()) {
+    Floats[1].makeZero(false);
+  } else {
+    // Floats[1] = (t - u) + tau
+    Status |= T.subtract(U, RM);
+    Status |= T.add(Tau, RM);
+    Floats[1] = T;
+  }
+  return (opStatus)Status;
 }
 
 APFloat::opStatus DoubleAPFloat::divide(const DoubleAPFloat &RHS,
                                         APFloat::roundingMode RM) {
   assert(Semantics == &semPPCDoubleDouble && "Unexpected Semantics");
   APFloat Tmp(semPPCDoubleDoubleLegacy, bitcastToAPInt());
   auto Ret =
       Tmp.divide(APFloat(semPPCDoubleDoubleLegacy, RHS.bitcastToAPInt()), RM);

llvm/unittests/ADT/APFloatTest.cpp

TEST(APFloatTest, PPCDoubleDoubleAddSpecial) {
   };
 
   for (auto Tp : Data) {
     uint64_t Op1[2], Op2[2];
     APFloat::fltCategory Expected;
     APFloat::roundingMode RM;
     std::tie(Op1[0], Op1[1], Op2[0], Op2[1], Expected, RM) = Tp;
 
+    {
       APFloat A1(APFloat::PPCDoubleDouble(), APInt(128, 2, Op1));
       APFloat A2(APFloat::PPCDoubleDouble(), APInt(128, 2, Op2));
       A1.add(A2, RM);
 
       EXPECT_EQ(Expected, A1.getCategory())
-        << formatv("({0:x} + {1:x}) + ({2:x} + {3:x})", Op1[0], Op1[1], Op2[0],
-                   Op2[1])
+          << formatv("({0:x} + {1:x}) + ({2:x} + {3:x})", Op1[0], Op1[1],
+                     Op2[0], Op2[1])
+                 .str();
+    }
+    {
+      APFloat A1(APFloat::PPCDoubleDouble(), APInt(128, 2, Op1));
+      APFloat A2(APFloat::PPCDoubleDouble(), APInt(128, 2, Op2));
+      A2.add(A1, RM);
+
+      EXPECT_EQ(Expected, A2.getCategory())
+          << formatv("({0:x} + {1:x}) + ({2:x} + {3:x})", Op2[0], Op2[1],
+                     Op1[0], Op1[1])
                  .str();
     }
   }
+}
 
 TEST(APFloatTest, PPCDoubleDoubleAdd) {
   using DataType = std::tuple<uint64_t, uint64_t, uint64_t, uint64_t, uint64_t,
                               uint64_t, APFloat::roundingMode>;
   DataType Data[] = {
       // (1 + 0) + (1e-105 + 0) = (1 + 1e-105)
       std::make_tuple(0x3ff0000000000000ull, 0, 0x3960000000000000ull, 0,
                       0x3ff0000000000000ull, 0x3960000000000000ull,
DataType Data[] = {
                       0x7c8ffffffffffffeull, APFloat::rmNearestTiesToEven),
   };
 
   for (auto Tp : Data) {
     uint64_t Op1[2], Op2[2], Expected[2];
     APFloat::roundingMode RM;
     std::tie(Op1[0], Op1[1], Op2[0], Op2[1], Expected[0], Expected[1], RM) = Tp;
 
+    {
       APFloat A1(APFloat::PPCDoubleDouble(), APInt(128, 2, Op1));
       APFloat A2(APFloat::PPCDoubleDouble(), APInt(128, 2, Op2));
       A1.add(A2, RM);
 
       EXPECT_EQ(Expected[0], A1.bitcastToAPInt().getRawData()[0])
-        << formatv("({0:x} + {1:x}) + ({2:x} + {3:x})", Op1[0], Op1[1], Op2[0],
-                   Op2[1])
+          << formatv("({0:x} + {1:x}) + ({2:x} + {3:x})", Op1[0], Op1[1],
+                     Op2[0], Op2[1])
                  .str();
       EXPECT_EQ(Expected[1], A1.bitcastToAPInt().getRawData()[1])
-        << formatv("({0:x} + {1:x}) + ({2:x} + {3:x})", Op1[0], Op1[1], Op2[0],
-                   Op2[1])
+          << formatv("({0:x} + {1:x}) + ({2:x} + {3:x})", Op1[0], Op1[1],
+                     Op2[0], Op2[1])
                  .str();
     }
+    {
+      APFloat A1(APFloat::PPCDoubleDouble(), APInt(128, 2, Op1));
+      APFloat A2(APFloat::PPCDoubleDouble(), APInt(128, 2, Op2));
+      A2.add(A1, RM);
+
+      EXPECT_EQ(Expected[0], A2.bitcastToAPInt().getRawData()[0])
+          << formatv("({0:x} + {1:x}) + ({2:x} + {3:x})", Op2[0], Op2[1],
+                     Op1[0], Op1[1])
+                 .str();
+      EXPECT_EQ(Expected[1], A2.bitcastToAPInt().getRawData()[1])
+          << formatv("({0:x} + {1:x}) + ({2:x} + {3:x})", Op2[0], Op2[1],
+                     Op1[0], Op1[1])
+                 .str();
+    }
+  }
 }
 
 TEST(APFloatTest, PPCDoubleDoubleSubtract) {
   using DataType = std::tuple<uint64_t, uint64_t, uint64_t, uint64_t, uint64_t,
                               uint64_t, APFloat::roundingMode>;
   DataType Data[] = {
       // (1 + 0) - (-1e-105 + 0) = (1 + 1e-105)
       std::make_tuple(0x3ff0000000000000ull, 0, 0xb960000000000000ull, 0,
EXPECT_EQ(Expected[0], A1.bitcastToAPInt().getRawData()[0])
                .str();
     EXPECT_EQ(Expected[1], A1.bitcastToAPInt().getRawData()[1])
         << formatv("({0:x} + {1:x}) - ({2:x} + {3:x})", Op1[0], Op1[1], Op2[0],
                    Op2[1])
                .str();
   }
 }
 
+TEST(APFloatTest, PPCDoubleDoubleMultiplySpecial) {
+  using DataType = std::tuple<uint64_t, uint64_t, uint64_t, uint64_t,
+                              APFloat::fltCategory, APFloat::roundingMode>;
+  DataType Data[] = {
+      // fcNaN * fcNaN = fcNaN
+      std::make_tuple(0x7ff8000000000000ull, 0, 0x7ff8000000000000ull, 0,
+                      APFloat::fcNaN, APFloat::rmNearestTiesToEven),
+      // fcNaN * fcZero = fcNaN
+      std::make_tuple(0x7ff8000000000000ull, 0, 0, 0, APFloat::fcNaN,
+                      APFloat::rmNearestTiesToEven),
+      // fcNaN * fcInfinity = fcNaN
+      std::make_tuple(0x7ff8000000000000ull, 0, 0x7ff0000000000000ull, 0,
+                      APFloat::fcNaN, APFloat::rmNearestTiesToEven),
+      // fcNaN * fcNormal = fcNaN
+      std::make_tuple(0x7ff8000000000000ull, 0, 0x3ff0000000000000ull, 0,
+                      APFloat::fcNaN, APFloat::rmNearestTiesToEven),
+      // fcInfinity * fcInfinity = fcInfinity
+      std::make_tuple(0x7ff0000000000000ull, 0, 0x7ff0000000000000ull, 0,
+                      APFloat::fcInfinity, APFloat::rmNearestTiesToEven),
+      // fcInfinity * fcZero = fcNaN
+      std::make_tuple(0x7ff0000000000000ull, 0, 0, 0, APFloat::fcNaN,
+                      APFloat::rmNearestTiesToEven),
+      // fcInfinity * fcNormal = fcInfinity
+      std::make_tuple(0x7ff0000000000000ull, 0, 0x3ff0000000000000ull, 0,
+                      APFloat::fcInfinity, APFloat::rmNearestTiesToEven),
+      // fcZero * fcZero = fcZero
+      std::make_tuple(0, 0, 0, 0, APFloat::fcZero,
+                      APFloat::rmNearestTiesToEven),
+      // fcZero * fcNormal = fcZero
+      std::make_tuple(0, 0, 0x3ff0000000000000ull, 0, APFloat::fcZero,
+                      APFloat::rmNearestTiesToEven),
+  };
+
+  for (auto Tp : Data) {
+    uint64_t Op1[2], Op2[2];
+    APFloat::fltCategory Expected;
+    APFloat::roundingMode RM;
+    std::tie(Op1[0], Op1[1], Op2[0], Op2[1], Expected, RM) = Tp;
+
+    {
+      APFloat A1(APFloat::PPCDoubleDouble(), APInt(128, 2, Op1));
+      APFloat A2(APFloat::PPCDoubleDouble(), APInt(128, 2, Op2));
+      A1.multiply(A2, RM);
+
+      EXPECT_EQ(Expected, A1.getCategory())
+          << formatv("({0:x} + {1:x}) * ({2:x} + {3:x})", Op1[0], Op1[1],
+                     Op2[0], Op2[1])
+                 .str();
+    }
+    {
+      APFloat A1(APFloat::PPCDoubleDouble(), APInt(128, 2, Op1));
+      APFloat A2(APFloat::PPCDoubleDouble(), APInt(128, 2, Op2));
+      A2.multiply(A1, RM);
+
+      EXPECT_EQ(Expected, A2.getCategory())
+          << formatv("({0:x} + {1:x}) * ({2:x} + {3:x})", Op2[0], Op2[1],
+                     Op1[0], Op1[1])
+                 .str();
+    }
+  }
+}
+
 TEST(APFloatTest, PPCDoubleDoubleMultiply) {
   using DataType = std::tuple<uint64_t, uint64_t, uint64_t, uint64_t, uint64_t,
                               uint64_t, APFloat::roundingMode>;
-  // TODO: Only a sanity check for now. Add more edge cases when the
-  // double-double algorithm is implemented.
   DataType Data[] = {
       // 1/3 * 3 = 1.0
       std::make_tuple(0x3fd5555555555555ull, 0x3c75555555555556ull,
                       0x4008000000000000ull, 0, 0x3ff0000000000000ull, 0,
                       APFloat::rmNearestTiesToEven),
+      // (1 + epsilon) * (1 + 0) = fcZero
+      std::make_tuple(0x3ff0000000000000ull, 0x0000000000000001ull,
+                      0x3ff0000000000000ull, 0, 0x3ff0000000000000ull,
+                      0x0000000000000001ull, APFloat::rmNearestTiesToEven),
+      // (1 + epsilon) * (1 + epsilon) = 1 + 2 * epsilon
+      std::make_tuple(0x3ff0000000000000ull, 0x0000000000000001ull,
+                      0x3ff0000000000000ull, 0x0000000000000001ull,
+                      0x3ff0000000000000ull, 0x0000000000000002ull,
+                      APFloat::rmNearestTiesToEven),
+      // -(1 + epsilon) * (1 + epsilon) = -1
+      std::make_tuple(0xbff0000000000000ull, 0x0000000000000001ull,
+                      0x3ff0000000000000ull, 0x0000000000000001ull,
+                      0xbff0000000000000ull, 0, APFloat::rmNearestTiesToEven),
+      // (0.5 + 0) * (1 + 2 * epsilon) = 0.5 + epsilon
+      std::make_tuple(0x3fe0000000000000ull, 0, 0x3ff0000000000000ull,
+                      0x0000000000000002ull, 0x3fe0000000000000ull,
+                      0x0000000000000001ull, APFloat::rmNearestTiesToEven),
+      // (0.5 + 0) * (1 + epsilon) = 0.5
+      std::make_tuple(0x3fe0000000000000ull, 0, 0x3ff0000000000000ull,
+                      0x0000000000000001ull, 0x3fe0000000000000ull, 0,
+                      APFloat::rmNearestTiesToEven),
+      // __LDBL_MAX__ * (1 + 1 << 106) = inf
+      std::make_tuple(0x7fefffffffffffffull, 0x7c8ffffffffffffeull,
+                      0x3ff0000000000000ull, 0x3950000000000000ull,
+                      0x7ff0000000000000ull, 0, APFloat::rmNearestTiesToEven),
+      // __LDBL_MAX__ * (1 + 1 << 107) > __LDBL_MAX__, but not inf, yes =_=|||
+      std::make_tuple(0x7fefffffffffffffull, 0x7c8ffffffffffffeull,
+                      0x3ff0000000000000ull, 0x3940000000000000ull,
+                      0x7fefffffffffffffull, 0x7c8fffffffffffffull,
+                      APFloat::rmNearestTiesToEven),
+      // __LDBL_MAX__ * (1 + 1 << 108) = __LDBL_MAX__
+      std::make_tuple(0x7fefffffffffffffull, 0x7c8ffffffffffffeull,
+                      0x3ff0000000000000ull, 0x3930000000000000ull,
+                      0x7fefffffffffffffull, 0x7c8ffffffffffffeull,
+                      APFloat::rmNearestTiesToEven),
   };
 
   for (auto Tp : Data) {
     uint64_t Op1[2], Op2[2], Expected[2];
     APFloat::roundingMode RM;
     std::tie(Op1[0], Op1[1], Op2[0], Op2[1], Expected[0], Expected[1], RM) = Tp;
 
+    {
       APFloat A1(APFloat::PPCDoubleDouble(), APInt(128, 2, Op1));
       APFloat A2(APFloat::PPCDoubleDouble(), APInt(128, 2, Op2));
       A1.multiply(A2, RM);
 
       EXPECT_EQ(Expected[0], A1.bitcastToAPInt().getRawData()[0])
-        << formatv("({0:x} + {1:x}) * ({2:x} + {3:x})", Op1[0], Op1[1], Op2[0],
-                   Op2[1])
+          << formatv("({0:x} + {1:x}) * ({2:x} + {3:x})", Op1[0], Op1[1],
+                     Op2[0], Op2[1])
                  .str();
       EXPECT_EQ(Expected[1], A1.bitcastToAPInt().getRawData()[1])
-        << formatv("({0:x} + {1:x}) * ({2:x} + {3:x})", Op1[0], Op1[1], Op2[0],
-                   Op2[1])
+          << formatv("({0:x} + {1:x}) * ({2:x} + {3:x})", Op1[0], Op1[1],
+                     Op2[0], Op2[1])
                  .str();
     }
+    {
+      APFloat A1(APFloat::PPCDoubleDouble(), APInt(128, 2, Op1));
+      APFloat A2(APFloat::PPCDoubleDouble(), APInt(128, 2, Op2));
+      A2.multiply(A1, RM);
+
+      EXPECT_EQ(Expected[0], A2.bitcastToAPInt().getRawData()[0])
+          << formatv("({0:x} + {1:x}) * ({2:x} + {3:x})", Op2[0], Op2[1],
+                     Op1[0], Op1[1])
+                 .str();
+      EXPECT_EQ(Expected[1], A2.bitcastToAPInt().getRawData()[1])
+          << formatv("({0:x} + {1:x}) * ({2:x} + {3:x})", Op2[0], Op2[1],
+                     Op1[0], Op1[1])
+                 .str();
+    }
+  }
 }
 
 TEST(APFloatTest, PPCDoubleDoubleDivide) {
   using DataType = std::tuple<uint64_t, uint64_t, uint64_t, uint64_t, uint64_t,
                               uint64_t, APFloat::roundingMode>;
   // TODO: Only a sanity check for now. Add more edge cases when the
   // double-double algorithm is implemented.
   DataType Data[] = {