llvm-project

1
//===-- Implementation of expk function ----------------------------------===//
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-exception
6
//
7
//===----------------------------------------------------------------------===//
8

9
#include "expk.h"
10
#include "src/__support/CPP/bit.h"
11
#include "src/__support/common.h"
12
#include "src/__support/fixed_point/fx_bits.h"
13

14
namespace LIBC_NAMESPACE {
15

16
namespace {
17

18
// Look up tables for exp(hi) and exp(mid).
19
// Generated with Sollya:
20
// > for i from 0 to 23 do {
21
//     hi = i - 11;
22
//     e_hi = nearestint(exp(hi) * 2^15) * 2^-15;
23
//     print(e_hi, "k,");
24
//   };
25
static constexpr accum EXP_HI[24] = {
26
    0x1p-15k,        0x1p-15k,         0x1p-13k,        0x1.6p-12k,
27
    0x1.ep-11k,      0x1.44p-9k,       0x1.bap-8k,      0x1.2cp-6k,
28
    0x1.97cp-5k,     0x1.153p-3k,      0x1.78b8p-2k,    0x1p0k,
29
    0x1.5bf1p1k,     0x1.d8e68p2k,     0x1.415e6p4k,    0x1.b4c9p5k,
30
    0x1.28d388p7k,   0x1.936dc6p8k,    0x1.1228858p10k, 0x1.749ea7cp11k,
31
    0x1.fa7157cp12k, 0x1.5829dcf8p14k, 0x1.d3c4489p15k, ACCUM_MAX,
32
};
33

34
// Generated with Sollya:
35
// > for i from 0 to 15 do {
36
//     m = i/16 - 0.0625;
37
//     e_m = nearestint(exp(m) * 2^15) * 2^-15;
38
//     print(e_m, "k,");
39
//   };
40
static constexpr accum EXP_MID[16] = {
41
    0x1.e0fcp-1k, 0x1p0k,      0x1.1082p0k, 0x1.2216p0k,
42
    0x1.34ccp0k,  0x1.48b6p0k, 0x1.5deap0k, 0x1.747ap0k,
43
    0x1.8c8p0k,   0x1.a612p0k, 0x1.c14cp0k, 0x1.de46p0k,
44
    0x1.fd1ep0k,  0x1.0efap1k, 0x1.2074p1k, 0x1.330ep1k,
45
};
46

47
} // anonymous namespace
48

49
LLVM_LIBC_FUNCTION(accum, expk, (accum x)) {
50
  using FXRep = fixed_point::FXRep<accum>;
51
  using StorageType = typename FXRep::StorageType;
52
  // Output overflow
53
  // > floor(log(2^16) * 2^15) * 2^-15
54
  if (LIBC_UNLIKELY(x >= 0x1.62e4p3k))
55
    return FXRep::MAX();
56
  // Lower bound where exp(x) -> 0:
57
  //   floor(log(2^-16) * 2^15) * 2^-15
58
  if (LIBC_UNLIKELY(x <= -0x1.62e44p3k))
59
    return FXRep::ZERO();
60

61
  // Current range of x:
62
  //   -0x1.62e4p3 <= x <= 0x1.62e3cp3
63
  // Range reduction:
64
  //   x = hi + mid + lo,
65
  // where:
66
  //   hi is an integer
67
  //   mid * 2^4 is an integer
68
  //   |lo| <= 2^-5.
69
  // Then exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo)
70
  //             ~ exp(hi) * exp(mid) * (1 + lo + lo^2 / 2)
71
  // with relative errors < |lo|^3/2 <= 2^-16.
72
  //   exp(hi) and exp(mid) are extracted from small lookup tables.
73

74
  // Round-to-nearest 1/16, tie-to-(+Int):
75
  constexpr accum ONE_THIRTY_SECOND = 0x1.0p-5k;
76
  // x_rounded = floor(x + 1/16).
77
  accum x_rounded = ((x + ONE_THIRTY_SECOND) >> (FXRep::FRACTION_LEN - 4))
78
                    << (FXRep::FRACTION_LEN - 4);
79
  accum lo = x - x_rounded;
80

81
  // Range of x_rounded:
82
  //   x_rounded >= floor((-0x1.62e4p3 + 0x1.0p-5) * 2^4) * 2^-4
83
  //              = -0x1.62p3 = -11.0625
84
  // To get the indices, we shift the values so that it start with 0.
85
  // Range of indices: 0 <= indices <= 355.
86
  StorageType indices = cpp::bit_cast<StorageType>((x_rounded + 0x1.62p3k) >>
87
                                                   (FXRep::FRACTION_LEN - 4));
88
  // So we have the following relation:
89
  //   indices = (hi + mid + 177/16) * 16
90
  // That implies:
91
  //   hi + mid = indices/16 - 11.0625
92
  // So for lookup tables, we can use the upper 4 bits to get:
93
  //   exp( floor(indices / 16) - 11 )
94
  // and lower 4 bits for:
95
  //   exp( (indices - floor(indices)) - 0.0625 )
96
  accum exp_hi = EXP_HI[indices >> 4];
97
  accum exp_mid = EXP_MID[indices & 0xf];
98
  // exp(x) ~ exp(hi) * exp(mid) * (1 + lo);
99
  accum l1 = 0x1.0p0k + (lo >> 1); // = 1 + lo / 2
100
  accum l2 = 0x1.0p0k + lo * l1;   // = 1 + lo * (1 + lo / 2) = 1 + lo + lo^2/2
101
  return (exp_hi * (exp_mid * l2));
102
}
103

104
} // namespace LIBC_NAMESPACE
105