[libc][math][c++23] Implement basic arithmetic operations for BFloat16 (#151228)

This PR implements addition, subtraction, multiplication and division
operations for BFloat16.

---------

Signed-off-by: krishna2803 <[email protected]>
Signed-off-by: Krishna Pandey <[email protected]>
Co-authored-by: OverMighty <[email protected]>

Author: krishna2803

libc/src/__support/FPUtil/CMakeLists.txt

@@ -285,6 +285,9 @@ add_header_library(
     libc.hdr.stdint_proxy
     libc.src.__support.CPP.bit
     libc.src.__support.CPP.type_traits
+    libc.src.__support.FPUtil.generic.add_sub
+    libc.src.__support.FPUtil.generic.div
+    libc.src.__support.FPUtil.generic.mul
     libc.src.__support.macros.config
     libc.src.__support.macros.properties.types
 )

libc/src/__support/FPUtil/bfloat16.h

@@ -15,6 +15,9 @@
 #include "src/__support/FPUtil/cast.h"
 #include "src/__support/FPUtil/comparison_operations.h"
 #include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/FPUtil/generic/add_sub.h"
+#include "src/__support/FPUtil/generic/div.h"
+#include "src/__support/FPUtil/generic/mul.h"
 #include "src/__support/macros/config.h"
 #include "src/__support/macros/properties/types.h"
 
@@ -81,6 +84,28 @@ struct BFloat16 {
   LIBC_INLINE bool operator>=(BFloat16 other) const {
     return fputil::greater_than_or_equals(*this, other);
   }
+
+  LIBC_INLINE constexpr BFloat16 operator-() const {
+    fputil::FPBits<bfloat16> result(*this);
+    result.set_sign(result.is_pos() ? Sign::NEG : Sign::POS);
+    return result.get_val();
+  }
+
+  LIBC_INLINE BFloat16 operator+(BFloat16 other) const {
+    return fputil::generic::add<BFloat16>(*this, other);
+  }
+
+  LIBC_INLINE BFloat16 operator-(BFloat16 other) const {
+    return fputil::generic::sub<BFloat16>(*this, other);
+  }
+
+  LIBC_INLINE BFloat16 operator*(BFloat16 other) const {
+    return fputil::generic::mul<bfloat16>(*this, other);
+  }
+
+  LIBC_INLINE BFloat16 operator/(BFloat16 other) const {
+    return fputil::generic::div<bfloat16>(*this, other);
+  }
 }; // struct BFloat16
 
 } // namespace fputil

libc/src/__support/FPUtil/cast.h

