[libc-commits] [libc] [libc][mathvec] Vectorise logf (PR #231190)

[Joseph Huber via libc-commits]

================
@@ -7,30 +7,118 @@
 //===----------------------------------------------------------------------===//
 ///
 /// \file
-/// This file contains the implemenation for single-precision SIMD log.
+/// This file contains the implementation for single-precision SIMD log.
 ///
 //===----------------------------------------------------------------------===//
 
 #ifndef LLVM_LIBC_SRC___SUPPORT_MATHVEC_LOGF_H
 #define LLVM_LIBC_SRC___SUPPORT_MATHVEC_LOGF_H
 
+#include "src/__support/CPP/bit.h"
 #include "src/__support/CPP/simd.h"
-#define LIBC_MATH_HAS_NO_ERRNO
-#define LIBC_MATH_HAS_NO_EXCEPT
-#define LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
-#include "src/__support/math/logf.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/common_constants.h"
+#include "src/__support/mathvec/common_constants.h"
 
 namespace LIBC_NAMESPACE_DECL {
 
 namespace mathvec {
 
 template <size_t N>
-LIBC_INLINE cpp::simd<float, N> logf(cpp::simd<float, N> x) {
-  return cpp::map(x, [](float a) { return math::logf(a); });
+LIBC_INLINE cpp::simd<float, N>
+inline_logf(cpp::simd<uint32_t, N> ux, cpp::simd<int32_t, N> exponent_bias) {
+  using namespace common_constants_internal;
+
+  // Extract the 7 leading bits of the mantissa
+  cpp::simd<uint32_t, N> index = (ux >> 16) & 0x7fU;
+
+  // Compute unbiased exponent
+  cpp::simd<int32_t, N> m = cpp::simd_cast<int32_t>(ux >> 23);
+  m += exponent_bias;
+
+  // Sets the exponent such that 1 <= u < 2
+  cpp::simd<float, N> u =
+      cpp::bit_cast<cpp::simd<float, N>>((ux & 0x007f'ffffU) | 0x3f80'0000U);
+
+  cpp::simd<double, N> log_r =
+      cpp::gather<cpp::simd<double, N>>(true, index, LOG_TABLE);
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT
+  cpp::simd<float, N> r = cpp::gather<cpp::simd<float, N>>(true, index, R);
+  cpp::simd<double, N> v = cpp::simd_cast<double>(
+      cpp::multiply_add(u, r, cpp::simd<float, N>(-1.0f)));
+#else
+  cpp::simd<double, N> r = cpp::gather<cpp::simd<double, N>>(true, index, RD);
+  cpp::simd<double, N> v = cpp::simd_cast<double>(u) * r - 1.0;
+#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
+
+  // Degree-5 polynomial approximation of log generated by Sollya with:
+  // > P = fpminimax(log(1 + x)/x, 4, [|1, D...|], [-2^-8, 2^-7]);
+  constexpr cpp::simd<double, N> c0 = -0x1.000000000fe63p-1;
+  constexpr cpp::simd<double, N> c1 = 0x1.555556e963c16p-2;
+  constexpr cpp::simd<double, N> c2 = -0x1.000028dedf986p-2;
+  constexpr cpp::simd<double, N> c3 = 0x1.966681bfda7f7p-3;
+  constexpr cpp::simd<double, N> log2 = 0x1.62e42fefa39efp-1;
+
+  // Evaluate log(1 + v) in double precision.
+  cpp::simd<double, N> v_squared = v * v;
+  cpp::simd<double, N> p2 = cpp::multiply_add(v, c3, c2);
+  cpp::simd<double, N> p1 = cpp::multiply_add(v, c1, c0);
+
+  // Accumulate the exponent independently of the polynomial's FMA chain.
+  cpp::simd<double, N> exponent = cpp::simd_cast<double>(m);
+  cpp::simd<double, N> p0 = cpp::multiply_add(exponent, log2, log_r) + v;
+  cpp::simd<double, N> y_d =
+      cpp::multiply_add(v_squared, cpp::multiply_add(v_squared, p2, p1), p0);
+  cpp::simd<float, N> y = cpp::simd_cast<float>(y_d);
+
+#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
----------------
jhuber6 wrote:

Could we combine these with `defined()`?

https://github.com/llvm/llvm-project/pull/231190