[libc] Add support for x86-64 targets that do not have FMA instructions.

[lntue]

Summary

Make FMA flag checks more accurate for x86-64 targets, and refactor
polyeval to use multiply and add instead when FMA instructions are not
available.

[lntue]

Fix missing endif directive and bazel BUILD file.

[michaelrj]

LGTM as is, I looked into if this patch allowed for reverting this one: https://reviews.llvm.org/D122689 and added my results as a comment.

libc/src/math/generic/CMakeLists.txt

517

On windows, removing the -mfma options in this cmake file with this patch doesn't break any of the tests except for this one. The expm1f tests fail with a ULP error, copied below.

[ RUN      ] LlvmLibcExpm1fTest.InFloatRange
C:/src/llvm-project/libc/test/src/math/expm1f_test.cpp:119: FAILURE
Failed to match __llvm_libc::expm1f(x) against __llvm_libc::testing::mpfr::get_mpfr_matcher<mpfr::Operation::Expm1>( x, __llvm_libc::expm1f(x), 0.5, mpfr::RoundingMode::TowardZero).
Match value not within tolerance value of MPFR result:
  Input decimal: -0.00000000000000000000000000000000000001175976677754
     Input bits: 0x80800d72, (S | E | M) = (1 | 0x01 | 0x000d72)

  Match decimal: -0.00000000000000000000000000000000000001175976677754
     Match bits: 0x80800d72, (S | E | M) = (1 | 0x01 | 0x000d72)

    MPFR result: -0.00000000000000000000000000000000000001175976537624
   MPFR rounded: 0x80800d71, (S | E | M) = (1 | 0x01 | 0x000d71)

      ULP error: 1.000000
[  FAILED  ] LlvmLibcExpm1fTest.InFloatRange
Ran 5 tests.  PASS: 4  FAIL: 1

[lntue]

libc/src/math/generic/CMakeLists.txt

517

This is expected. There will be followup changes to take care of different exceptional values from math functions when FMA instructions are not supported.

[sivachandra]

libc/src/__support/architectures.h

41

Nit: a more indicative name is probably LIBC_TARGET_HAS_FMA? Drop the LLVM_ prefix also.

[lntue]

Change the macro name to LIBC_TARGET_HAS_FMA.

libc/src/__support/FPUtil/CMakeLists.txt

 #TODO(michaelrj): separate the fputil library into individual targets.
 add_header_library(
   fputil
   HDRS
     FEnvImpl.h
     BasicOperations.h
     DivisionAndRemainderOperations.h
     FloatProperties.h
     FPBits.h
     BasicOperations.h
     ManipulationFunctions.h
     NearestIntegerOperations.h
     NormalFloat.h
     PlatformDefs.h
-    PolyEval.h
     UInt.h
     XFloat.h
   DEPENDS
     libc.include.math
     libc.include.errno
     libc.include.fenv
     libc.src.__support.common
     libc.src.__support.CPP.type_traits
     libc.src.__support.CPP.bit
     libc.src.errno.errno
 )
 
 add_header_library(
   sqrt
   HDRS
     sqrt.h
   DEPENDS
     .fputil
     libc.src.__support.FPUtil.generic.sqrt
 )
 
+add_header_library(
+  fma
+  HDRS
+    FMA.h
+  DEPENDS
+    .fputil
+    libc.src.__support.FPUtil.generic.fma
+)
+
+add_header_library(
+  multiply_add
+  HDRS
+    multiply_add.h
+  DEPENDS
+    .fma
+)
+
+add_header_library(
+  polyeval
+  HDRS
+    PolyEval.h
+  DEPENDS
+    .multiply_add
+)
+
 add_subdirectory(generic)

libc/src/__support/FPUtil/FMA.h

 //===-- Common header for FMA implementations -------------------*- 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_FMA_H
 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_FMA_H
 
 #include "src/__support/architectures.h"
 
+#if defined(LIBC_TARGET_HAS_FMA)
+
 #if defined(LLVM_LIBC_ARCH_X86_64)
 #include "x86_64/FMA.h"
 #elif defined(LLVM_LIBC_ARCH_AARCH64)
 #include "aarch64/FMA.h"
+#endif
+
 #else
