1 //===-- Implementation of exphk 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 "exphk.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 89 do { 21 // hi = floor(i/8) - 5; 22 // m = i/8 - floor(i/8) - 0.5; 23 // e_hi = nearestint(exp(hi) * 2^7) * 2^-7; 24 // e_mid = nearestint(exp(m) * 2^7) * 2^-7; 25 // print(hi, e_hi, m, e_mid); 26 // }; 27 // Notice that when i = 88 and 89, e_hi will overflow short accum range. 28 static constexpr short accum EXP_HI[12] = { 29 0x1.0p-7hk, 0x1.0p-6hk, 0x1.8p-5hk, 0x1.1p-3hk, 0x1.78p-2hk, 0x1.0p0hk, 30 0x1.5cp1hk, 0x1.d9p2hk, 0x1.416p4hk, 0x1.b4dp5hk, 0x1.28d4p7hk, SACCUM_MAX, 31 }; 32 33 static constexpr short accum EXP_MID[8] = { 34 0x1.38p-1hk, 0x1.6p-1hk, 0x1.9p-1hk, 0x1.c4p-1hk, 35 0x1.0p0hk, 0x1.22p0hk, 0x1.48p0hk, 0x1.74p0hk, 36 }; 37 38 } // anonymous namespace 39 40 LLVM_LIBC_FUNCTION(short accum, exphk, (short accum x)) { 41 using FXRep = fixed_point::FXRep<short accum>; 42 using StorageType = typename FXRep::StorageType; 43 // Output overflow 44 if (LIBC_UNLIKELY(x >= 0x1.64p2hk)) 45 return FXRep::MAX(); 46 // Lower bound where exp(x) -> 0: 47 // floor(log(2^-8) * 2^7) * 2^-7 48 if (LIBC_UNLIKELY(x <= -0x1.63p2hk)) 49 return FXRep::ZERO(); 50 51 // Current range of x: 52 // -0x1.628p2 <= x <= 0x1.638p2 53 // Range reduction: 54 // x = hi + mid + lo, 55 // where: 56 // hi is an integer 57 // mid * 2^3 is an integer 58 // |lo| <= 2^-4. 59 // Then exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo) 60 // ~ exp(hi) * exp(mid) * (1 + lo) 61 // with relative errors < |lo|^2 <= 2^-8. 62 // exp(hi) and exp(mid) are extracted from small lookup tables. 63 64 // Round-to-nearest 1/8, tie-to-(+Int): 65 constexpr short accum ONE_SIXTEENTH = 0x1.0p-4hk; 66 // x_rounded = floor(x + 1/16). 67 short accum x_rounded = ((x + ONE_SIXTEENTH) >> (FXRep::FRACTION_LEN - 3)) 68 << (FXRep::FRACTION_LEN - 3); 69 short accum lo = x - x_rounded; 70 71 // Range of x_rounded: 72 // x_rounded >= floor((-0x1.628p2 + 0x1.0p-4) * 2^3) * 2^-3 73 // = -0x1.6p2 = -5.5 74 // To get the indices, we shift the values so that it start with 0. 75 // Range of indices: 0 <= indices <= 89 76 StorageType indices = cpp::bit_cast<StorageType>((x_rounded + 0x1.6p2hk) >> 77 (FXRep::FRACTION_LEN - 3)); 78 // So we have the following relation: 79 // indices = (hi + mid + 44/8) * 8 80 // That implies: 81 // hi + mid = indices/8 - 5.5 82 // So for lookup tables, we can use the upper 4 bits to get: 83 // exp( floor(indices / 8) - 5 ) 84 // and lower 3 bits for: 85 // exp( (indices - floor(indices)) - 0.5 ) 86 short accum exp_hi = EXP_HI[indices >> 3]; 87 short accum exp_mid = EXP_MID[indices & 0x7]; 88 // exp(x) ~ exp(hi) * exp(mid) * (1 + lo); 89 return (exp_hi * (exp_mid * (0x1.0p0hk + lo))); 90 } 91 92 } // namespace LIBC_NAMESPACE 93