libc/src/stdfix/exphk.cpp

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