xref: /aosp_15_r20/external/llvm-libc/src/math/generic/expxf16.h (revision 71db0c75aadcf003ffe3238005f61d7618a3fead)

//===-- Common utilities for half-precision exponential functions ---------===//
//
// 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_MATH_GENERIC_EXPXF16_H
#define LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H

#include "src/__support/CPP/array.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/cast.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/macros/attributes.h"
#include "src/__support/macros/config.h"
#include <stdint.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;
};

LIBC_INLINE 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};
}

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > for i from 0 to 7 do printsingle(round(2^(i * 2^-3), SG, RN));
constexpr cpp::array<uint32_t, 8> EXP2_MID_BITS = {
    0x3f80'0000U, 0x3f8b'95c2U, 0x3f98'37f0U, 0x3fa5'fed7U,
    0x3fb5'04f3U, 0x3fc5'672aU, 0x3fd7'44fdU, 0x3fea'c0c7U,
};

LIBC_INLINE ExpRangeReduction exp2_range_reduction(float16 x) {
  // For -25 < x < 16, to compute 2^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,
  //   2^x = 2^(hi + mid + lo) = 2^hi * 2^mid * 2^lo.
  // We store 2^mid in the lookup table EXP2_MID_BITS, and compute 2^hi * 2^mid
  // by adding hi to the exponent field of 2^mid.  2^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);
  unsigned x_hi = static_cast<unsigned>(x_hi_mid) >> 3;
  unsigned x_mid = static_cast<unsigned>(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);

  uint32_t exp2_hi_mid_bits =
      EXP2_MID_BITS[x_mid] +
      static_cast<uint32_t>(x_hi << fputil::FPBits<float>::FRACTION_LEN);
  float exp2_hi_mid = fputil::FPBits<float>(exp2_hi_mid_bits).get_val();
  // Degree-3 minimax polynomial generated by Sollya with the following
  // commands:
  //   > display = hexadecimal;
  //   > P = fpminimax((2^x - 1)/x, 2, [|SG...|], [-2^-4, 2^-4]);
  //   > 1 + x * P;
  float exp2_lo = fputil::polyeval(lo, 0x1p+0f, 0x1.62e43p-1f, 0x1.ec0aa6p-3f,
                                   0x1.c6b4a6p-5f);
  return {exp2_hi_mid, exp2_lo};
}

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > round(log2(10), SG, RN);
static constexpr float LOG2F_10 = 0x1.a934fp+1f;

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > round(log10(2), SG, RN);
static constexpr float LOG10F_2 = 0x1.344136p-2f;

LIBC_INLINE ExpRangeReduction exp10_range_reduction(float16 x) {
  // For -8 < x < 5, to compute 10^x, we perform the following range reduction:
  // find hi, mid, lo, such that:
  //   x = (hi + mid) * log2(10) + 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,
  //   10^x = 10^(hi + mid + lo) = 2^((hi + mid) * log2(10)) + 10^lo
  // We store 2^mid in the lookup table EXP2_MID_BITS, and compute 2^hi * 2^mid
  // by adding hi to the exponent field of 2^mid.  10^lo is computed using a
  // degree-4 minimax polynomial generated by Sollya.

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

  uint32_t exp2_hi_mid_bits =
      EXP2_MID_BITS[x_mid] +
      static_cast<uint32_t>(x_hi << fputil::FPBits<float>::FRACTION_LEN);
  float exp2_hi_mid = fputil::FPBits<float>(exp2_hi_mid_bits).get_val();
  // Degree-4 minimax polynomial generated by Sollya with the following
  // commands:
  //   > display = hexadecimal;
  //   > P = fpminimax((10^x - 1)/x, 3, [|SG...|], [-2^-4, 2^-4]);
  //   > 1 + x * P;
  float exp10_lo = fputil::polyeval(lo, 0x1p+0f, 0x1.26bb14p+1f, 0x1.53526p+1f,
                                    0x1.04b434p+1f, 0x1.2bcf9ep+0f);
  return {exp2_hi_mid, exp10_lo};
}

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > round(log2(exp(1)), SG, RN);
static constexpr float LOG2F_E = 0x1.715476p+0f;

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > round(log(2), SG, RN);
static constexpr float LOGF_2 = 0x1.62e43p-1f;

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > for i from 0 to 31 do printsingle(round(2^(i * 2^-5), SG, RN));
static constexpr cpp::array<uint32_t, 32> EXP2_MID_5_BITS = {
    0x3f80'0000U, 0x3f82'cd87U, 0x3f85'aac3U, 0x3f88'980fU, 0x3f8b'95c2U,
    0x3f8e'a43aU, 0x3f91'c3d3U, 0x3f94'f4f0U, 0x3f98'37f0U, 0x3f9b'8d3aU,
    0x3f9e'f532U, 0x3fa2'7043U, 0x3fa5'fed7U, 0x3fa9'a15bU, 0x3fad'583fU,
    0x3fb1'23f6U, 0x3fb5'04f3U, 0x3fb8'fbafU, 0x3fbd'08a4U, 0x3fc1'2c4dU,
    0x3fc5'672aU, 0x3fc9'b9beU, 0x3fce'248cU, 0x3fd2'a81eU, 0x3fd7'44fdU,
    0x3fdb'fbb8U, 0x3fe0'ccdfU, 0x3fe5'b907U, 0x3fea'c0c7U, 0x3fef'e4baU,
    0x3ff5'257dU, 0x3ffa'83b3U,
};

