91 lines · c
1//===-- Common utils for expf16 functions -----------------------*- C++ -*-===//2//3// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.4// See https://llvm.org/LICENSE.txt for license information.5// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_UTILS_H10#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_UTILS_H11 12#include "include/llvm-libc-macros/float16-macros.h"13 14#ifdef LIBC_TYPES_HAS_FLOAT1615 16#include "src/__support/CPP/array.h"17#include "src/__support/FPUtil/PolyEval.h"18#include "src/__support/FPUtil/nearest_integer.h"19#include "src/__support/macros/properties/types.h"20 21namespace LIBC_NAMESPACE_DECL {22 23// Generated by Sollya with the following commands:24// > display = hexadecimal;25// > for i from -18 to 12 do print(round(exp(i), SG, RN));26static constexpr cpp::array<float, 31> EXP_HI = {27 0x1.05a628p-26f, 0x1.639e32p-25f, 0x1.e355bcp-24f, 0x1.4875cap-22f,28 0x1.be6c7p-21f, 0x1.2f6054p-19f, 0x1.9c54c4p-18f, 0x1.183542p-16f,29 0x1.7cd79cp-15f, 0x1.02cf22p-13f, 0x1.5fc21p-12f, 0x1.de16bap-11f,30 0x1.44e52p-9f, 0x1.b993fep-8f, 0x1.2c155cp-6f, 0x1.97db0cp-5f,31 0x1.152aaap-3f, 0x1.78b564p-2f, 0x1p+0f, 0x1.5bf0a8p+1f,32 0x1.d8e64cp+2f, 0x1.415e5cp+4f, 0x1.b4c902p+5f, 0x1.28d38ap+7f,33 0x1.936dc6p+8f, 0x1.122886p+10f, 0x1.749ea8p+11f, 0x1.fa7158p+12f,34 0x1.5829dcp+14f, 0x1.d3c448p+15f, 0x1.3de166p+17f,35};36 37// Generated by Sollya with the following commands:38// > display = hexadecimal;39// > for i from 0 to 7 do print(round(exp(i * 2^-3), SG, RN));40static constexpr cpp::array<float, 8> EXP_MID = {41 0x1p+0f, 0x1.221604p+0f, 0x1.48b5e4p+0f, 0x1.747a52p+0f,42 0x1.a61298p+0f, 0x1.de455ep+0f, 0x1.0ef9dcp+1f, 0x1.330e58p+1f,43};44 45struct ExpRangeReduction {46 float exp_hi_mid;47 float exp_lo;48};49 50[[maybe_unused]] LIBC_INLINE static ExpRangeReduction51exp_range_reduction(float16 x) {52 // For -18 < x < 12, to compute exp(x), we perform the following range53 // reduction: find hi, mid, lo, such that:54 // x = hi + mid + lo, in which55 // hi is an integer,56 // mid * 2^3 is an integer,57 // -2^(-4) <= lo < 2^(-4).58 // In particular,59 // hi + mid = round(x * 2^3) * 2^(-3).60 // Then,61 // exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo).62 // We store exp(hi) and exp(mid) in the lookup tables EXP_HI and EXP_MID63 // respectively. exp(lo) is computed using a degree-3 minimax polynomial64 // generated by Sollya.65 66 float xf = x;67 float kf = fputil::nearest_integer(xf * 0x1.0p+3f);68 int x_hi_mid = static_cast<int>(kf);69 int x_hi = x_hi_mid >> 3;70 int x_mid = x_hi_mid & 0x7;71 // lo = x - (hi + mid) = round(x * 2^3) * (-2^(-3)) + x72 float lo = fputil::multiply_add(kf, -0x1.0p-3f, xf);73 74 float exp_hi = EXP_HI[x_hi + 18];75 float exp_mid = EXP_MID[x_mid];76 // Degree-3 minimax polynomial generated by Sollya with the following77 // commands:78 // > display = hexadecimal;79 // > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-4, 2^-4]);80 // > 1 + x * P;81 float exp_lo =82 fputil::polyeval(lo, 0x1p+0f, 0x1p+0f, 0x1.001p-1f, 0x1.555ddep-3f);83 return {exp_hi * exp_mid, exp_lo};84}85 86} // namespace LIBC_NAMESPACE_DECL87 88#endif // LIBC_TYPES_HAS_FLOAT1689 90#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPF16_UTILS_H91