[libc] Use intrinsics for x86-64 fma and optimize PolyEval for x86-64 with degree 3 & 5 polynomials.

[lntue]

Summary

• Use intrinsics for x86-64 fma
• Optimize PolyEval for x86-64 with degree 3 & 5 polynomials.
• There might be a slight loss of accuracy compared to Horner's scheme due to usages of higher powers x^2 and x^3 in the computations.

[sivachandra]

The subject line and description are confusing. It seems to me like you have added specializations for 3rd and 5th degree polynomials. But, the description has "3-6 polynomials".

libc/src/__support/FPUtil/PolyEval.h

36

Couple of things:

1. We want individual header files to be self-contained. So, this include should be moved to the header files where it is required.
2. Because of the same self-contained headers requirement, the entities defined in a header file should be in the __llvm_libc namespace explicitly as opposed to ending up in a namespace when included like this here. So, these includes should be moved outside of the namespace here.

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

1 ↗

(On Diff #392768)

Fix line length but not sure if the description is correct also.

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

1 ↗

(On Diff #392768)

Fix line length but not sure if the description is correct also.

libc/src/math/generic/CMakeLists.txt

503 ↗

(On Diff #392768)

Can you add a comment explaining what in the implementation would be affected by this? AFAICT, we are either using the fma instructions directly, or are calling the fma related builtins. So, it is not clear to me as to why this should be required.

libc/test/src/math/expm1f_test.cpp

111

If there is a possibility of trading off between performance and accuracy, we should provide build time switches for users to pick one or the other based on their needs. By default we want to "err" on the side of providing more accurate implementations over providing faster but less accurate implementations.

[sivachandra]

libc/src/math/generic/CMakeLists.txt

503 ↗

(On Diff #392768)

So I learn't that the fma related builtins require this option. The correct way to do this would be to do it this way:

1. Make a separate helper library of x86 builtin calls.
2. List -mfma as an INTERFACE compile option for that library.

Then targets depending on the helper library will automatically get the -mfma compile option.

[lntue]

libc/test/src/math/expm1f_test.cpp

111

We're going to have a correctly rounded version for expm1f soon so I'm not worried about this regression yet.

[lntue]

[libc] Optimize PolyEval for x86-64 with degree 3-6 polynomials.

[lntue]

[libc] Optimize PolyEval for x86-64 with degree 3-6 polynomials.

[lntue]

[libc] Optimize PolyEval for x86-64 with degree 3-6 polynomials.

[lntue]

[libc] Optimize PolyEval for x86-64 with degree 3-6 polynomials.

[lntue]

[libc] Optimize PolyEval for x86-64 with degree 3-6 polynomials.

[sivachandra]

Is the subject line and description correct?

libc/src/__support/FPUtil/PolyEval.h

40

Can these be in just one file say, PolyEvalSpecializations.h?

libc/test/src/math/expm1f_test.cpp

111

Do we know why the accuracy dropped? If yes, can we add the reason to the commit description?

[lntue]

[libc] Use intrinsic for x86-64 fma and optimize PolyEval for x86-64 with degree 3 & 5 polynomials.

[lntue]

libc/test/src/math/expm1f_test.cpp

111

Added a possible reason to the patch's summary.

[lntue]

[libc] Use intrinsics for x86-64 fma and optimize PolyEval for x86-64 with degree 3 & 5 polynomials.

libc/src/__support/FPUtil/PolyEval.h

 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 fma(x, polyeval(x, a...), a0);
 }
 
 } // namespace fputil
 } // namespace __llvm_libc
 
+#ifdef LLVM_LIBC_ARCH_X86_64
+
+#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>

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

 
 #include "src/__support/architectures.h"
 
 #if !defined(LLVM_LIBC_ARCH_X86)
 #error "Invalid include"
 #endif
 
 #include "src/__support/CPP/TypeTraits.h"
+#include <immintrin.h>
 
 namespace __llvm_libc {
 namespace fputil {
 
 template <typename T>
-static inline cpp::EnableIfType<cpp::IsSame<T, float>::Value, T> fma(T x, T y,
-                                                                     T z) {
-  float result = x;
-  __asm__ __volatile__("vfmadd213ss %x2, %x1, %x0"
-                       : "+x"(result)
-                       : "x"(y), "x"(z));
+__attribute__((target(
+    "fma"))) static inline cpp::EnableIfType<cpp::IsSame<T, float>::Value, T>
+fma(T x, T y, T z) {
+  float result;
+  __m128 xmm = _mm_load_ss(&x);
+  __m128 ymm = _mm_load_ss(&y);
+  __m128 zmm = _mm_load_ss(&z);
+  __m128 r = _mm_fmadd_ss(xmm, ymm, zmm);
+  _mm_store_ss(&result, r);
   return result;
 }
 
 template <typename T>
-static inline cpp::EnableIfType<cpp::IsSame<T, double>::Value, T> fma(T x, T y,
-                                                                      T z) {
-  double result = x;
-  __asm__ __volatile__("vfmadd213sd %x2, %x1, %x0"
-                       : "+x"(result)
-                       : "x"(y), "x"(z));
+__attribute__((target(
+    "fma"))) static inline cpp::EnableIfType<cpp::IsSame<T, double>::Value, T>
+fma(T x, T y, T z) {
+  double result;
+  __m128d xmm = _mm_load_sd(&x);
+  __m128d ymm = _mm_load_sd(&y);
+  __m128d zmm = _mm_load_sd(&z);
+  __m128d r = _mm_fmadd_sd(xmm, ymm, zmm);
+  _mm_store_sd(&result, r);
   return result;
 }
 
 } // namespace fputil
 } // namespace __llvm_libc
 
 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_FMA_H

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

+//===-- Optimized PolyEval implementations for x86_64 --------- 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_X86_64_POLYEVAL_H
+#define LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_POLYEVAL_H
+
+#include "src/__support/architectures.h"
+
+#if !defined(LLVM_LIBC_ARCH_X86_64)
+#error "Invalid include"
+#endif
+
+#include <immintrin.h>
+
+namespace __llvm_libc {
+namespace fputil {
+
+// Cubic polynomials:
+//   polyeval(x, a0, a1, a2, a3) = a3*x^3 + a2*x^2 + a1*x + a0
+template <>
+__attribute__((target("fma"))) inline float
+polyeval(float x, float a0, float a1, float a2, float a3) {
+  __m128 xmm = _mm_set1_ps(x);
+  __m128 a13 = _mm_set_ps(0.0f, x, a3, a1);
+  __m128 a02 = _mm_set_ps(0.0f, 0.0f, a2, a0);
+  // r = (0, x^2, a3*x + a2, a1*x + a0)
+  __m128 r = _mm_fmadd_ps(a13, xmm, a02);
+  // result = (a3*x + a2) * x^2 + (a1*x + a0)
+  return fma(r[2], r[1], r[0]);
+}
+
+template <>
+__attribute__((target("fma"))) inline double
+polyeval(double x, double a0, double a1, double a2, double a3) {
+  __m256d xmm = _mm256_set1_pd(x);
+  __m256d a13 = _mm256_set_pd(0.0, x, a3, a1);
+  __m256d a02 = _mm256_set_pd(0.0, 0.0, a2, a0);
+  // r = (0, x^2, a3*x + a2, a1*x + a0)
+  __m256d r = _mm256_fmadd_pd(a13, xmm, a02);
+  // result = (a3*x + a2) * x^2 + (a1*x + a0)
+  return fma(r[2], r[1], r[0]);
+}
+
+// Quintic polynomials:
+//   polyeval(x, a0, a1, a2, a3, a4, a5) = a5*x^5 + a4*x^4 + a3*x^3 + a2*x^2 +
+//                                         + a1*x + a0
+template <>
+__attribute__((target("fma"))) inline float
+polyeval(float x, float a0, float a1, float a2, float a3, float a4, float a5) {
+  __m128 xmm = _mm_set1_ps(x);
+  __m128 a25 = _mm_set_ps(0.0f, x, a5, a2);
+  __m128 a14 = _mm_set_ps(0.0f, 0.0f, a4, a1);
+  __m128 a03 = _mm_set_ps(0.0f, 0.0f, a3, a0);
+  // r1 = (0, x^2, a5*x + a4, a2*x + a1)
+  __m128 r1 = _mm_fmadd_ps(a25, xmm, a14);
+  // r2 = (0, x^3, (a5*x + a4)*x + a3, (a2*x + a1)*x + a0
+  __m128 r2 = _mm_fmadd_ps(r1, xmm, a03);
+  // result = ((a5*x + a4)*x + a3) * x^3 + ((a2*x + a1)*x + a0)
+  return fma(r2[2], r2[1], r2[0]);
+}
+
+template <>
+__attribute__((target("fma"))) inline double
+polyeval(double x, double a0, double a1, double a2, double a3, double a4,
+         double a5) {
+  __m256d xmm = _mm256_set1_pd(x);
+  __m256d a25 = _mm256_set_pd(0.0, x, a5, a2);
+  __m256d a14 = _mm256_set_pd(0.0, 0.0, a4, a1);
+  __m256d a03 = _mm256_set_pd(0.0, 0.0, a3, a0);
+  // r1 = (0, x^2, a5*x + a4, a2*x + a1)
+  __m256d r1 = _mm256_fmadd_pd(a25, xmm, a14);
+  // r2 = (0, x^3, (a5*x + a4)*x + a3, (a2*x + a1)*x + a0
+  __m256d r2 = _mm256_fmadd_pd(r1, xmm, a03);
+  // result = ((a5*x + a4)*x + a3) * x^3 + ((a2*x + a1)*x + a0)
+  return fma(r2[2], r2[1], r2[0]);
+}
+
+} // namespace fputil
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_X86_64_POLYEVAL_H

libc/test/src/math/expm1f_test.cpp

for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) {
     float result = __llvm_libc::expm1f(x);
 
     // If the computation resulted in an error or did not produce valid result
     // in the single-precision floating point range, then ignore comparing with
     // MPFR result as MPFR can still produce valid results because of its
     // wider precision.
     if (isnan(result) || isinf(result) || errno != 0)
       continue;
-    ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.5);
+    ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 2.2);
   }
 }