[libc] Add implementation of expm1f.

[lntue]

Summary

Use expm1f(x) = exp(x) - 1 for |x| > ln(2).
For |x| <= ln(2), divide it into 3 subintervals: [-ln2, -1/8], [-1/8, 1/8], [1/8, ln2]
and use a degree-6 polynomial approximation generated by Sollya's fpminmax for each interval.
Errors < 1.5 ULPs when we use fma to evaluate the polynomials.

[sivachandra]

LGTM but I have added a comment to add a sollya script. Good to go after that.

libc/src/math/generic/expm1f.cpp

28

I think we should add a sollya script for this under a directory named utils/mathtools so that future tuning can have a starting point to play with. And then, instead of inline comments like this, you can point to that sollya script at the top of this function.

libc/test/src/math/expm1f_test.cpp

31

Nitty nit: The name should probably be LlvmLibcExpm1fTest (ie. with the Exmpm1 capitalized.)

[lntue]

Add more detail about generating polynomial coefficients with Sollya to utils/mathtools/expm1f.sollya.
Add differential and perf tests.

[lntue]

Fix unit test name.

[sivachandra]

fma and fmaf have been enabled on the bots so this is ready to go.

libc/utils/mathtools/expm1f.sollya

1 ↗

(On Diff #341417)

Add LLVM license header.

9 ↗

(On Diff #341417)

Do we really need the > in the script?

38 ↗

(On Diff #341417)

Can we use print commands to make this nicer?

[lntue]

[libc] Add implementation of expm1f.

libc/config/linux/aarch64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.cosf
     libc.src.math.expf
     libc.src.math.exp2f
+    libc.src.math.expm1f
     libc.src.math.fabs
     libc.src.math.fabsf
     libc.src.math.fabsl
     libc.src.math.fdim
     libc.src.math.fdimf
     libc.src.math.fdiml
     libc.src.math.floor
     libc.src.math.floorf

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.copysignf
     libc.src.math.copysignl
     libc.src.math.ceil
     libc.src.math.ceilf
     libc.src.math.ceill
     libc.src.math.cosf
     libc.src.math.expf
     libc.src.math.exp2f
+    libc.src.math.expm1f
     libc.src.math.fabs
     libc.src.math.fabsf
     libc.src.math.fabsl
     libc.src.math.fdim
     libc.src.math.fdimf
     libc.src.math.fdiml
     libc.src.math.floor
     libc.src.math.floorf

libc/spec/stdc.td

HeaderSpec Math = HeaderSpec<
           FunctionSpec<"modff", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatPtr>]>,
           FunctionSpec<"modfl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoublePtr>]>,
 
           FunctionSpec<"cosf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"sinf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
 
           FunctionSpec<"expf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"exp2f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+          FunctionSpec<"expm1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
 
           FunctionSpec<"remainderf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
           FunctionSpec<"remainder", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
           FunctionSpec<"remainderl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
 
           FunctionSpec<"remquof", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>, ArgSpec<IntPtr>]>,
           FunctionSpec<"remquo", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>, ArgSpec<IntPtr>]>,
           FunctionSpec<"remquol", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>, ArgSpec<IntPtr>]>,

libc/src/math/CMakeLists.txt

 add_math_entrypoint_object(copysignl)
 
 add_math_entrypoint_object(cosf)
 
 add_math_entrypoint_object(expf)
 
 add_math_entrypoint_object(exp2f)
 
+add_math_entrypoint_object(expm1f)
+
 add_math_entrypoint_object(fabs)
 add_math_entrypoint_object(fabsf)
 add_math_entrypoint_object(fabsl)
 
 add_math_entrypoint_object(fdim)
 add_math_entrypoint_object(fdimf)
 add_math_entrypoint_object(fdiml)
 

libc/src/math/expm1f.h