+// FMA instructions are not available
 #include "generic/FMA.h"
 #include "src/__support/CPP/TypeTraits.h"
 
 namespace __llvm_libc {
 namespace fputil {
 
 // We have a generic implementation available only for single precision fma as
 // we restrict it to float values for now.

libc/src/__support/FPUtil/PolyEval.h

 //===-- Common header for PolyEval implementations --------------*- 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_POLYEVAL_H
 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_POLYEVAL_H
 
-#include "src/__support/CPP/TypeTraits.h"
-#include "src/__support/architectures.h"
+#include "multiply_add.h"
 
 // Evaluate polynomial using Horner's Scheme:
 // With polyeval(x, a_0, a_1, ..., a_n) = a_n * x^n + ... + a_1 * x + a_0, we
 // evaluated it as:  a_0 + x * (a_1 + x * ( ... (a_(n-1) + x * a_n) ... ) ) ).
-// We will use fma instructions if available.
+// We will use FMA instructions if available.
 // Example: to evaluate x^3 + 2*x^2 + 3*x + 4, call
 //   polyeval( x, 4.0, 3.0, 2.0, 1.0 )
 
-#if defined(LLVM_LIBC_ARCH_X86_64) || defined(LLVM_LIBC_ARCH_AARCH64)
-#include "FMA.h"
-
 namespace __llvm_libc {
 namespace fputil {
 
 template <typename T> static inline T polyeval(T x, T a0) { return a0; }
 
 template <typename T, typename... Ts>
 INLINE_FMA static inline T polyeval(T x, T a0, Ts... a) {
-  return fma(x, polyeval(x, a...), a0);
+  return multiply_add(x, polyeval(x, a...), a0);
 }
 
 } // namespace fputil
 } // namespace __llvm_libc
 
-#ifdef LLVM_LIBC_ARCH_X86_64
-
-// [DISABLED] There is a regression with using vectorized version for polyeval
-// compared to the naive Horner's scheme with fma.  Need further investigation
-// #include "x86_64/PolyEval.h"
-
-#endif // LLVM_LIBC_ARCH_X86_64
-
-#else
-
-namespace __llvm_libc {
-namespace fputil {
-
-template <typename T> static inline T polyeval(T x, T a0) { return a0; }
-
-template <typename T, typename... Ts>
-static inline T polyeval(T x, T a0, Ts... a) {
-  return x * polyeval(x, a...) + a0;
-}
-
-} // namespace fputil
-} // namespace __llvm_libc
-
-#endif
-
-#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_FMA_H
+#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_POLYEVAL_H

libc/src/__support/FPUtil/aarch64/FMA.h

 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_AARCH64_FMA_H
 
 #include "src/__support/architectures.h"
 
 #if !defined(LLVM_LIBC_ARCH_AARCH64)
 #error "Invalid include"
 #endif
 
+#if !defined(LIBC_TARGET_HAS_FMA)
+#error "FMA instructions are not supported"
+#endif
+
 #include "src/__support/CPP/TypeTraits.h"
 
 namespace __llvm_libc {
 namespace fputil {
 
 template <typename T>
 cpp::EnableIfType<cpp::IsSame<T, float>::Value, T> fma(T x, T y, T z) {
   float result;

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

 add_header_library(
   sqrt
   HDRS
     sqrt.h
     sqrt_80_bit_long_double.h
 )
+
+add_header_library(
+  fma
+  HDRS
+    FMA.h
+)

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

 //===-- Common header for FMA implementations -------------------*- 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_GENERIC_FMA_H
 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_GENERIC_FMA_H
 
 #include "src/__support/CPP/TypeTraits.h"
+#include "src/__support/FPUtil/FPBits.h"
 
 namespace __llvm_libc {
 namespace fputil {
 namespace generic {
 
 template <typename T>
 static inline cpp::EnableIfType<cpp::IsSame<T, float>::Value, T> fma(T x, T y,
                                                                      T z) {

libc/src/__support/FPUtil/multiply_add.h

+//===-- Common header for multiply-add implementations ----------*- 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_MULTIPLY_ADD_H
+#define LLVM_LIBC_SRC_SUPPORT_FPUTIL_MULTIPLY_ADD_H
+
+#include "src/__support/architectures.h"
+
+namespace __llvm_libc {
+namespace fputil {
+
+// Implement a simple wrapper for multiply-add operation:
+//   multiply_add(x, y, z) = x*y + z
+// which uses FMA instructions to speed up if available.
+
+template <typename T> static inline T multiply_add(T x, T y, T z) {
+  return x * y + z;
+}
+
+#if defined(LIBC_TARGET_HAS_FMA)
+// FMA instructions are available.
+#include "FMA.h"
+
+template <> inline float multiply_add<float>(float x, float y, float z) {
+  return fma(x, y, z);
+}
+
+template <> inline double multiply_add<double>(double x, double y, double z) {
+  return fma(x, y, z);
+}
+#endif // LIBC_TARGET_HAS_FMA
+
+} // namespace fputil
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_MULTIPLY_ADD_H

libc/src/__support/FPUtil/x86_64/FMA.h

 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_FMA_H
 
 #include "src/__support/architectures.h"
 
 #if !defined(LLVM_LIBC_ARCH_X86_64)
 #error "Invalid include"
 #endif
 
+#if !defined(LIBC_TARGET_HAS_FMA)
+#error "FMA instructions are not supported"
+#endif
+
 #include "src/__support/CPP/TypeTraits.h"
 #include <immintrin.h>
 
 namespace __llvm_libc {
 namespace fputil {
 
 template <typename T>
 INLINE_FMA static inline cpp::EnableIfType<cpp::IsSame<T, float>::Value, T>

libc/src/__support/architectures.h

 #if defined(__aarch64__) || defined(__arm64__) || defined(_M_ARM64)
 #define LLVM_LIBC_ARCH_AARCH64
 #endif
 
 #if (defined(LLVM_LIBC_ARCH_AARCH64) || defined(LLVM_LIBC_ARCH_ARM))
 #define LLVM_LIBC_ARCH_ANY_ARM
 #endif
 
-#if defined(LLVM_LIBC_ARCH_X86_64)
+#if defined(LLVM_LIBC_ARCH_AARCH64)
+#define LIBC_TARGET_HAS_FMA
+#elif defined(LLVM_LIBC_ARCH_X86_64)
+#if (defined(__AVX2__) || defined(__FMA__))
+#define LIBC_TARGET_HAS_FMA
+#endif
+#endif
+
+#if (defined(LLVM_LIBC_ARCH_X86_64) && defined(LIBC_TARGET_HAS_FMA))
 #define INLINE_FMA __attribute__((target("fma")))
 #else
 #define INLINE_FMA
 #endif // LLVM_LIBC_ARCH_X86_64
 
 #endif // LLVM_LIBC_SUPPORT_ARCHITECTURES_H

libc/src/math/CMakeLists.txt

 add_entrypoint_object(
   fmaf
   SRCS
     fmaf.cpp
   HDRS
     fmaf.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.fma
   COMPILE_OPTIONS
-    -O2
+    -O3
     -mfma
 )
 
 add_entrypoint_object(
   fma
   SRCS
     fma.cpp
   HDRS
     fma.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.fma
   COMPILE_OPTIONS
-    -O2
+    -O3
     -mfma
 )
 
 add_math_entrypoint_object(ceil)
 add_math_entrypoint_object(ceilf)
 add_math_entrypoint_object(ceill)
 
 add_math_entrypoint_object(copysign)

libc/src/math/generic/CMakeLists.txt

add_entrypoint_object(
   expf
   SRCS
     expf.cpp
   HDRS
     ../expf.h
   DEPENDS
     .common_constants
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
     -mfma
 )
 
 add_entrypoint_object(
   exp2f
   SRCS
     exp2f.cpp
   HDRS
     ../exp2f.h
   DEPENDS
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
     -mfma
 )
 
 add_entrypoint_object(
   expm1f
   SRCS
     expm1f.cpp
   HDRS
     ../expm1f.h
   DEPENDS
     .common_constants
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.polyeval
     libc.include.math
   COMPILE_OPTIONS
     -O3
     -mfma
 )
 
 add_entrypoint_object(
   copysign
   SRCS
     copysign.cpp
   HDRS
     ../copysign.h
add_entrypoint_object(
   log10f
   SRCS
     log10f.cpp
   HDRS
     ../log10f.h
   DEPENDS
     .common_constants
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.polyeval
   COMPILE_OPTIONS
     -O3
     -mfma
 )
 
 add_entrypoint_object(
   log1pf
   SRCS
     log1pf.cpp
   HDRS
     ../log1pf.h
   DEPENDS
     .common_constants
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.polyeval
   COMPILE_OPTIONS
     -O3
     -mfma
 )
 
 add_entrypoint_object(
   log2f
   SRCS
     log2f.cpp
   HDRS
     ../log2f.h
   DEPENDS
     .common_constants
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.polyeval
     COMPILE_OPTIONS
     -O3
     -mfma
 )
 
 add_entrypoint_object(
   logf
   SRCS
     logf.cpp
   HDRS
     ../logf.h
   DEPENDS
     .common_constants
     libc.src.__support.FPUtil.fputil
+    libc.src.__support.FPUtil.multiply_add
+    libc.src.__support.FPUtil.polyeval
   COMPILE_OPTIONS
     -O3
     -mfma
 )
 
 add_entrypoint_object(
   logb
   SRCS

libc/src/math/generic/expm1f.cpp

if (x_abs < 0x3300'0000U) {
       //                               < 2^-25
       //                               < epsilon(1)/2.
       // So the correctly rounded values of expm1(x) are:
       //   = x + eps(x) if rounding mode = FE_UPWARD,
       //                   or (rounding mode = FE_TOWARDZERO and x is negative),
       //   = x otherwise.
       // To simplify the rounding decision and make it more efficient, we use
       //   fma(x, x, x) ~ x + x^2 instead.
-      return fputil::fma(x, x, x);
+      return fputil::multiply_add(x, x, x);
     }
 
     // 2^-25 <= |x| < 2^-4
     double xd = static_cast<double>(x);
     double xsq = xd * xd;
     // Degree-8 minimax polynomial generated by Sollya with:
     // > display = hexadecimal;
     // > P = fpminimax((expm1(x) - x)/x^2, 6, [|D...|], [-2^-4, 2^-4]);
     double r =
         fputil::polyeval(xd, 0x1p-1, 0x1.55555555557ddp-3, 0x1.55555555552fap-5,
                          0x1.111110fcd58b7p-7, 0x1.6c16c1717660bp-10,
                          0x1.a0241f0006d62p-13, 0x1.a01e3f8d3c06p-16);
-    return static_cast<float>(fputil::fma(r, xsq, xd));
+    return static_cast<float>(fputil::multiply_add(r, xsq, xd));
   }
 
   // For -18 < x < 89, to compute expm1(x), we perform the following range
   // reduction: find hi, mid, lo such that:
   //   x = hi + mid + lo, in which
   //     hi is an integer,
   //     mid * 2^7 is an integer
   //     -2^(-8) <= lo < 2^-8.
LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
   // Degree-4 minimax polynomial generated by Sollya with the following
   // commands:
   //   > display = hexadecimal;
   //   > Q = fpminimax(expm1(x)/x, 3, [|D...|], [-2^-8, 2^-8]);
   //   > Q;
   double exp_lo =
       fputil::polyeval(xd, 0x1.0p0, 0x1.ffffffffff777p-1, 0x1.000000000071cp-1,
                        0x1.555566668e5e7p-3, 0x1.55555555ef243p-5);
-  return static_cast<float>(fputil::fma(exp_hi_mid, exp_lo, -1.0));
+  return static_cast<float>(fputil::multiply_add(exp_hi_mid, exp_lo, -1.0));
 }
 
 } // namespace __llvm_libc

libc/src/math/generic/log10f.cpp

LLVM_LIBC_FUNCTION(float, log10f, (float x)) {
   int f_index = xbits.get_mantissa() >> 16;
 
   FPBits f = xbits;
   f.bits &= ~0x0000'FFFF;
 
   double d = static_cast<float>(xbits) - static_cast<float>(f);
   d *= ONE_OVER_F[f_index];
 
-  double extra_factor = fputil::fma(m, LOG10_2, LOG10_F[f_index]);
+  double extra_factor = fputil::multiply_add(m, LOG10_2, LOG10_F[f_index]);
 
   double r = fputil::polyeval(d, extra_factor, 0x1.bcb7b1526e4c5p-2,
                               -0x1.bcb7b1518a5e9p-3, 0x1.287a72a6f716p-3,
                               -0x1.bcadb40b85565p-4, 0x1.5e0bc97f97e22p-4);
 
   return static_cast<float>(r);
 }
 
 } // namespace __llvm_libc

libc/src/math/generic/log1pf.cpp

INLINE_FMA static inline float log(double x) {
 
   FPBits f = xbits;
   // Clear the lowest 45 bits.
   f.bits &= ~0x0000'1FFF'FFFF'FFFFULL;
 
   double d = static_cast<double>(xbits) - static_cast<double>(f);
   d *= ONE_OVER_F[f_index];
 
-  double extra_factor = fputil::fma(m, LOG_2, LOG_F[f_index]);
+  double extra_factor = fputil::multiply_add(m, LOG_2, LOG_F[f_index]);
 
   double r = fputil::polyeval(d, extra_factor, 0x1.fffffffffffacp-1,
                               -0x1.fffffffef9cb2p-2, 0x1.5555513bc679ap-2,
                               -0x1.fff4805ea441p-3, 0x1.930180dbde91ap-3);
 
   return static_cast<float>(r);
 }
 
case 0xbb0ec8c4U: // x = -0x1.1d9188p-9f
     break;
   }
 
   double r;
   // Polymial generated with Sollya:
   // > fpminimax(log(1 + x)/x, 5, [|D...|], [-2^-8; 2^-8]);
   r = fputil::polyeval(xd, -0x1p-1, 0x1.5555555515551p-2, -0x1.ffffffff82bdap-3,
                        0x1.999b33348d3aep-3, -0x1.5556cae3adcc3p-3);
-  return static_cast<float>(fputil::fma(r, xd * xd, xd));
+  return static_cast<float>(fputil::multiply_add(r, xd * xd, xd));
 }
 
 } // namespace __llvm_libc

libc/src/math/generic/logf.cpp

LLVM_LIBC_FUNCTION(float, logf, (float x)) {
 
   FPBits f = xbits;
   f.bits &= ~0x0000'FFFF;
 
   double d = static_cast<float>(xbits) - static_cast<float>(f);
   d *= ONE_OVER_F[f_index];
 
   double extra_factor =
-      fputil::fma(static_cast<double>(m), LOG_2, LOG_F[f_index]);
+      fputil::multiply_add(static_cast<double>(m), LOG_2, LOG_F[f_index]);
 
   double r = __llvm_libc::fputil::polyeval(
       d, extra_factor, 0x1.fffffffffffacp-1, -0x1.fffffffef9cb2p-2,
       0x1.5555513bc679ap-2, -0x1.fff4805ea441p-3, 0x1.930180dbde91ap-3);
 
   return static_cast<float>(r);
 }
 
 } // namespace __llvm_libc

utils/bazel/llvm-project-overlay/libc/BUILD.bazel

deps = [
         ":__support_common",
         ":__support_cpp_bit",
         ":__support_cpp_type_traits",
         ":__support_fputil",
         ":libc_root",
     ],
 )
 
+fma_common_hdrs = [
+    "src/__support/FPUtil/FMA.h",
+    "src/__support/FPUtil/generic/FMA.h",
+]
+
+fma_hdrs = selects.with_or({
+    "//conditions:default": fma_common_hdrs,
+    PLATFORM_CPU_X86_64: fma_common_hdrs + [
+        "src/__support/FPUtil/x86_64/FMA.h",
+    ],
+    PLATFORM_CPU_ARM64: fma_common_hdrs + [
+        "src/__support/FPUtil/aarch64/FMA.h",
+    ],
+})
+
+cc_library(
+    name = "__support_fputil_fma",
+    hdrs = fma_hdrs,
+    deps = [
+        ":__support_common",
+        ":__support_cpp_bit",
+        ":__support_cpp_type_traits",
+        ":__support_fputil",
+        ":libc_root",
+    ],
+)
+
+cc_library(
+    name = "__support_fputil_multiply_add",
+    hdrs = [
+        "src/__support/FPUtil/multiply_add.h",
+    ],
+    deps = [
+        ":__support_common",
+        ":__support_fputil_fma",
+    ],
+)
+
+cc_library(
+    name = "__support_fputil_polyeval",
+    hdrs = [
+        "src/__support/FPUtil/PolyEval.h",
+    ],
+    deps = [
+        ":__support_fputil_multiply_add",
+    ],
+)
+
 ################################ fenv targets ################################
 
 libc_function(
     name = "fetestexcept",
     srcs = ["src/fenv/fetestexcept.cpp"],
     hdrs = ["src/fenv/fetestexcept.h"],
     deps = [
         ":__support_common",