@@ -27,47 +27,47 @@ LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<OutType> &&
                                        OutType>
 cast(InType x) {
   // Casting to the same type is a no-op.
-  if constexpr (cpp::is_same_v<InType, OutType>)
+  if constexpr (cpp::is_same_v<InType, OutType>) {
     return x;
-
-  // bfloat16 is always defined (for now)
-  if constexpr (cpp::is_same_v<OutType, bfloat16> ||
-                cpp::is_same_v<InType, bfloat16>
+  } else {
+    if constexpr (cpp::is_same_v<OutType, bfloat16> ||
+                  cpp::is_same_v<InType, bfloat16>
 #if defined(LIBC_TYPES_HAS_FLOAT16) && !defined(__LIBC_USE_FLOAT16_CONVERSION)
-                || cpp::is_same_v<OutType, float16> ||
-                cpp::is_same_v<InType, float16>
+                  || cpp::is_same_v<OutType, float16> ||
+                  cpp::is_same_v<InType, float16>
 #endif
-  ) {
-    using InFPBits = FPBits<InType>;
-    using InStorageType = typename InFPBits::StorageType;
-    using OutFPBits = FPBits<OutType>;
-    using OutStorageType = typename OutFPBits::StorageType;
+    ) {
+      using InFPBits = FPBits<InType>;
+      using InStorageType = typename InFPBits::StorageType;
+      using OutFPBits = FPBits<OutType>;
+      using OutStorageType = typename OutFPBits::StorageType;
 
-    InFPBits x_bits(x);
+      InFPBits x_bits(x);
 
-    if (x_bits.is_nan()) {
-      if (x_bits.is_signaling_nan()) {
-        raise_except_if_required(FE_INVALID);
-        return OutFPBits::quiet_nan().get_val();
-      }
+      if (x_bits.is_nan()) {
+        if (x_bits.is_signaling_nan()) {
+          raise_except_if_required(FE_INVALID);
+          return OutFPBits::quiet_nan().get_val();
+        }
 
-      InStorageType x_mant = x_bits.get_mantissa();
-      if (InFPBits::FRACTION_LEN > OutFPBits::FRACTION_LEN)
-        x_mant >>= InFPBits::FRACTION_LEN - OutFPBits::FRACTION_LEN;
-      return OutFPBits::quiet_nan(x_bits.sign(),
-                                  static_cast<OutStorageType>(x_mant))
-          .get_val();
-    }
+        InStorageType x_mant = x_bits.get_mantissa();
+        if (InFPBits::FRACTION_LEN > OutFPBits::FRACTION_LEN)
+          x_mant >>= InFPBits::FRACTION_LEN - OutFPBits::FRACTION_LEN;
+        return OutFPBits::quiet_nan(x_bits.sign(),
+                                    static_cast<OutStorageType>(x_mant))
+            .get_val();
+      }
 
-    if (x_bits.is_inf())
-      return OutFPBits::inf(x_bits.sign()).get_val();
+      if (x_bits.is_inf())
+        return OutFPBits::inf(x_bits.sign()).get_val();
 
-    constexpr size_t MAX_FRACTION_LEN =
-        cpp::max(OutFPBits::FRACTION_LEN, InFPBits::FRACTION_LEN);
-    DyadicFloat<cpp::bit_ceil(MAX_FRACTION_LEN)> xd(x);
-    return xd.template as<OutType, /*ShouldSignalExceptions=*/true>();
-  } else {
-    return static_cast<OutType>(x);
+      constexpr size_t MAX_FRACTION_LEN =
+          cpp::max(OutFPBits::FRACTION_LEN, InFPBits::FRACTION_LEN);
+      DyadicFloat<cpp::bit_ceil(MAX_FRACTION_LEN)> xd(x);
+      return xd.template as<OutType, /*ShouldSignalExceptions=*/true>();
+    } else {
+      return static_cast<OutType>(x);
+    }
   }
 }
 

libc/src/__support/FPUtil/dyadic_float.h

@@ -576,7 +576,7 @@ LIBC_INLINE constexpr DyadicFloat<Bits> quick_mul(const DyadicFloat<Bits> &a,
     // Check the leading bit directly, should be faster than using clz in
     // normalize().
     if (result.mantissa.val[DyadicFloat<Bits>::MantissaType::WORD_COUNT - 1] >>
-            63 ==
+            (DyadicFloat<Bits>::MantissaType::WORD_SIZE - 1) ==
         0)
       result.shift_left(1);
   } else {

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

@@ -68,6 +68,7 @@ add_header_library(
     libc.src.__support.FPUtil.rounding_mode
     libc.src.__support.macros.attributes
     libc.src.__support.macros.optimization
+    libc.src.__support.macros.properties.types
 )
 
 add_header_library(
@@ -77,6 +78,7 @@ add_header_library(
   DEPENDS
     libc.hdr.errno_macros
     libc.hdr.fenv_macros
+    libc.src.__support.CPP.algorithm
     libc.src.__support.CPP.bit
     libc.src.__support.CPP.type_traits
     libc.src.__support.FPUtil.basic_operations

libc/src/__support/FPUtil/generic/add_sub.h

@@ -104,13 +104,22 @@ add_or_sub(InType x, InType y) {
         }
       }
 
-      // volatile prevents Clang from converting tmp to OutType and then
-      // immediately back to InType before negating it, resulting in double
-      // rounding.
-      volatile InType tmp = y;
-      if constexpr (IsSub)
-        tmp = -tmp;
-      return cast<OutType>(tmp);
+      if constexpr (cpp::is_same_v<InType, bfloat16> &&
+                    cpp::is_same_v<OutType, bfloat16>) {
+        OutFPBits y_bits(y);
+        if constexpr (IsSub)
+          y_bits.set_sign(y_bits.sign().negate());
+        return y_bits.get_val();
+      } else {
+
+        // volatile prevents Clang from converting tmp to OutType and then
+        // immediately back to InType before negating it, resulting in double
+        // rounding.
+        volatile InType tmp = y;
+        if constexpr (IsSub)
+          tmp = -tmp;
+        return cast<OutType>(tmp);
+      }
     }
 
     if (y_bits.is_zero())

libc/src/__support/FPUtil/generic/div.h

@@ -11,6 +11,7 @@
 
 #include "hdr/errno_macros.h"
 #include "hdr/fenv_macros.h"
+#include "src/__support/CPP/algorithm.h"
 #include "src/__support/CPP/bit.h"
 #include "src/__support/CPP/type_traits.h"
 #include "src/__support/FPUtil/BasicOperations.h"
@@ -34,8 +35,9 @@ div(InType x, InType y) {
   using OutStorageType = typename OutFPBits::StorageType;
   using InFPBits = FPBits<InType>;
   using InStorageType = typename InFPBits::StorageType;
-  using DyadicFloat =
-      DyadicFloat<cpp::bit_ceil(static_cast<size_t>(InFPBits::SIG_LEN + 1))>;
+  using DyadicFloat = DyadicFloat<cpp::max(
+      static_cast<size_t>(16),
+      cpp::bit_ceil(static_cast<size_t>(InFPBits::SIG_LEN + 1)))>;
 
   InFPBits x_bits(x);
   InFPBits y_bits(y);

libc/test/src/math/exhaustive/CMakeLists.txt

@@ -567,3 +567,75 @@ add_fp_unittest(
   LINK_LIBRARIES
     -lpthread
 )
+
+add_fp_unittest(
+  bfloat16_add_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    bfloat16_add_test.cpp
+  COMPILE_OPTIONS
+    ${libc_opt_high_flag}
+  DEPENDS
+    .exhaustive_test
+    libc.src.__support.FPUtil.bfloat16
+    libc.src.__support.FPUtil.fp_bits
+  LINK_LIBRARIES
+    -lpthread
+)
+
+add_fp_unittest(
+  bfloat16_div_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    bfloat16_div_test.cpp
+  COMPILE_OPTIONS
+    ${libc_opt_high_flag}
+  DEPENDS
+    .exhaustive_test
+    libc.src.__support.FPUtil.bfloat16
+    libc.src.__support.FPUtil.fp_bits
+  LINK_LIBRARIES
+    -lpthread
+)
+
+add_fp_unittest(
+  bfloat16_mul_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    bfloat16_mul_test.cpp
+  COMPILE_OPTIONS
+    ${libc_opt_high_flag}
+  DEPENDS
+    .exhaustive_test
+    libc.src.__support.FPUtil.bfloat16
+    libc.src.__support.FPUtil.fp_bits
+  LINK_LIBRARIES
+    -lpthread
+)
+
+add_fp_unittest(
+  bfloat16_sub_test
+  NO_RUN_POSTBUILD
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    bfloat16_sub_test.cpp
+  COMPILE_OPTIONS
+    ${libc_opt_high_flag}
+  DEPENDS
+    .exhaustive_test
+    libc.src.__support.FPUtil.bfloat16
+    libc.src.__support.FPUtil.fp_bits
+  LINK_LIBRARIES
+    -lpthread
+)

libc/test/src/math/exhaustive/bfloat16_add_test.cpp

@@ -0,0 +1,65 @@
+//===-- Exhaustive tests for bfloat16 addition ----------------------------===//
+//
+// 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 "exhaustive_test.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/bfloat16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "utils/MPFRWrapper/MPCommon.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+using LIBC_NAMESPACE::fputil::BFloat16;
+
+static BFloat16 add_func(BFloat16 x, BFloat16 y) { return x + y; }
+
+struct Bfloat16AddChecker : public virtual LIBC_NAMESPACE::testing::Test {
+  using FloatType = BFloat16;
+  using FPBits = LIBC_NAMESPACE::fputil::FPBits<bfloat16>;
+  using StorageType = typename FPBits::StorageType;
+
+  uint64_t check(uint16_t x_start, uint16_t x_stop, uint16_t y_start,
+                 uint16_t y_stop, mpfr::RoundingMode rounding) {
+    mpfr::ForceRoundingMode r(rounding);
+    if (!r.success)
+      return true;
+    uint16_t xbits = x_start;
+    uint64_t failed = 0;
+    do {
+      BFloat16 x = FPBits(xbits).get_val();
+      uint16_t ybits = xbits;
+      do {
+        BFloat16 y = FPBits(ybits).get_val();
+        mpfr::BinaryInput<BFloat16> input{x, y};
+        bool correct = TEST_MPFR_MATCH_ROUNDING_SILENTLY(
+            mpfr::Operation::Add, input, add_func(x, y), 0.5, rounding);
+        failed += (!correct);
+      } while (ybits++ < y_stop);
+    } while (xbits++ < x_stop);
+    return failed;
+  }
+};
+
+using LlvmLibcBfloat16ExhaustiveAddTest =
+    LlvmLibcExhaustiveMathTest<Bfloat16AddChecker, 1 << 2>;
+
+// range: [0, inf]
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x7f80U;
+
+// range: [-0, -inf]
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xff80U;
+
+TEST_F(LlvmLibcBfloat16ExhaustiveAddTest, PositiveRange) {
+  test_full_range_all_roundings(POS_START, POS_STOP, POS_START, POS_STOP);
+}
+
+TEST_F(LlvmLibcBfloat16ExhaustiveAddTest, NegativeRange) {
+  test_full_range_all_roundings(NEG_START, NEG_STOP, NEG_START, NEG_STOP);
+}