+//===-- Implementation header for expm1f ------------------------*- C++ -*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_EXPM1F_H
+#define LLVM_LIBC_SRC_MATH_EXPM1F_H
+
+namespace __llvm_libc {
+
+float expm1f(float x);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_EXPM1F_H

libc/src/math/generic/CMakeLists.txt

HDRS
     ../exp2f.h
   DEPENDS
     .exp_utils
     .math_utils
     libc.include.math
 )
 
 add_entrypoint_object(
+  expm1f
+  SRCS
+    expm1f.cpp
+  HDRS
+    ../expm1f.h
+  DEPENDS
+    libc.include.math
+    libc.src.math.expf
+    libc.src.math.fabsf
+)
+
+add_entrypoint_object(
   copysign
   SRCS
     copysign.cpp
   HDRS
     ../copysign.h
   DEPENDS
     libc.utils.FPUtil.fputil
   COMPILE_OPTIONS

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/common.h"
+#include "src/math/expf.h"
+#include "utils/FPUtil/BasicOperations.h"
+#include "utils/FPUtil/PolyEval.h"
+
+namespace __llvm_libc {
+
+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) {
+    // We divide [-Ln2; Ln2] into 3 subintervals [-Ln2; -1/8], [-1/8; 1/8],
+    // [1/8; Ln2]. And we use a degree-5 polynomial to approximate exp(x) - 1 in
+    // each interval. The coefficients were generated by Sollya's fpminmax.
+    if (abs_x < 0.25f) {
+      // From Sollya command:
+      // > fpminimax(f, [|1, ..., 6|], [|1, 24...|], [-1/4, 1/4], relative);
+      // x * (1.0f + x * (0.5f + x * (0.166666328907012939453125f + x *
+      // (4.1666679084300994873046875e-2f + x * (8.350677788257598876953125e-3f
+      // + x * 1.39005668461322784423828125e-3f)))));
+      return x * __llvm_libc::fputil::polyeval(
+                     x, 1.0f, 0.5f, 0.166666328907012939453125f,
+                     4.1666679084300994873046875e-2f,
+                     8.350677788257598876953125e-3f,
+                     1.39005668461322784423828125e-3f);
+    }
+    if (x >= 0.25f) {
+      // From Sollya command:
+      // > fpminimax(f, [|0, ..., 5|], [|24...|], [1/8, log(2)], relative);
+      // 8.6965019363560713827610015869140625e-7f + x *
+      // (0.999985158443450927734375f + x * (0.500105917453765869140625f + x *
+      // (0.16625757515430450439453125f + x * (4.25875484943389892578125e-2f + x
+      // * (7.1335560642182826995849609375e-3f + x
+      // * 2.208019606769084930419921875e-3f)))));
+      return __llvm_libc::fputil::polyeval(
+          x, 8.6965019363560713827610015869140625e-7f,
+          0.999985158443450927734375f, 0.500105917453765869140625f,
+          0.16625757515430450439453125f, 4.25875484943389892578125e-2f,
+          7.1335560642182826995849609375e-3f,
+          2.208019606769084930419921875e-3f);
+    }
+    // From Sollya command:
+    // > fpminimax(f, [|0, ..., 5|], [|24...|], [-log(2), -1/4], relative);
+    // -4.05176479034707881510257720947265625e-7f + x *
+    // (0.999993026256561279296875f + x * (0.4999495446681976318359375f + x *
+    // (0.166467368602752685546875f + x * (4.12006340920925140380859375e-2f + x
+    // * (7.68242776393890380859375e-3f + x
+    // * 8.735431474633514881134033203125e-4f)))));
+    return __llvm_libc::fputil::polyeval(
+        x, -4.05176479034707881510257720947265625e-7f,
+        0.999993026256561279296875f, 0.4999495446681976318359375f,
+        0.166467368602752685546875f, 4.12006340920925140380859375e-2f,
+        7.68242776393890380859375e-3f, 8.735431474633514881134033203125e-4f);
+  }
+  // When |x| > Ln2, catastrophic cancellation does not occur with the
+  // subtruction expf(x) - 1.0f.
+  return expf(x) - 1.0f;
+}
+
+} // namespace __llvm_libc

libc/test/src/math/CMakeLists.txt

add_fp_unittest(
   SRCS
     fma_test.cpp
   DEPENDS
     libc.include.math
     libc.src.math.fma
     libc.utils.FPUtil.fputil
 )
 
+add_fp_unittest(
+  expm1f_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    expm1f_test.cpp
+  DEPENDS
+    libc.include.errno
+    libc.include.math
+    libc.src.math.expm1f
+    libc.utils.FPUtil.fputil
+)
+
 add_subdirectory(generic)
 add_subdirectory(exhaustive)