// This function correctly calculates sinh(x) and cosh(x) by calculating exp(x)
// and exp(-x) simultaneously.
// To compute e^x, we perform the following range reduction:
// find hi, mid, lo such that:
//   x = (hi + mid) * log(2) + lo, in which
//     hi is an integer,
//     0 <= mid * 2^5 < 32 is an integer
//     -2^(-5) <= lo * log2(e) <= 2^-5.
// In particular,
//   hi + mid = round(x * log2(e) * 2^5) * 2^(-5).
// Then,
//   e^x = 2^(hi + mid) * e^lo = 2^hi * 2^mid * e^lo.
// We store 2^mid in the lookup table EXP2_MID_5_BITS, and compute 2^hi * 2^mid
// by adding hi to the exponent field of 2^mid.
// e^lo is computed using a degree-3 minimax polynomial generated by Sollya:
//   e^lo ~ P(lo)
//        = 1 + lo + c2 * lo^2 + ... + c5 * lo^5
//        = (1 + c2*lo^2 + c4*lo^4) + lo * (1 + c3*lo^2 + c5*lo^4)
//        = P_even + lo * P_odd
// To compute e^(-x), notice that:
//   e^(-x) = 2^(-(hi + mid)) * e^(-lo)
//          ~ 2^(-(hi + mid)) * P(-lo)
//          = 2^(-(hi + mid)) * (P_even - lo * P_odd)
// So:
//   sinh(x) = (e^x - e^(-x)) / 2
//           ~ 0.5 * (2^(hi + mid) * (P_even + lo * P_odd) -
//                    2^(-(hi + mid)) * (P_even - lo * P_odd))
//           = 0.5 * (P_even * (2^(hi + mid) - 2^(-(hi + mid))) +
//                    lo * P_odd * (2^(hi + mid) + 2^(-(hi + mid))))
// And similarly:
//   cosh(x) = (e^x + e^(-x)) / 2
//           ~ 0.5 * (P_even * (2^(hi + mid) + 2^(-(hi + mid))) +
//                    lo * P_odd * (2^(hi + mid) - 2^(-(hi + mid))))
// The main point of these formulas is that the expensive part of calculating
// the polynomials approximating lower parts of e^x and e^(-x) is shared and
// only done once.
template <bool IsSinh> LIBC_INLINE float16 eval_sinh_or_cosh(float16 x) {
  float xf = x;
  float kf = fputil::nearest_integer(xf * (LOG2F_E * 0x1.0p+5f));
  int x_hi_mid_p = static_cast<int>(kf);
  int x_hi_mid_m = -x_hi_mid_p;

  unsigned x_hi_p = static_cast<unsigned>(x_hi_mid_p) >> 5;
  unsigned x_hi_m = static_cast<unsigned>(x_hi_mid_m) >> 5;
  unsigned x_mid_p = static_cast<unsigned>(x_hi_mid_p) & 0x1f;
  unsigned x_mid_m = static_cast<unsigned>(x_hi_mid_m) & 0x1f;

  uint32_t exp2_hi_mid_bits_p =
      EXP2_MID_5_BITS[x_mid_p] +
      static_cast<uint32_t>(x_hi_p << fputil::FPBits<float>::FRACTION_LEN);
  uint32_t exp2_hi_mid_bits_m =
      EXP2_MID_5_BITS[x_mid_m] +
      static_cast<uint32_t>(x_hi_m << fputil::FPBits<float>::FRACTION_LEN);
  // exp2_hi_mid_p = 2^(hi + mid)
  float exp2_hi_mid_p = fputil::FPBits<float>(exp2_hi_mid_bits_p).get_val();
  // exp2_hi_mid_m = 2^(-(hi + mid))
  float exp2_hi_mid_m = fputil::FPBits<float>(exp2_hi_mid_bits_m).get_val();

  // exp2_hi_mid_sum = 2^(hi + mid) + 2^(-(hi + mid))
  float exp2_hi_mid_sum = exp2_hi_mid_p + exp2_hi_mid_m;
  // exp2_hi_mid_diff = 2^(hi + mid) - 2^(-(hi + mid))
  float exp2_hi_mid_diff = exp2_hi_mid_p - exp2_hi_mid_m;

  // lo = x - (hi + mid) = round(x * log2(e) * 2^5) * log(2) * (-2^(-5)) + x
  float lo = fputil::multiply_add(kf, LOGF_2 * -0x1.0p-5f, xf);
  float lo_sq = lo * lo;

  // Degree-3 minimax polynomial generated by Sollya with the following
  // commands:
  //   > display = hexadecimal;
  //   > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-5, 2^-5]);
  //   > 1 + x * P;
  constexpr cpp::array<float, 4> COEFFS = {0x1p+0f, 0x1p+0f, 0x1.0004p-1f,
                                           0x1.555778p-3f};
  float half_p_odd =
      fputil::polyeval(lo_sq, COEFFS[1] * 0.5f, COEFFS[3] * 0.5f);
  float half_p_even =
      fputil::polyeval(lo_sq, COEFFS[0] * 0.5f, COEFFS[2] * 0.5f);

  // sinh(x) = lo * (0.5 * P_odd * (2^(hi + mid) + 2^(-(hi + mid)))) +
  //                (0.5 * P_even * (2^(hi + mid) - 2^(-(hi + mid))))
  if constexpr (IsSinh)
    return fputil::cast<float16>(fputil::multiply_add(
        lo, half_p_odd * exp2_hi_mid_sum, half_p_even * exp2_hi_mid_diff));
  // cosh(x) = lo * (0.5 * P_odd * (2^(hi + mid) - 2^(-(hi + mid)))) +
  //                (0.5 * P_even * (2^(hi + mid) + 2^(-(hi + mid))))
  return fputil::cast<float16>(fputil::multiply_add(
      lo, half_p_odd * exp2_hi_mid_diff, half_p_even * exp2_hi_mid_sum));
}

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > for i from 0 to 31 do print(round(log(1 + i * 2^-5), SG, RN));
constexpr cpp::array<float, 32> LOGF_F = {
    0x0p+0f,        0x1.f829bp-6f,  0x1.f0a30cp-5f, 0x1.6f0d28p-4f,
    0x1.e27076p-4f, 0x1.29553p-3f,  0x1.5ff308p-3f, 0x1.9525aap-3f,
    0x1.c8ff7cp-3f, 0x1.fb9186p-3f, 0x1.1675cap-2f, 0x1.2e8e2cp-2f,
    0x1.4618bcp-2f, 0x1.5d1bdcp-2f, 0x1.739d8p-2f,  0x1.89a338p-2f,
    0x1.9f323ep-2f, 0x1.b44f78p-2f, 0x1.c8ff7cp-2f, 0x1.dd46ap-2f,
    0x1.f128f6p-2f, 0x1.02552ap-1f, 0x1.0be72ep-1f, 0x1.154c3ep-1f,
    0x1.1e85f6p-1f, 0x1.2795e2p-1f, 0x1.307d74p-1f, 0x1.393e0ep-1f,
    0x1.41d8fep-1f, 0x1.4a4f86p-1f, 0x1.52a2d2p-1f, 0x1.5ad404p-1f,
};

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > for i from 0 to 31 do print(round(log2(1 + i * 2^-5), SG, RN));
constexpr cpp::array<float, 32> LOG2F_F = {
    0x0p+0f,        0x1.6bad38p-5f, 0x1.663f7p-4f,  0x1.08c588p-3f,
    0x1.5c01a4p-3f, 0x1.acf5e2p-3f, 0x1.fbc16cp-3f, 0x1.24407ap-2f,
    0x1.49a784p-2f, 0x1.6e221cp-2f, 0x1.91bba8p-2f, 0x1.b47ecp-2f,
    0x1.d6753ep-2f, 0x1.f7a856p-2f, 0x1.0c105p-1f,  0x1.1bf312p-1f,
    0x1.2b8034p-1f, 0x1.3abb4p-1f,  0x1.49a784p-1f, 0x1.584822p-1f,
    0x1.66a008p-1f, 0x1.74b1fep-1f, 0x1.82809ep-1f, 0x1.900e62p-1f,
    0x1.9d5dap-1f,  0x1.aa709p-1f,  0x1.b74948p-1f, 0x1.c3e9cap-1f,
    0x1.d053f6p-1f, 0x1.dc899ap-1f, 0x1.e88c6cp-1f, 0x1.f45e08p-1f,
};

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > for i from 0 to 31 do print(round(log10(1 + i * 2^-5), SG, RN));
constexpr cpp::array<float, 32> LOG10F_F = {
    0x0p+0f,        0x1.b5e908p-7f, 0x1.af5f92p-6f, 0x1.3ed11ap-5f,
    0x1.a30a9ep-5f, 0x1.02428cp-4f, 0x1.31b306p-4f, 0x1.5fe804p-4f,
    0x1.8cf184p-4f, 0x1.b8de4ep-4f, 0x1.e3bc1ap-4f, 0x1.06cbd6p-3f,
    0x1.1b3e72p-3f, 0x1.2f3b6ap-3f, 0x1.42c7e8p-3f, 0x1.55e8c6p-3f,
    0x1.68a288p-3f, 0x1.7af974p-3f, 0x1.8cf184p-3f, 0x1.9e8e7cp-3f,
    0x1.afd3e4p-3f, 0x1.c0c514p-3f, 0x1.d1653p-3f,  0x1.e1b734p-3f,
    0x1.f1bdeep-3f, 0x1.00be06p-2f, 0x1.087a08p-2f, 0x1.101432p-2f,
    0x1.178da6p-2f, 0x1.1ee778p-2f, 0x1.2622bp-2f,  0x1.2d404cp-2f,
};

