libc/src/math/generic/expm1f.cpp

//===-- Implementation of expm1f function ---------------------------------===//
//
// 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
//
//===----------------------------------------------------------------------===//

#include "src/math/expm1f.h"
#include "src/__support/FPUtil/BasicOperations.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/common.h"
#include "src/math/expf.h"

namespace __llvm_libc {

// When |x| > Ln2, catastrophic cancellation does not occur with the
// subtraction expf(x) - 1.0f, so we use it to compute expm1f(x).
//
// We divide [-Ln2; Ln2] into 3 subintervals [-Ln2; -1/8], [-1/8; 1/8],
// [1/8; Ln2]. And we use a degree-6 polynomial to approximate exp(x) - 1 in
// each interval. The coefficients were generated by Sollya's fpminmax.
//
// See libc/utils/mathtools/expm1f.sollya for more detail.
LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
  const float Ln2 =
      0.69314718055994530941723212145817656807550013436025f; // For C++17:
                                                             // 0x1.62e'42ffp-1
  float abs_x = __llvm_libc::fputil::abs(x);

  if (abs_x <= Ln2) {
    if (abs_x <= 0.125f) {
      return x * __llvm_libc::fputil::polyeval(
                     x, 1.0f, 0.5f, 0.16666664183139801025390625f,
                     4.1666664183139801025390625e-2f,
                     8.3379410207271575927734375e-3f,
                     1.3894210569560527801513671875e-3f);
    }
    if (x > 0.125f) {
      return __llvm_libc::fputil::polyeval(
          x, 1.23142086749794543720781803131103515625e-7f,
          0.9999969005584716796875f, 0.500031292438507080078125f,
          0.16650259494781494140625f, 4.21491153538227081298828125e-2f,
          7.53940828144550323486328125e-3f,
          2.05591344274580478668212890625e-3f);
    }
    return __llvm_libc::fputil::polyeval(
        x, -6.899231408397099585272371768951416015625e-8f,
        0.999998271465301513671875f, 0.4999825656414031982421875f,
        0.16657467186450958251953125f, 4.1390590369701385498046875e-2f,
        7.856394164264202117919921875e-3f,
        9.380675037391483783721923828125e-4f);
  }
  return expf(x) - 1.0f;
}

} // namespace __llvm_libc