libc/src/__support/math/powf_small_tables.h

//===-- Implementation header for powf using less memory --------*- 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_MATH_POWF_SMALL_TABLES_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_POWF_SMALL_TABLES_H

#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"

namespace LIBC_NAMESPACE_DECL {

namespace math {

namespace powf_internal {

LIBC_INLINE LIBC_CONSTEXPR float powf_small_tables(float x, int ex,
                                                   uint64_t sign, float y) {
  using FloatBits = fputil::FPBits<float>;
  using DoubleBits = fputil::FPBits<double>;

  constexpr double ONE_OVER_SQRT2 = 0x1.6a09e667f3bcdp-1;

  // x^y = 2^( y * log2(x) )
  //     = 2^( y * ( e_x + log2(m_x) ) )
  // First we compute log2(x) = e_x + log2(m_x)
  uint32_t x_u = FloatBits(x).uintval();

  double yd = static_cast<double>(y);

  // Extract exponent field of x.
  ex += (x_u >> FloatBits::FRACTION_LEN);
  double e_x = static_cast<double>(ex);

  // Add the hidden bit to the mantissa.
  // 1 <= m_x < 2
  uint32_t x_mant = (x_u & FloatBits::FRACTION_MASK);
  double m_x = static_cast<double>(cpp::bit_cast<float>(x_mant | 0x3f800000));
  // Reduce to 1 <= mx <= sqrt(2).
  if (x_mant > 0x0045'04f3) {
    e_x += 0.5;
    m_x *= ONE_OVER_SQRT2;
  }
  // 0 <= dx <= sqrt(2) - 1.
  double dx = m_x - 1.0;

  // Degree-13 polynomial approximation:
  //   dx * P(dx) ~ log2(1 + dx)
  // Generated by Sollya with:
  // > P = fpminimax(log2(1 + x)/x, 13, [|D...|], [0, sqrt(2) - 1]);
  // > dirtyinfnorm((log2(1 + x) - x*P)/log2(1 + x), [0, sqrt(2) - 1]);
  //   0x1.b2d...p-53
  constexpr double LOG2_COEFFS[] = {
      0x1.71547652b82fdp0,   -0x1.71547652b7a2ap-1, 0x1.ec709dc2edfa6p-2,
      -0x1.71547626a9d98p-2, 0x1.2776bf5f6f40ep-2,  -0x1.ec6fbbf289ce3p-3,
      0x1.a60bf904470a7p-3,  -0x1.70ef61b01fc1ep-3, 0x1.45d3270454507p-3,
      -0x1.1c5fc05b06e8fp-3, 0x1.d0f57944a937fp-4,  -0x1.413e22be24d32p-4,
      0x1.3c84b66491ccp-5,   -0x1.3df9cfe5e602ep-7};

  double dx2 = dx * dx;
  double c0 = fputil::multiply_add(dx, LOG2_COEFFS[1], LOG2_COEFFS[0]);
  double c1 = fputil::multiply_add(dx, LOG2_COEFFS[3], LOG2_COEFFS[2]);
  double c2 = fputil::multiply_add(dx, LOG2_COEFFS[5], LOG2_COEFFS[4]);
  double c3 = fputil::multiply_add(dx, LOG2_COEFFS[7], LOG2_COEFFS[6]);
  double c4 = fputil::multiply_add(dx, LOG2_COEFFS[9], LOG2_COEFFS[8]);
  double c5 = fputil::multiply_add(dx, LOG2_COEFFS[11], LOG2_COEFFS[10]);
  double c6 = fputil::multiply_add(dx, LOG2_COEFFS[13], LOG2_COEFFS[12]);

  double dx4 = dx2 * dx2;
  double d0 = fputil::multiply_add(dx2, c1, c0);
  double d1 = fputil::multiply_add(dx2, c3, c2);
  double d2 = fputil::multiply_add(dx2, c5, c4);

  double p = fputil::polyeval(dx4, d0, d1, d2, c6);
  // u ~ y * log2(x).
  double u = yd * fputil::multiply_add(dx, p, e_x);

  double hi = fputil::nearest_integer(u);
  double lo = u - hi;
  int e_hi = static_cast<int>(hi) + DoubleBits::EXP_BIAS;
  double exp_hi = cpp::bit_cast<double>(
      (static_cast<uint64_t>(e_hi) << DoubleBits::FRACTION_LEN) | sign);
  // Degree-6 polynomial approximation P(lo6) ~ 2^(lo6 / 2^6) = 2^(lo).
  // Generated by Sollya with:
  // > P = fpminimax(2^x, 6, [|1, D...|], [-0.5, 0.5]);
  // > dirtyinfnorm(2^x - P, [-0.5, 0.5]);
  // 0x1.5f7...p-29
  constexpr double EXP2_COEFFS[] = {
      0x1.62e430c7b13a8p-1, 0x1.ebfbdd2f82f6fp-3, 0x1.c6aed4f186f34p-5,
      0x1.3b2c96c9aa336p-7, 0x1.5f4553ff53f9p-10, 0x1.4278e5fa9de78p-13};

  double lo2 = lo * lo;
  double f0 = fputil::multiply_add(lo, EXP2_COEFFS[1], EXP2_COEFFS[0]);
  double f1 = fputil::multiply_add(lo, EXP2_COEFFS[3], EXP2_COEFFS[2]);
  double f2 = fputil::multiply_add(lo, EXP2_COEFFS[5], EXP2_COEFFS[4]);

  double pp = fputil::polyeval(lo2, f0, f1, f2);

  double r = fputil::multiply_add(lo, pp, 1.0);

  double result = r * exp_hi;

  return static_cast<float>(result);
}

} // namespace powf_internal
} // namespace math
} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_POWF_SMALL_TABLES_H