Index: llvm/trunk/lib/Support/APFloat.cpp =================================================================== --- llvm/trunk/lib/Support/APFloat.cpp +++ llvm/trunk/lib/Support/APFloat.cpp @@ -3895,6 +3895,7 @@ 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, @@ -4037,12 +4038,88 @@ 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, Index: llvm/trunk/unittests/ADT/APFloatTest.cpp =================================================================== --- llvm/trunk/unittests/ADT/APFloatTest.cpp +++ llvm/trunk/unittests/ADT/APFloatTest.cpp @@ -3203,14 +3203,26 @@ 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]) - .str(); + { + 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]) + .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(); + } } } @@ -3255,18 +3267,34 @@ 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]) - .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(); + { + 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]) + .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(); + } + { + 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(); + } } } @@ -3304,16 +3332,111 @@ } } +TEST(APFloatTest, PPCDoubleDoubleMultiplySpecial) { + using DataType = std::tuple; + 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; - // 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) { @@ -3321,18 +3444,34 @@ 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]) - .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(); + { + 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]) + .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(); + } + { + 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(); + } } }