libc/test/src/math/exhaustive/CMakeLists.txt

 add_libc_exhaustive_testsuite(libc_math_exhaustive_tests)
 
 add_fp_unittest(
   sqrtf_test
   NEED_MPFR
   SUITE
     libc_math_exhaustive_tests
   SRCS
     sqrtf_test.cpp
   DEPENDS
     libc.include.math
     libc.src.math.sqrtf
     libc.utils.FPUtil.fputil
 )
+
+
+add_fp_unittest(
+  expm1f_test
+  NEED_MPFR
+  SUITE
+    libc_math_exhaustive_tests
+  SRCS
+    expm1f_test.cpp
+  DEPENDS
+    libc.include.math
+    libc.src.math.expm1f
+    libc.utils.FPUtil.fputil
+)
 No newline at end of file

libc/test/src/math/exhaustive/expm1f_test.cpp

+//===-- Exhaustive test for expm1f
+//-----------------------------------------===//
+//
+// 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 "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include <math.h>
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+TEST(LlvmLibcExpm1fExhaustiveTest, AllValues) {
+  uint32_t bits = 0;
+  do {
+    FPBits x(bits);
+    if (!x.isInfOrNaN() && float(x) < 88.70f) {
+      ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, float(x),
+                        __llvm_libc::expm1f(float(x)), 1.5);
+    }
+  } while (bits++ < 0xffff'ffffU);
+}

libc/test/src/math/expm1f_test.cpp

+//===-- Unittests for expm1f
+//------------------------------------------------===//
+//
+// 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 "utils/FPUtil/BitPatterns.h"
+#include "utils/FPUtil/ClassificationFunctions.h"
+#include "utils/FPUtil/FloatOperations.h"
+#include "utils/FPUtil/FloatProperties.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+using __llvm_libc::fputil::isNegativeQuietNaN;
+using __llvm_libc::fputil::isQuietNaN;
+using __llvm_libc::fputil::valueAsBits;
+using __llvm_libc::fputil::valueFromBits;
+
+using BitPatterns = __llvm_libc::fputil::BitPatterns<float>;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+TEST(LlvmLibcexpm1fTest, SpecialNumbers) {
+  errno = 0;
+
+  EXPECT_TRUE(
+      isQuietNaN(__llvm_libc::expm1f(valueFromBits(BitPatterns::aQuietNaN))));
+  EXPECT_EQ(errno, 0);
+
+  EXPECT_TRUE(isNegativeQuietNaN(
+      __llvm_libc::expm1f(valueFromBits(BitPatterns::aNegativeQuietNaN))));
+  EXPECT_EQ(errno, 0);
+
+  EXPECT_TRUE(isQuietNaN(
+      __llvm_libc::expm1f(valueFromBits(BitPatterns::aSignallingNaN))));
+  EXPECT_EQ(errno, 0);
+
+  EXPECT_TRUE(isNegativeQuietNaN(
+      __llvm_libc::expm1f(valueFromBits(BitPatterns::aNegativeSignallingNaN))));
+  EXPECT_EQ(errno, 0);
+
+  EXPECT_EQ(BitPatterns::inf,
+            valueAsBits(__llvm_libc::expm1f(valueFromBits(BitPatterns::inf))));
+  EXPECT_EQ(errno, 0);
+
+  EXPECT_EQ(BitPatterns::negOne, valueAsBits(__llvm_libc::expm1f(
+                                     valueFromBits(BitPatterns::negInf))));
+  EXPECT_EQ(errno, 0);
+
+  EXPECT_EQ(BitPatterns::zero,
+            valueAsBits(__llvm_libc::expm1f(valueFromBits(BitPatterns::zero))));
+  EXPECT_EQ(errno, 0);
+
+  EXPECT_EQ(BitPatterns::negZero, valueAsBits(__llvm_libc::expm1f(
+                                      valueFromBits(BitPatterns::negZero))));
+  EXPECT_EQ(errno, 0);
+}
+
+TEST(LlvmLibcexpm1fTest, Overflow) {
+  errno = 0;
+  EXPECT_EQ(BitPatterns::inf,
+            valueAsBits(__llvm_libc::expm1f(valueFromBits(0x7f7fffffU))));
+  EXPECT_EQ(errno, ERANGE);
+
+  errno = 0;
+  EXPECT_EQ(BitPatterns::inf,
+            valueAsBits(__llvm_libc::expm1f(valueFromBits(0x42cffff8U))));
+  EXPECT_EQ(errno, ERANGE);
+
+  errno = 0;
+  EXPECT_EQ(BitPatterns::inf,
+            valueAsBits(__llvm_libc::expm1f(valueFromBits(0x42d00008U))));
+  EXPECT_EQ(errno, ERANGE);
+}
+
+TEST(LlvmLibcexpm1fTest, Underflow) {
+  errno = 0;
+  EXPECT_EQ(BitPatterns::negOne,
+            valueAsBits(__llvm_libc::expm1f(valueFromBits(0xff7fffffU))));
+  EXPECT_EQ(errno, ERANGE);
+
+  errno = 0;
+  EXPECT_EQ(BitPatterns::negOne,
+            valueAsBits(__llvm_libc::expm1f(valueFromBits(0xc2cffff8U))));
+  EXPECT_EQ(errno, ERANGE);
+
+  errno = 0;
+  EXPECT_EQ(BitPatterns::negOne,
+            valueAsBits(__llvm_libc::expm1f(valueFromBits(0xc2d00008U))));
+  EXPECT_EQ(errno, ERANGE);
+}
+
+// Test with inputs which are the borders of underflow/overflow but still
+// produce valid results without setting errno.
+TEST(LlvmLibcexpm1fTest, Borderline) {
+  float x;
+
+  errno = 0;
+  x = valueFromBits(0x42affff8U);
+  ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.0);
+  EXPECT_EQ(errno, 0);
+
+  x = valueFromBits(0x42b00008U);
+  ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.0);
+  EXPECT_EQ(errno, 0);
+
+  x = valueFromBits(0xc2affff8U);
+  ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.0);
+  EXPECT_EQ(errno, 0);
+
+  x = valueFromBits(0xc2b00008U);
+  ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 1.0);
+  EXPECT_EQ(errno, 0);
+}
+
+TEST(LlvmLibcexpm1fTest, InFloatRange) {
+  constexpr uint32_t count = 1000000;
+  constexpr uint32_t step = UINT32_MAX / count;
+  for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) {
+    float x = valueFromBits(v);
+    if (isnan(x) || isinf(x))
+      continue;
+    errno = 0;
+    float result = __llvm_libc::expm1f(x);
+
+    // If the computation resulted in an error or did not produce valid result
+    // in the single-precision floating point range, then ignore comparing with
+    // MPFR result as MPFR can still produce valid results because of its
+    // wider precision.
+    if (isnan(result) || isinf(result) || errno != 0)
+      continue;
+    ASSERT_MPFR_MATCH(mpfr::Operation::Expm1, x, __llvm_libc::expm1f(x), 2.0);
+  }
+}

