libc/src/__support/math/expf16_utils.h

//===-- Common utils for expf16 functions -----------------------*- 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_EXPF16_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_UTILS_H

#include "include/llvm-libc-macros/float16-macros.h"

#ifdef LIBC_TYPES_HAS_FLOAT16

#include "src/__support/CPP/array.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/macros/properties/types.h"

namespace LIBC_NAMESPACE_DECL {

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > for i from -18 to 12 do print(round(exp(i), SG, RN));
static constexpr cpp::array<float, 31> EXP_HI = {
    0x1.05a628p-26f, 0x1.639e32p-25f, 0x1.e355bcp-24f, 0x1.4875cap-22f,
    0x1.be6c7p-21f,  0x1.2f6054p-19f, 0x1.9c54c4p-18f, 0x1.183542p-16f,
    0x1.7cd79cp-15f, 0x1.02cf22p-13f, 0x1.5fc21p-12f,  0x1.de16bap-11f,
    0x1.44e52p-9f,   0x1.b993fep-8f,  0x1.2c155cp-6f,  0x1.97db0cp-5f,
    0x1.152aaap-3f,  0x1.78b564p-2f,  0x1p+0f,         0x1.5bf0a8p+1f,
    0x1.d8e64cp+2f,  0x1.415e5cp+4f,  0x1.b4c902p+5f,  0x1.28d38ap+7f,
    0x1.936dc6p+8f,  0x1.122886p+10f, 0x1.749ea8p+11f, 0x1.fa7158p+12f,
    0x1.5829dcp+14f, 0x1.d3c448p+15f, 0x1.3de166p+17f,
};

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > for i from 0 to 7 do print(round(exp(i * 2^-3), SG, RN));
static constexpr cpp::array<float, 8> EXP_MID = {
    0x1p+0f,        0x1.221604p+0f, 0x1.48b5e4p+0f, 0x1.747a52p+0f,
    0x1.a61298p+0f, 0x1.de455ep+0f, 0x1.0ef9dcp+1f, 0x1.330e58p+1f,
};

struct ExpRangeReduction {
  float exp_hi_mid;
  float exp_lo;
};

static constexpr ExpRangeReduction exp_range_reduction(float16 x) {
  // For -18 < x < 12, to compute exp(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^3 is an integer,
  //     -2^(-4) <= lo < 2^(-4).
  // In particular,
  //   hi + mid = round(x * 2^3) * 2^(-3).
  // Then,
  //   exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo).
  // We store exp(hi) and exp(mid) in the lookup tables EXP_HI and EXP_MID
  // respectively.  exp(lo) is computed using a degree-3 minimax polynomial
  // generated by Sollya.

  float xf = x;
  float kf = fputil::nearest_integer(xf * 0x1.0p+3f);
  int x_hi_mid = static_cast<int>(kf);
  int x_hi = x_hi_mid >> 3;
  int x_mid = x_hi_mid & 0x7;
  // lo = x - (hi + mid) = round(x * 2^3) * (-2^(-3)) + x
  float lo = fputil::multiply_add(kf, -0x1.0p-3f, xf);

  float exp_hi = EXP_HI[x_hi + 18];
  float exp_mid = EXP_MID[x_mid];
  // Degree-3 minimax polynomial generated by Sollya with the following
  // commands:
  //   > display = hexadecimal;
  //   > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-4, 2^-4]);
  //   > 1 + x * P;
  float exp_lo =
      fputil::polyeval(lo, 0x1p+0f, 0x1p+0f, 0x1.001p-1f, 0x1.555ddep-3f);
  return {exp_hi * exp_mid, exp_lo};
}

} // namespace LIBC_NAMESPACE_DECL

#endif // LIBC_TYPES_HAS_FLOAT16

#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_UTILS_H