[libc] Implement a high-precision floating point class.

[lntue]

Summary

Implement a high-precision floating point class using UInt<> as its
mantissa. This will be used in accurate pass for double precision math
functions.

[michaelrj]

Overall this looks promising, but I'd like to see some unit tests for the big float class if possible.

libc/src/__support/FPUtil/big_float.h

46 ↗

(On Diff #470951)

this could probably be templated by using FloatProperties<T>::MantissaWidth

75 ↗

(On Diff #470951)

same as above, I think you can make this either an operator or a templated to_float<T>

87 ↗

(On Diff #470951)

did you mean to have this twice?

209 ↗

(On Diff #470951)

nit: please fix

libc/src/__support/UInt.h

28–31

The pair class should probably be in a separate file.

592–602

is it possible to have this as <T> instead of specifying 128, 196, and 256? Same below.

[sivachandra]

There is a similar class by name XFloat but it does not have addition operation: https://github.com/llvm/llvm-project/blob/main/libc/src/__support/FPUtil/XFloat.h
Remove it if you don't like it.

libc/src/__support/FPUtil/big_float.h

30 ↗

(On Diff #470951)

Explain why this is better or preferred over the conventional exponent.

42 ↗

(On Diff #470951)

You can drop this assumption by using NormalFloat.

180 ↗

(On Diff #470951)

You should explain why these functiona have quick_ prefix in their name, and why they are not member operators.

libc/src/__support/UInt.h

28–31

Agree - you can choose to add a cpp:pair class.

61

These functions just seem to me like specializations of the UInt ful_add and ful_mul. If you really need these separately, you need to add some commentary to justify.

590

There could be a symbolic constant for this number.

[tschuett]

libc/src/__support/FPUtil/big_float.h

22 ↗

(On Diff #470951)

What is libc`s stand on nested namespaces?

[orex]

libc/src/__support/FPUtil/big_float.h

32 ↗

(On Diff #470951)

static_assert in UInt is stricter than this one (https://github.com/llvm/llvm-project/blob/7f93ae808634e33e4dc9bce753c909aa5f9a6eb4/libc/src/__support/UInt.h#L25) Probably you should remove or adjust this?

55 ↗

(On Diff #470951)

Will it be better to put M_LENGTH >= MANTISSA_LENGTH? For "highly templated" cases it can be easy to generalize such conversion.

75 ↗

(On Diff #470951)

General idea. Probably it is good to have conversion to double double? It can be precise for 64 bits mantissa?

84 ↗

(On Diff #470951)

The if is not needed. The operations in the if statement probably can be even faster than check and conditional jump?

90 ↗

(On Diff #470951)

Nit: Probably it is good to have 53 54 and 55 like some constexpr?

120 ↗

(On Diff #470951)

I'm not an expert in C++ standard, but pair with one template argument looks suspicious for me. Does C++ standard saying how it should be expanded?

123 ↗

(On Diff #470951)

Is this static cast by standard drop higher part of the variable?

129 ↗

(On Diff #470951)

Am I understand right, that rounding and sign of zero in this function is not important?

136 ↗

(On Diff #470951)

Probably to put this normalization to dedicated function?

1. It is quite general.
2. Remove the code duplication.

[orex]

libc/src/__support/FPUtil/big_float.h

120 ↗

(On Diff #470951)

The comment is not valid. Sorry.

[lntue]

libc/src/__support/FPUtil/big_float.h

22 ↗

(On Diff #470951)

I think we are totally fine with nested namespaces. This is mostly my part of copy-pasting from other headers in the same directory.

32 ↗

(On Diff #470951)

I was debating about letting this class using other uint types for the mantissa. But it's probably better to just let the UInt class taking care of it. Removed this static_assert.

42 ↗

(On Diff #470951)

I was thinking about using NormalFloat but there are some difference in the expectation between them that will need some extra work to be done in order to merge this and NormalFloat:

• currently NormalFloat expected to use floating point type as the template parameter, while we want to explicitly control the bits

84 ↗

(On Diff #470951)

The current <<= op in UInt is non-trivial and more expensive than a conditional jump. Other option is to put the check for trivial case into UInt itself.

123 ↗

(On Diff #470951)

Yes, it is specified for unsigned integers.

129 ↗

(On Diff #470951)

Yes, for now its main purpose is to aid the computation of polynomial approximations. So we don't pay too much attention on signed zeroes and rounding, which in turns giving us a tiny bit of performance gain.

libc/src/__support/UInt.h

28–31

I was debating about it, because the standard std::pair has members named first and second, while I want to named them hi and lo to signify the high part and the low part of the numbers. I'm not sure what's the best way to have such distinction yet. Maybe rename this to cpp::number_pair or something similar?

[lntue]

Change the name of the class to DyadicFloat, add normalization, and add tests.

[lntue]

libc/src/__support/UInt.h

28–31

These are refactored and cleaned up in https://reviews.llvm.org/D137871 and https://reviews.llvm.org/D138182

61

The optimizations and updates for UInt add and mul are done in https://reviews.llvm.org/D137871 and https://reviews.llvm.org/D138541

[michaelrj]

LGTM with minor nits

libc/src/__support/FPUtil/dyadic_float.h

14

I think dyadic_float.h should include float_properties.h directly.

libc/src/__support/UInt.h

368–369

shouldn't this be above the ifdef __SIZEOF_INT128__?

[sivachandra]

Mostly LGTM. I am not yet accepting as I have asked for definitions.

libc/src/__support/FPUtil/dyadic_float.h

2

Will this class be ever used from outside of the math directory? If no, then it probably should live in src/math.

22

This sentence feels incomplete.

25

One can read the code but couple of things should be described here I think:

1. That bits is the number of bits of mantissa and that it should be a multiple of 64.
2. Define the meaning of normalization.
3. Add a note next to the constructors that the object created is normalized as per the above definition.

27

Type name should be MantissaType.

35

This should like have additional constraints. For example, should this allow a quad-precision number to be converted to DyadicFloat<64>?

74

Why is this assumption required?

109

Add definition of "Quick Add".

110

Again, why?

160

Same comments as those for quick_add.

[lntue]

Address comments.

libc/src/__support/FPUtil/dyadic_float.h

2

I plan to use this to replace NormalFloat and XFloat classes, along with their usages in DivisionAndRemainderOperations or NearestIntegerOperations. Once the replacement completes, we will move them to src/math together.

35

Added checks for mantissa lengths.

74

Added explanation.

[lntue]

libc/src/__support/UInt.h

368–369

I intentially put it only for non 128-bit / without builtin. Since the built-in 128-bit shifts are quite efficient, so adding an extra check for the special case s == 0 does not improve the performance much, but actually reduce the throughput due to extra branching. For bigger size, the extra branching cost is much smaller compared with the work saving from avoiding all the shifting loops.

[sivachandra]

libc/src/__support/FPUtil/dyadic_float.h

47

This can lead to truncation if Bits is less than MANTISSA_WIDTH. Should we add a static check here, may be extending the above enable_if:

template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T> && FloatProperties<T>::MANTISSA_WIDTH <= Bits, int> = 0>

121

Add a mathematical expression which illustrates what actually is quick_add doing. Something like:

quick_add.exponent = max(...)
// aligning exponents - explain why
quick_add.mantissa = a.mantissa + b.mantissa;

174

Same mathematical explanation as that for quick_add.

libc/src/__support/FPUtil/CMakeLists.txt

add_header_library(
   HDRS
     nearest_integer.h
   DEPENDS
     libc.src.__support.common
   FLAGS
     ROUND_OPT
 )
 
+add_header_library(
+  dyadic_float
+  HDRS
+    dyadic_float.h
+  DEPENDS
+    .float_properties
+    .fp_bits
+    .multiply_add
+    libc.src.__support.common
+    libc.src.__support.uint
+)
+
 add_subdirectory(generic)

libc/src/__support/FPUtil/dyadic_float.h

+//===-- A class to store high precision floating point numbers --*- C++ -*-===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_SUPPORT_FPUTIL_DYADIC_FLOAT_H
+#define LLVM_LIBC_SRC_SUPPORT_FPUTIL_DYADIC_FLOAT_H
+
+#include "FPBits.h"
+#include "FloatProperties.h"
+#include "multiply_add.h"
+#include "src/__support/CPP/type_traits.h"
+#include "src/__support/UInt.h"
+
+#include <stddef.h>
+
+namespace __llvm_libc::fputil {
+
+// A generic class to perform comuptations of high precision floating points.
+// We store the value in dyadic format, including 3 fields:
+//   sign    : boolean value - false means positive, true means negative
+//   exponent: the exponent value of the least significant bit of the mantissa.
+//   mantissa: unsigned integer of length `Bits`.
+// So the real value that is stored is:
+//   real value = (-1)^sign * 2^exponent * (mantissa as unsigned integer)
+// The stored data is normal if for non-zero mantissa, the leading bit is 1.
+// The outputs of the constructors and most functions will be normalized.
+// To simplify and improve the efficiency, many functions will assume that the
+// inputs are normal.
+template <size_t Bits> struct DyadicFloat {
+  using MantissaType = __llvm_libc::cpp::UInt<Bits>;
+
+  bool sign = false;
+  int exponent = 0;
+  MantissaType mantissa = MantissaType(0);
+
+  DyadicFloat() = default;
+
+  template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
+  DyadicFloat(T x) {
+    FPBits<T> x_bits(x);
+    sign = x_bits.get_sign();
+    exponent = x_bits.get_exponent() - FloatProperties<T>::MANTISSA_WIDTH;
+    mantissa = MantissaType(x_bits.get_explicit_mantissa());
+    normalize();
+  }
+
+  DyadicFloat(bool s, int e, MantissaType m)
+      : sign(s), exponent(e), mantissa(m) {
+    normalize();
+  };
+
+  // Normalizing the mantissa, bringing the leading 1 bit to the most
+  // significant bit.
+  DyadicFloat &normalize() {
+    if (!mantissa.is_zero()) {
+      int shift_length = static_cast<int>(mantissa.clz());
+      exponent -= shift_length;
+      mantissa.shift_left(static_cast<size_t>(shift_length));
+    }
+    return *this;
+  }
+
+  // Used for aligning exponents.  Output might not be normalized.
+  DyadicFloat &shift_left(int shift_length) {
+    exponent -= shift_length;
+    mantissa <<= static_cast<size_t>(shift_length);
+    return *this;
+  }
+
+  // Used for aligning exponents.  Output might not be normalized.
+  DyadicFloat &shift_right(int shift_length) {
+    exponent += shift_length;
+    mantissa >>= static_cast<size_t>(shift_length);
+    return *this;
+  }
+
+  // Assume that it is already normalized and output is also normal.
+  // Output is rounded correctly with respect to the current rounding mode.
+  // TODO(lntue): Test or add support for denormal output.
+  // TODO(lntue): Test or add specialization for x86 long double.
+  template <typename T, typename = cpp::enable_if_t<
+                            cpp::is_floating_point_v<T> &&
+                                (FloatProperties<T>::MANTISSA_WIDTH <= Bits),
+                            void>>
+  explicit operator T() const {
+    // TODO(lntue): Do we need to treat signed zeros properly?
+    if (mantissa.is_zero())
+      return 0.0;
+
+    // Assume that it is normalized, and output is also normal.
+    constexpr size_t PRECISION = FloatProperties<T>::MANTISSA_WIDTH + 1;
+    using output_bits_t = typename FPBits<T>::UIntType;
+
+    MantissaType m_hi(mantissa >> (Bits - PRECISION));
+    auto d_hi = FPBits<T>::create_value(
+        sign, exponent + (Bits - 1) + FloatProperties<T>::EXPONENT_BIAS,
+        output_bits_t(m_hi) & FloatProperties<T>::MANTISSA_MASK);
+
+    const MantissaType ROUND_MASK = MantissaType(1) << (Bits - PRECISION - 1);
+    const MantissaType STICKY_MASK = ROUND_MASK - MantissaType(1);
+
+    bool round_bit = !(mantissa & ROUND_MASK).is_zero();
+    bool sticky_bit = !(mantissa & STICKY_MASK).is_zero();
+    int round_and_sticky = int(round_bit) * 2 + int(sticky_bit);
+    auto d_lo = FPBits<T>::create_value(sign,
+                                        exponent + (Bits - PRECISION - 2) +
+                                            FloatProperties<T>::EXPONENT_BIAS,
+                                        output_bits_t(0));
+
+    // Still correct without FMA instructions if `d_lo` is not underflow.
+    return multiply_add(d_lo.get_val(), T(round_and_sticky), d_hi.get_val());
+  }
+};
+
+// Quick add - Add 2 dyadic floats with rounding toward 0 and then normalize the
+// output.  The absolute errors compared to the mathematical sum is bounded by:
+//   | quick_add(a, b) - (a + b) | < MSB(a + b) * 2^(-Bits + 2),
+// i.e., errors are up to 2 ULPs.
+// Assume inputs are normalized (by constructors or other functions) so that we
+// don't need to normalize the inputs again in this function.  If the inputs are
+// not normalized, the results might lose precision significantly.
+template <size_t Bits>
+constexpr DyadicFloat<Bits> quick_add(DyadicFloat<Bits> a,
+                                      DyadicFloat<Bits> b) {
+  if (unlikely(a.mantissa.is_zero()))
+    return b;
+  if (unlikely(b.mantissa.is_zero()))
+    return a;
+
+  // Align exponents
+  if (a.exponent > b.exponent)
+    b.shift_right(a.exponent - b.exponent);
+  else if (b.exponent > a.exponent)
+    a.shift_right(b.exponent - a.exponent);
+
+  DyadicFloat<Bits> result;
+
+  if (a.sign == b.sign) {
+    // Addition
+    result.sign = a.sign;
+    result.exponent = a.exponent;
+    result.mantissa = a.mantissa;
+    if (result.mantissa.add(b.mantissa)) {
+      // Mantissa addition overflow.
+      result.shift_right(1);
+      result.mantissa.val[DyadicFloat<Bits>::MantissaType::WordCount - 1] |=
+          (uint64_t(1) << 63);
+    }
+    // Result is already normalized.
+    return result;
+  }
+
+  // Subtraction
+  if (a.mantissa >= b.mantissa) {
+    result.sign = a.sign;
+    result.exponent = a.exponent;
+    result.mantissa = a.mantissa - b.mantissa;
+  } else {
+    result.sign = b.sign;
+    result.exponent = b.exponent;
+    result.mantissa = b.mantissa - a.mantissa;
+  }
+
+  return result.normalize();
+}
+
+// Quick Mul - Slightly less accurate but efficient multiplication of 2 dyadic
+// floats with rounding toward 0 and then normalize the output.  The errors
+// compared to the mathematical product is bounded by:
+//   2 * errors of quick_mul_hi = 2 * (UInt<Bits>::WordCount - 1) in ULPs.
+// Assume inputs are normalized (by constructors or other functions) so that we
+// don't need to normalize the inputs again in this function.  If the inputs are
+// not normalized, the results might lose precision significantly.
+template <size_t Bits>
+constexpr DyadicFloat<Bits> quick_mul(DyadicFloat<Bits> a,
+                                      DyadicFloat<Bits> b) {
+  DyadicFloat<Bits> result;
+  result.sign = (a.sign != b.sign);
+  result.exponent = a.exponent + b.exponent + int(Bits);
+
+  if (!(a.mantissa.is_zero() || b.mantissa.is_zero())) {
+    result.mantissa = a.mantissa.quick_mul_hi(b.mantissa);
+    // Check the leading bit directly, should be faster than using clz in
+    // normalize().
+    if (result.mantissa.val[DyadicFloat<Bits>::MantissaType::WordCount - 1] >>
+            63 ==
+        0)
+      result.shift_left(1);
+  } else {
+    result.mantissa = (typename DyadicFloat<Bits>::MantissaType)(0);
+  }
+  return result;
+}
+
+} // namespace __llvm_libc::fputil
+
+#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_DYADIC_FLOAT_H

libc/src/__support/UInt.h

 //===-- A class to manipulate wide integers. --------------------*- C++ -*-===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef LLVM_LIBC_SRC_SUPPORT_UINT_H
 #define LLVM_LIBC_SRC_SUPPORT_UINT_H
 
 #include "src/__support/CPP/array.h"
 #include "src/__support/CPP/limits.h"
 #include "src/__support/CPP/optional.h"
 #include "src/__support/CPP/type_traits.h"
 #include "src/__support/builtin_wrappers.h"
+#include "src/__support/common.h"
 #include "src/__support/integer_utils.h"
 #include "src/__support/number_pair.h"
 
 #include <stddef.h> // For size_t
 #include <stdint.h>
 
 namespace __llvm_libc::cpp {
 
 template <size_t Bits> struct UInt {
 
   static_assert(Bits > 0 && Bits % 64 == 0,
                 "Number of bits in UInt should be a multiple of 64.");
   static constexpr size_t WordCount = Bits / 64;
   uint64_t val[WordCount];
 
   static constexpr uint64_t MASK32 = 0xFFFFFFFFu;
 
   static constexpr uint64_t low(uint64_t v) { return v & MASK32; }
   static constexpr uint64_t high(uint64_t v) { return (v >> 32) & MASK32; }
 
   constexpr UInt() {}
 
if (OtherBits >= Bits) {
       for (; i < WordCount; ++i)
         val[i] = 0;
     }
   }
 
   // Construct a UInt from a C array.
   template <size_t N, enable_if_t<N <= WordCount, int> = 0>
   constexpr UInt(const uint64_t (&nums)[N]) {
     size_t min_wordcount = N < WordCount ? N : WordCount;
     size_t i = 0;
     for (; i < min_wordcount; ++i)
       val[i] = nums[i];
 
     // If nums doesn't completely fill val, then fill the rest with zeroes.
     for (; i < WordCount; ++i)
       val[i] = 0;
   }
template <size_t Bits> struct UInt {
   }
 
   UInt<Bits> &operator=(const UInt<Bits> &other) {
     for (size_t i = 0; i < WordCount; ++i)
       val[i] = other.val[i];
     return *this;
   }
 
+  constexpr bool is_zero() const {
+    for (size_t i = 0; i < WordCount; ++i) {
+      if (val[i] != 0)
+        return false;
+    }
+    return true;
+  }
+
   // Add x to this number and store the result in this number.
   // Returns the carry value produced by the addition operation.
   constexpr uint64_t add(const UInt<Bits> &x) {
     SumCarry<uint64_t> s{0, 0};
     for (size_t i = 0; i < WordCount; ++i) {
       s = add_with_carry(val[i], x.val[i], s.carry);
       val[i] = s.sum;
     }
if constexpr (Bits == 128) {
       }
       __uint128_t tmp = __uint128_t(val[0]) + (__uint128_t(val[1]) << 64);
       tmp <<= s;
       val[0] = uint64_t(tmp);
       val[1] = uint64_t(tmp >> 64);
       return;
     }
 #endif                           // __SIZEOF_INT128__
+    if (unlikely(s == 0))
+      return;
+
     const size_t drop = s / 64;  // Number of words to drop
     const size_t shift = s % 64; // Bits to shift in the remaining words.
     size_t i = WordCount;
 
     if (drop < WordCount) {
       i = WordCount - 1;
       size_t j = WordCount - 1 - drop;
       for (; j > 0; --i, --j) {
if constexpr (Bits == 128) {
       __uint128_t tmp = __uint128_t(val[0]) + (__uint128_t(val[1]) << 64);
       tmp >>= s;
       val[0] = uint64_t(tmp);
       val[1] = uint64_t(tmp >> 64);
       return;
     }
 #endif // __SIZEOF_INT128__
 
+    if (unlikely(s == 0))
+      return;
     const size_t drop = s / 64;  // Number of words to drop
     const size_t shift = s % 64; // Bit shift in the remaining words.
 
     size_t i = 0;
 
     if (drop < WordCount) {
       if (shift > 0) {
         for (size_t j = drop; j < WordCount - 1; ++i, ++j) {
#endif // __SIZEOF_INT128__
   const uint64_t *data() const { return val; }
 };
 
 // Provides limits of UInt<128>.
 template <> class numeric_limits<UInt<128>> {
 public:
   static constexpr UInt<128> max() { return ~UInt<128>(0); }
   static constexpr UInt<128> min() { return 0; }
 };
 
 // Provides is_integral of UInt<128>, UInt<192>, UInt<256>.
 template <size_t Bits> struct is_integral<UInt<Bits>> : public cpp::true_type {
   static_assert(Bits > 0 && Bits % 64 == 0,
                 "Number of bits in UInt should be a multiple of 64.");
 };
 
 // Provides is_unsigned of UInt<128>, UInt<192>, UInt<256>.
 template <size_t Bits> struct is_unsigned<UInt<Bits>> : public cpp::true_type {
   static_assert(Bits > 0 && Bits % 64 == 0,
                 "Number of bits in UInt should be a multiple of 64.");
 };
 
 } // namespace __llvm_libc::cpp
 
 #endif // LLVM_LIBC_SRC_SUPPORT_UINT_H

libc/test/src/__support/CMakeLists.txt

add_custom_command(TARGET libc_str_to_float_comparison_test
                    COMMAND $<TARGET_FILE:libc_str_to_float_comparison_test> ${float_test_file}
                    DEPENDS ${float_test_file}
                    COMMENT "Test the strtof and strtod implementations against precomputed results." 
                    VERBATIM)
 
 add_subdirectory(CPP)
 add_subdirectory(File)
 add_subdirectory(OSUtil)
+add_subdirectory(FPUtil)

libc/test/src/__support/FPUtil/CMakeLists.txt

+add_libc_testsuite(libc_fputil_unittests)
+
+add_fp_unittest(
+  dyadic_float_test
+  NEED_MPFR
+  SUITE
+    libc_fputil_unittests
+  SRCS
+    dyadic_float_test.cpp
+  DEPENDS
+    libc.src.__support.FPUtil.dyadic_float
+)

libc/test/src/__support/FPUtil/dyadic_float_test.cpp

+//===-- Unittests for the DyadicFloat class -------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/UInt.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+
+using Float128 = __llvm_libc::fputil::DyadicFloat<128>;
+using Float192 = __llvm_libc::fputil::DyadicFloat<192>;
+using Float256 = __llvm_libc::fputil::DyadicFloat<256>;
+
+TEST(LlvmLibcDyadicFloatTest, BasicConversions) {
+  Float128 x(/*sign*/ false, /*exponent*/ 0,
+             /*mantissa*/ Float128::MantissaType(1));
+  volatile float xf = float(x);
+  volatile double xd = double(x);
+  ASSERT_FP_EQ(1.0f, xf);
+  ASSERT_FP_EQ(1.0, xd);
+
+  Float128 y(0x1.0p-53);
+  volatile float yf = float(y);
+  volatile double yd = double(y);
+  ASSERT_FP_EQ(0x1.0p-53f, yf);
+  ASSERT_FP_EQ(0x1.0p-53, yd);
+
+  Float128 z = quick_add(x, y);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(xf + yf, float(z));
+  EXPECT_FP_EQ_ALL_ROUNDING(xd + yd, double(z));
+}
+
+TEST(LlvmLibcDyadicFloatTest, QuickAdd) {
+  Float192 x(/*sign*/ false, /*exponent*/ 0,
+             /*mantissa*/ Float192::MantissaType(0x123456));
+  volatile double xd = double(x);
+  ASSERT_FP_EQ(0x1.23456p20, xd);
+
+  Float192 y(0x1.abcdefp-20);
+  volatile double yd = double(y);
+  ASSERT_FP_EQ(0x1.abcdefp-20, yd);
+
+  Float192 z = quick_add(x, y);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(xd + yd, (volatile double)(z));
+}
+
+TEST(LlvmLibcDyadicFloatTest, QuickMul) {
+  Float256 x(/*sign*/ false, /*exponent*/ 0,
+             /*mantissa*/ Float256::MantissaType(0x123456));
+  volatile double xd = double(x);
+  ASSERT_FP_EQ(0x1.23456p20, xd);
+
+  Float256 y(0x1.abcdefp-25);
+  volatile double yd = double(y);
+  ASSERT_FP_EQ(0x1.abcdefp-25, yd);
+
+  Float256 z = quick_mul(x, y);
+
+  EXPECT_FP_EQ_ALL_ROUNDING(xd * yd, double(z));
+}

libc/utils/UnitTest/CMakeLists.txt

 add_library(
   LibcFPTestHelpers
     FPExceptMatcher.cpp
     FPExceptMatcher.h
     FPMatcher.cpp
     FPMatcher.h
 )
 target_include_directories(LibcFPTestHelpers PUBLIC ${LIBC_SOURCE_DIR})
-target_link_libraries(LibcFPTestHelpers LibcUnitTest)
+target_link_libraries(LibcFPTestHelpers LibcUnitTest libc_test_utils)
 add_dependencies(
   LibcFPTestHelpers
   LibcUnitTest
   libc.utils.UnitTest.string_utils
   libc.src.__support.FPUtil.fp_bits
   libc.src.__support.FPUtil.fenv_impl
 )
 

libc/utils/UnitTest/FPMatcher.h

-//===-- TestMatchers.h ------------------------------------------*- C++ -*-===//
+//===-- FPMatchers.h --------------------------------------------*- C++ -*-===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #ifndef LLVM_LIBC_UTILS_UNITTEST_FPMATCHER_H
 #define LLVM_LIBC_UTILS_UNITTEST_FPMATCHER_H
 
 #include "src/__support/FPUtil/FEnvImpl.h"
 #include "src/__support/FPUtil/FPBits.h"
 #include "utils/UnitTest/Test.h"
+#include "utils/testutils/RoundingModeUtils.h"
 
 #include <errno.h>
 #include <math.h>
 
 namespace __llvm_libc {
 namespace fputil {
 namespace testing {
 

 
 #define ASSERT_FP_EXCEPTION(expected)                                          \
   do {                                                                         \
     if (math_errhandling & MATH_ERREXCEPT) {                                   \
       ASSERT_EQ(__llvm_libc::fputil::test_except(FE_ALL_EXCEPT), expected);    \
     }                                                                          \
   } while (0)
 
+#define EXPECT_FP_EQ_ALL_ROUNDING(expected, actual)                            \
+  do {                                                                         \
+    using namespace __llvm_libc::testutils;                                    \
+    ForceRoundingMode __r1(RoundingMode::Nearest);                             \
+    EXPECT_FP_EQ((expected), (actual));                                        \
+    ForceRoundingMode __r2(RoundingMode::Upward);                              \
+    EXPECT_FP_EQ((expected), (actual));                                        \
+    ForceRoundingMode __r3(RoundingMode::Downward);                            \
+    EXPECT_FP_EQ((expected), (actual));                                        \
+    ForceRoundingMode __r4(RoundingMode::TowardZero);                          \
+    EXPECT_FP_EQ((expected), (actual));                                        \
+  } while (0)
+
 #endif // LLVM_LIBC_UTILS_UNITTEST_FPMATCHER_H