// Generated by Sollya with the following commands:
//   > display = hexadecimal;
//   > for i from 0 to 31 do print(round(1 / (1 + i * 2^-5), SG, RN));
constexpr cpp::array<float, 32> ONE_OVER_F_F = {
    0x1p+0f,        0x1.f07c2p-1f,  0x1.e1e1e2p-1f, 0x1.d41d42p-1f,
    0x1.c71c72p-1f, 0x1.bacf92p-1f, 0x1.af286cp-1f, 0x1.a41a42p-1f,
    0x1.99999ap-1f, 0x1.8f9c18p-1f, 0x1.861862p-1f, 0x1.7d05f4p-1f,
    0x1.745d18p-1f, 0x1.6c16c2p-1f, 0x1.642c86p-1f, 0x1.5c9882p-1f,
    0x1.555556p-1f, 0x1.4e5e0ap-1f, 0x1.47ae14p-1f, 0x1.414142p-1f,
    0x1.3b13b2p-1f, 0x1.3521dp-1f,  0x1.2f684cp-1f, 0x1.29e412p-1f,
    0x1.24924ap-1f, 0x1.1f7048p-1f, 0x1.1a7b96p-1f, 0x1.15b1e6p-1f,
    0x1.111112p-1f, 0x1.0c9714p-1f, 0x1.08421p-1f,  0x1.041042p-1f,
};

} // namespace LIBC_NAMESPACE_DECL

#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H