libc/src/__support/mathvec/expf_utils.h

//===----------------------------------------------------------------------===//
//
// 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
//
//===----------------------------------------------------------------------===//
///
/// \file
/// This file contains utility functions for single-precision SIMD exp/2/10.
///
//===----------------------------------------------------------------------===//

#ifndef LLVM_LIBC_SRC___SUPPORT_MATHVEC_EXP_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATHVEC_EXP_UTILS_H

#include "src/__support/CPP/simd.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/mathvec/common_constants.h"

namespace LIBC_NAMESPACE_DECL {

namespace mathvec {

template <size_t N>
LIBC_INLINE static cpp::simd<double, N> exp_lookup(cpp::simd<uint64_t, N> u) {
  cpp::simd<uint64_t, N> index = u & cpp::simd<uint64_t, N>(0x3f);
  cpp::simd<uint64_t, N> mantissa =
      cpp::gather<cpp::simd<uint64_t, N>>(true, index, EXP_MANTISSA);
  cpp::simd<uint64_t, N> exponent = (u >> 6) << 52;
  cpp::simd<uint64_t, N> result = mantissa | exponent;
  return cpp::bit_cast<cpp::simd<double, N>>(result);
}

template <size_t N>
LIBC_INLINE static cpp::simd<double, N> eval_exp(cpp::simd<double, N> r,
                                                 cpp::simd<double, N> z) {
  // Coefficients of exp approximation, generated by Sollya with:
  // poly = 1 + x;
  // for i from 2 to 5 do {
  //   r = remez(exp(x)-poly(x), 5-i, [-log(2)/128;log(2)/128], x^i, 1e-10);
  //   c = coeff(roundcoefficients(r, [|D ...|]), 0);
  //   poly = poly + x^i*c;
  //   c;
  // };
  constexpr cpp::simd<double, N> c0 = 0x1.fffffffffdbcep-2;
  constexpr cpp::simd<double, N> c1 = 0x1.55555555543c2p-3;
  constexpr cpp::simd<double, N> c2 = 0x1.555573c64f2e3p-5;
  constexpr cpp::simd<double, N> c3 = 0x1.111126b4eff73p-7;

  // y = exp(r) - 1 ~= r + C0 r^2 + C1 r^3 + C2 r^4 + C3 r^5.
  cpp::simd<double, N> r2 = r * r;
  cpp::simd<double, N> p01 = cpp::multiply_add(c1, r, c0);
  cpp::simd<double, N> p23 = cpp::multiply_add(c3, r, c2);
  cpp::simd<double, N> p04 = cpp::multiply_add(p23, r2, p01);
  cpp::simd<double, N> y = cpp::multiply_add(p04, r2, r);

  // Table lookup for 2^n, where n is a multiple of 1/64
  cpp::simd<uint64_t, N> u = cpp::bit_cast<cpp::simd<uint64_t, N>>(z);
  cpp::simd<double, N> s = exp_lookup(u);

  // e^x = 2^n * exp(r)
  // Since y = exp(r) - 1, e^x = 2^n * (1 + y)
  // Or as an FMA: e^x = 2^n + (2^n * y)
  return cpp::multiply_add(y, s, s);
}

} // namespace mathvec

} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATHVEC_EXP_UTILS_H