libc/utils/FPUtil/BitPatterns.h

template <> struct BitPatterns<float> {
 
   static constexpr BitsType inf = 0x7f800000U;
   static constexpr BitsType negInf = 0xff800000U;
 
   static constexpr BitsType zero = 0x0;
   static constexpr BitsType negZero = 0x80000000U;
 
   static constexpr BitsType one = 0x3f800000U;
+  static constexpr BitsType negOne = 0xbf800000U;
 
   // Examples of quiet NAN.
   static constexpr BitsType aQuietNaN = 0x7fc00000U;
   static constexpr BitsType aNegativeQuietNaN = 0xffc00000U;
 
   // Examples of signalling NAN.
   static constexpr BitsType aSignallingNaN = 0x7f800001U;
   static constexpr BitsType aNegativeSignallingNaN = 0xff800001U;

libc/utils/FPUtil/CMakeLists.txt

HDRS
     FEnv.h
     FloatOperations.h
     FloatProperties.h
     FPBits.h
     BasicOperations.h
     ManipulationFunctions.h
     NearestIntegerOperations.h
     NormalFloat.h
+    PolyEval.h
   DEPENDS
     libc.include.math
     libc.include.errno
     libc.include.fenv
     libc.src.__support.common
     libc.utils.CPP.standalone_cpp
 )
 

libc/utils/FPUtil/PolyEval.h

