[libc][math][c23] fmul correcly rounded to all rounding modes (#91537)

Job Henandez Lara · web-flow · commit 263be9fb0085 · 2024-06-08T15:07:27.000-04:00
This is an implementation of floating point multiplication:

It will consist of 
   - `double x double -&gt; float`
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
@@ -394,6 +394,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fminimum_mag_num
     libc.src.math.fminimum_mag_numf
     libc.src.math.fminimum_mag_numl
+    libc.src.math.fmul
     libc.src.math.fmod
     libc.src.math.fmodf
     libc.src.math.fmodl
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
@@ -261,6 +261,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fminimum_mag_num
     libc.src.math.fminimum_mag_numf
     libc.src.math.fminimum_mag_numl
+    libc.src.math.fmul
     libc.src.math.fmod
     libc.src.math.fmodf
     libc.src.math.frexp
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
@@ -402,6 +402,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fminimum_mag_num
     libc.src.math.fminimum_mag_numf
     libc.src.math.fminimum_mag_numl
+    libc.src.math.fmul
     libc.src.math.fmod
     libc.src.math.fmodf
     libc.src.math.fmodl
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
@@ -421,6 +421,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fminimum_mag_num
     libc.src.math.fminimum_mag_numf
     libc.src.math.fminimum_mag_numl
+    libc.src.math.fmul
     libc.src.math.fmod
     libc.src.math.fmodf
     libc.src.math.fmodl
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
@@ -180,6 +180,7 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.fminimum_mag_num
     libc.src.math.fminimum_mag_numf
     libc.src.math.fminimum_mag_numl
+    libc.src.math.fmul
     libc.src.math.fmod
     libc.src.math.fmodf
     libc.src.math.fmodl
diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst
@@ -158,7 +158,7 @@ Basic Operations
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | fmod             | |check|          | |check|         | |check|                | |check|              | |check|                | 7.12.10.1              | F.10.7.1                   |
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| fmul             | N/A              |                 |                        | N/A                  |                        | 7.12.14.3              | F.10.11                    |
+| fmul             | N/A              | |check|         |                        | N/A                  |                        | 7.12.14.3              | F.10.11                    |
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
 | frexp            | |check|          | |check|         | |check|                |                      | |check|                | 7.12.6.7               | F.10.3.7                   |
 +------------------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
@@ -472,6 +472,9 @@ def StdC : StandardSpec<"stdc"> {
           GuardedFunctionSpec<"fminimum_mag_numf16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>, ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
           GuardedFunctionSpec<"fminimum_mag_numf128", RetValSpec<Float128Type>, [ArgSpec<Float128Type>, ArgSpec<Float128Type>], "LIBC_TYPES_HAS_FLOAT128">,
 
+          FunctionSpec<"fmul", RetValSpec<FloatType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
+
+
           FunctionSpec<"fma", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"fmaf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>, ArgSpec<FloatType>]>,
 
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
@@ -180,6 +180,8 @@ add_math_entrypoint_object(fminimum_mag_numl)
 add_math_entrypoint_object(fminimum_mag_numf16)
 add_math_entrypoint_object(fminimum_mag_numf128)
 
+add_math_entrypoint_object(fmul)
+
 add_math_entrypoint_object(fmod)
 add_math_entrypoint_object(fmodf)
 add_math_entrypoint_object(fmodl)
diff --git a/libc/src/math/fmul.h b/libc/src/math/fmul.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for fmul --------------------------*- 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_MATH_FMUL_H
+#define LLVM_LIBC_SRC_MATH_FMUL_H
+
+namespace LIBC_NAMESPACE {
+
+float fmul(double x, double y);
+
+} // namespace LIBC_NAMESPACE
+
+#endif // LLVM_LIBC_SRC_MATH_FMUL_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
@@ -2354,6 +2354,22 @@ add_entrypoint_object(
     -O3
 )
 
+add_entrypoint_object(
+  fmul
+  SRCS
+    fmul.cpp
+  HDRS
+    ../fmul.h
+  DEPENDS
+    libc.src.__support.FPUtil.basic_operations
+    libc.src.__support.uint128
+    libc.src.__support.CPP.bit
+    libc.src.__support.FPUtil.fp_bits
+    libc.src.__support.FPUtil.rounding_mode
+  COMPILE_OPTIONS
+    -O3
+)
+
 add_entrypoint_object(
   sqrt
   SRCS
diff --git a/libc/src/math/generic/fmul.cpp b/libc/src/math/generic/fmul.cpp
@@ -0,0 +1,128 @@
+//===-- Implementation of fmul function------------------------------------===//
+//
+// 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/math/fmul.h"
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/BasicOperations.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/common.h"
+#include "src/__support/uint128.h"
+
+namespace LIBC_NAMESPACE {
+
+LLVM_LIBC_FUNCTION(float, fmul, (double x, double y)) {
+  auto x_bits = fputil::FPBits<double>(x);
+
+  auto y_bits = fputil::FPBits<double>(y);
+
+  auto output_sign = (x_bits.sign() != y_bits.sign()) ? Sign::NEG : Sign::POS;
+
+  if (LIBC_UNLIKELY(x_bits.is_inf_or_nan() || y_bits.is_inf_or_nan() ||
+                    x_bits.is_zero() || y_bits.is_zero())) {
+    if (x_bits.is_nan())
+      return static_cast<float>(x);
+    if (y_bits.is_nan())
+      return static_cast<float>(y);
+    if (x_bits.is_inf())
+      return y_bits.is_zero()
+                 ? fputil::FPBits<float>::quiet_nan().get_val()
+                 : fputil::FPBits<float>::inf(output_sign).get_val();
+    if (y_bits.is_inf())
+      return x_bits.is_zero()
+                 ? fputil::FPBits<float>::quiet_nan().get_val()
+                 : fputil::FPBits<float>::inf(output_sign).get_val();
+    // Now either x or y is zero, and the other one is finite.
+    return fputil::FPBits<float>::zero(output_sign).get_val();
+  }
+
+  uint64_t mx, my;
+
+  // Get mantissa and append the hidden bit if needed.
+  mx = x_bits.get_explicit_mantissa();
+  my = y_bits.get_explicit_mantissa();
+
+  // Get the corresponding biased exponent.
+  int ex = x_bits.get_explicit_exponent();
+  int ey = y_bits.get_explicit_exponent();
+
+  // Count the number of leading zeros of the explicit mantissas.
+  int nx = cpp::countl_zero(mx);
+  int ny = cpp::countl_zero(my);
+  // Shift the leading 1 bit to the most significant bit.
+  mx <<= nx;
+  my <<= ny;
+
+  // Adjust exponent accordingly: If x or y are normal, we will only need to
+  // shift by (exponent length + sign bit = 11 bits. If x or y are denormal, we
+  // will need to shift more than 11 bits.
+  ex -= (nx - 11);
+  ey -= (ny - 11);
+
+  UInt128 product = static_cast<UInt128>(mx) * static_cast<UInt128>(my);
+  int32_t dm1;
+  uint64_t highs, lows;
+  uint64_t g, hight, lowt;
+  uint32_t m;
+  uint32_t b;
+  int c;
+
+  highs = static_cast<uint64_t>(product >> 64);
+  c = static_cast<int>(highs >= 0x8000000000000000);
+  lows = static_cast<uint64_t>(product);
+
+  lowt = (lows != 0);
+
+  dm1 = ex + ey + c + fputil::FPBits<float>::EXP_BIAS;
+
+  int round_mode = fputil::quick_get_round();
+  if (dm1 >= 255) {
+    if ((round_mode == FE_TOWARDZERO) ||
+        (round_mode == FE_UPWARD && output_sign.is_neg()) ||
+        (round_mode == FE_DOWNWARD && output_sign.is_pos())) {
+      return fputil::FPBits<float>::max_normal(output_sign).get_val();
+    }
+    return fputil::FPBits<float>::inf().get_val();
+  } else if (dm1 <= 0) {
+
+    int m_shift = 40 + c - dm1;
+    int g_shift = m_shift - 1;
+    int h_shift = 64 - g_shift;
+    m = (m_shift >= 64) ? 0 : static_cast<uint32_t>(highs >> m_shift);
+
+    g = g_shift >= 64 ? 0 : (highs >> g_shift) & 1;
+    hight = h_shift >= 64 ? highs : (highs << h_shift) != 0;
+
+    dm1 = 0;
+  } else {
+    m = static_cast<uint32_t>(highs >> (39 + c));
+    g = (highs >> (38 + c)) & 1;
+    hight = (highs << (26 - c)) != 0;
+  }
+
+  if (round_mode == FE_TONEAREST) {
+    b = g && ((hight && lowt) || ((m & 1) != 0));
+  } else if ((output_sign.is_neg() && round_mode == FE_DOWNWARD) ||
+             (output_sign.is_pos() && round_mode == FE_UPWARD)) {
+    b = (g == 0 && (hight && lowt) == 0) ? 0 : 1;
+  } else {
+    b = 0;
+  }
+
+  uint32_t exp16 = (dm1 << 23);
+
+  uint32_t m2 = m & fputil::FPBits<float>::FRACTION_MASK;
+
+  uint32_t result = (exp16 + m2) + b;
+
+  auto result_bits = fputil::FPBits<float>(result);
+  result_bits.set_sign(output_sign);
+  return result_bits.get_val();
+}
+
+} // namespace LIBC_NAMESPACE
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
@@ -2400,6 +2400,19 @@ add_fp_unittest(
     libc.src.__support.FPUtil.fp_bits
 )
 
+add_fp_unittest(
+  fmul_test
+  SUITE
+    libc-math-smoke-tests
+  SRCS
+    fmul_test.cpp
+  HDRS
+    FMulTest.h
+  DEPENDS
+    libc.src.math.fmul
+    libc.src.__support.FPUtil.fp_bits
+)
+
 add_fp_unittest(
   sqrtf_test
   SUITE
diff --git a/libc/test/src/math/smoke/FMulTest.h b/libc/test/src/math/smoke/FMulTest.h
@@ -0,0 +1,104 @@
+//===-- Utility class to test fmul[f|l] ---------------------*- 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_TEST_SRC_MATH_SMOKE_FMULTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_SMOKE_FMULTEST_H
+
+#include "test/UnitTest/FEnvSafeTest.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+template <typename T, typename R>
+class FmulTest : public LIBC_NAMESPACE::testing::FEnvSafeTest {
+
+  DECLARE_SPECIAL_CONSTANTS(T)
+
+public:
+  typedef T (*FMulFunc)(R, R);
+
+  void testMul(FMulFunc func) {
+
+    EXPECT_FP_EQ_ALL_ROUNDING(T(15.0), func(3.0, 5.0));
+    EXPECT_FP_EQ_ALL_ROUNDING(T(0x1.0p-130), func(0x1.0p1, 0x1.0p-131));
+    EXPECT_FP_EQ_ALL_ROUNDING(T(0x1.0p-127), func(0x1.0p2, 0x1.0p-129));
+    EXPECT_FP_EQ_ALL_ROUNDING(T(1.0), func(1.0, 1.0));
+
+    EXPECT_FP_EQ_ALL_ROUNDING(T(0.0), func(-0.0, -0.0));
+    EXPECT_FP_EQ_ALL_ROUNDING(T(-0.0), func(0.0, -0.0));
+    EXPECT_FP_EQ_ALL_ROUNDING(T(-0.0), func(-0.0, 0.0));
+
+    EXPECT_FP_EQ_ROUNDING_NEAREST(inf, func(0x1.0p100, 0x1.0p100));
+    EXPECT_FP_EQ_ROUNDING_UPWARD(inf, func(0x1.0p100, 0x1.0p100));
+    EXPECT_FP_EQ_ROUNDING_DOWNWARD(max_normal, func(0x1.0p100, 0x1.0p100));
+    EXPECT_FP_EQ_ROUNDING_TOWARD_ZERO(max_normal, func(0x1.0p100, 0x1.0p100));
+
+    EXPECT_FP_EQ_ROUNDING_NEAREST(
+        0x1p0, func(1.0, 1.0 + 0x1.0p-128 + 0x1.0p-149 + 0x1.0p-150));
+    EXPECT_FP_EQ_ROUNDING_DOWNWARD(
+        0x1p0, func(1.0, 1.0 + 0x1.0p-128 + 0x1.0p-149 + 0x1.0p-150));
+    EXPECT_FP_EQ_ROUNDING_TOWARD_ZERO(
+        0x1p0, func(1.0, 1.0 + 0x1.0p-128 + 0x1.0p-149 + 0x1.0p-150));
+    EXPECT_FP_EQ_ROUNDING_UPWARD(
+        0x1p0, func(1.0, 1.0 + 0x1.0p-128 + 0x1.0p-149 + 0x1.0p-150));
+
+    EXPECT_FP_EQ_ROUNDING_NEAREST(
+        0x1.0p-128f + 0x1.0p-148f,
+        func(1.0, 0x1.0p-128 + 0x1.0p-149 + 0x1.0p-150));
+    EXPECT_FP_EQ_ROUNDING_UPWARD(
+        0x1.0p-128f + 0x1.0p-148f,
+        func(1.0, 0x1.0p-128 + 0x1.0p-149 + 0x1.0p-150));
+    EXPECT_FP_EQ_ROUNDING_DOWNWARD(
+        0x1.0p-128f + 0x1.0p-149f,
+        func(1.0, 0x1.0p-128 + 0x1.0p-149 + 0x1.0p-150));
+    EXPECT_FP_EQ_ROUNDING_TOWARD_ZERO(
+        0x1.0p-128f + 0x1.0p-149f,
+        func(1.0, 0x1.0p-128 + 0x1.0p-149 + 0x1.0p-150));
+  }
+
+  void testSpecialInputs(FMulFunc func) {
+    EXPECT_FP_EQ_ALL_ROUNDING(inf, func(inf, 0x1.0p-129));
+    EXPECT_FP_EQ_ALL_ROUNDING(inf, func(0x1.0p-129, inf));
+    EXPECT_FP_EQ_ALL_ROUNDING(inf, func(inf, 2.0));
+    EXPECT_FP_EQ_ALL_ROUNDING(inf, func(3.0, inf));
+    EXPECT_FP_EQ_ALL_ROUNDING(0.0, func(0.0, 0.0));
+
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(neg_inf, aNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(aNaN, neg_inf));
+    EXPECT_FP_EQ_ALL_ROUNDING(inf, func(neg_inf, neg_inf));
+
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(0.0, neg_inf));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(neg_inf, 0.0));
+
+    EXPECT_FP_EQ_ALL_ROUNDING(neg_inf, func(neg_inf, 1.0));
+    EXPECT_FP_EQ_ALL_ROUNDING(neg_inf, func(1.0, neg_inf));
+
+    EXPECT_FP_EQ_ALL_ROUNDING(neg_inf, func(neg_inf, 0x1.0p-129));
+    EXPECT_FP_EQ_ALL_ROUNDING(neg_inf, func(0x1.0p-129, neg_inf));
+
+    EXPECT_FP_EQ_ALL_ROUNDING(0.0, func(0.0, 0x1.0p-129));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(inf, 0.0));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(0.0, inf));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(0.0, aNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(2.0, aNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(0x1.0p-129, aNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(inf, aNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(aNaN, aNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(0.0, sNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(2.0, sNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(0x1.0p-129, sNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(inf, sNaN));
+    EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(sNaN, sNaN));
+  }
+};
+
+#define LIST_FMUL_TESTS(T, R, func)                                            \
+  using LlvmLibcFmulTest = FmulTest<T, R>;                                     \
+  TEST_F(LlvmLibcFmulTest, Mul) { testMul(&func); }                            \
+  TEST_F(LlvmLibcFmulTest, NaNInf) { testSpecialInputs(&func); }
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_SMOKE_FMULTEST_H
diff --git a/libc/test/src/math/smoke/fmul_test.cpp b/libc/test/src/math/smoke/fmul_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for fmul-------------------------------------------------===//
+//
+// 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 "FMulTest.h"
+
+#include "src/math/fmul.h"
+
+LIST_FMUL_TESTS(float, double, LIBC_NAMESPACE::fmul)