+//===-- Common header for PolyEval implementations --------------*- C++ -*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_UTILS_FPUTIL_POLYEVAL_H
+#define LLVM_LIBC_UTILS_FPUTIL_POLYEVAL_H
+
+#include "utils/CPP/TypeTraits.h"
+
+// Evaluate polynomial using Horner's Scheme:
+// With polyeval(x, a_0, a_1, ..., a_n) = a_n * x^n + ... + a_1 * x + a_0, we
+// evaluated it as:  a_0 + x * (a_1 + x * ( ... (a_(n-1) + x * a_n) ... ) ) ).
+// We will use fma instructions if available.
+// Example: to evaluate x^3 + 2*x^2 + 3*x + 4, call
+//   polyeval( x, 4.0, 3.0, 2.0, 1.0 )
+
+#if defined(__x86_64__) || defined(__aarch64__)
+#include "FMA.h"
+
+namespace __llvm_libc {
+namespace fputil {
+
+template <typename T> static inline T polyeval(T x, T a0) { return a0; }
+
+template <typename T, typename... Ts>
+static inline T polyeval(T x, T a0, Ts... a) {
+  return fma(x, polyeval(x, a...), a0);
+}
+
+} // namespace fputil
+} // namespace __llvm_libc
+
+#else
+
+namespace __llvm_libc {
+namespace fputil {
+
+template <typename T> static inline T polyeval(T x, T a0) { return a0; }
+
+template <typename T, typename... Ts>
+static inline T polyeval(T x, T a0, Ts... a) {
+  return x * polyeval(x, a...) + a0;
+}
+
+} // namespace fputil
+} // namespace __llvm_libc
+
+#endif
+
+#endif // LLVM_LIBC_UTILS_FPUTIL_FMA_H

libc/utils/FPUtil/generic/FMA.h

static inline cpp::EnableIfType<cpp::IsSame<T, float>::Value, T> fma(T x, T y,
 
   return static_cast<float>(static_cast<double>(bit_sum));
 }
 
 } // namespace generic
 } // namespace fputil
 } // namespace __llvm_libc
 
-#endif // Generic fma implementations
-
 #endif // LLVM_LIBC_UTILS_FPUTIL_GENERIC_FMA_H

libc/utils/MPFRWrapper/MPFRUtils.h

enum class Operation : int {
   // and produce a single floating point number as output. The input
   // and output floating point numbers are of the same kind.
   BeginUnaryOperationsSingleOutput,
   Abs,
   Ceil,
   Cos,
   Exp,
   Exp2,
+  Expm1,
   Floor,
   Round,
   Sin,
   Sqrt,
   Trunc,
   EndUnaryOperationsSingleOutput,
 
   // Operations which take a single floating point nubmer as input

libc/utils/MPFRWrapper/MPFRUtils.cpp

 #ifdef CUSTOM_MPFR_INCLUDER
 // Some downstream repos are monoliths carrying MPFR sources in their third
 // party directory. In such repos, including the MPFR header as
 // `#include <mpfr.h>` is either disallowed or not possible. If that is the
 // case, a file named `CustomMPFRIncluder.h` should be added through which the
 // MPFR header can be included in manner allowed in that repo.
 #include "CustomMPFRIncluder.h"
 #else
 #include <mpfr.h>
 #endif
 
 template <typename T> using FPBits = __llvm_libc::fputil::FPBits<T>;
 
 namespace __llvm_libc {
 namespace testing {
 namespace mpfr {
 
public:
   }
 
   MPFRNumber exp2() const {
     MPFRNumber result;
     mpfr_exp2(result.value, value, MPFR_RNDN);
     return result;
   }
 
+  MPFRNumber expm1() const {
+    MPFRNumber result;
+    mpfr_expm1(result.value, value, MPFR_RNDN);
+    return result;
+  }
+
   MPFRNumber floor() const {
     MPFRNumber result;
     mpfr_floor(result.value, value);
     return result;
   }
 
   MPFRNumber frexp(int &exp) {
     MPFRNumber result;
unaryOperation(Operation op, InputType input) {
   case Operation::Ceil:
     return mpfrInput.ceil();
   case Operation::Cos:
     return mpfrInput.cos();
   case Operation::Exp:
     return mpfrInput.exp();
   case Operation::Exp2:
     return mpfrInput.exp2();
+  case Operation::Expm1:
+    return mpfrInput.expm1();
   case Operation::Floor:
     return mpfrInput.floor();
   case Operation::Round:
     return mpfrInput.round();
   case Operation::Sin:
     return mpfrInput.sin();
   case Operation::Sqrt:
     return mpfrInput.sqrt();