diff --git a/libc/src/math/generic/expxf.h b/libc/src/math/generic/expxf.h
new file mode 100644
--- /dev/null
+++ b/libc/src/math/generic/expxf.h
@@ -0,0 +1,75 @@
+//===-- Single-precision general exp 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H
+#define LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H
+
+#include "common_constants.h" // Lookup tables EXP_M
+#include "math_utils.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/common.h"
+#include <src/__support/FPUtil/NearestIntegerOperations.h>
+
+#include <errno.h>
+
+namespace __llvm_libc {
+
+// The algorithm represents exp(x) as
+//   exp(x) = 2^(ln(2) * i) * 2^(ln(2) * j / NUM_P )) * exp(dx)
+// where i integer value, j integer in range [-NUM_P/2, NUM_P/2).
+// 2^(ln(2) * j / NUM_P )) is a table values: 1.0 + EXP_M
+// exp(dx) calculates by taylor expansion.
+
+// Inversion of ln(2). Multiplication by EXP_num_p due to sampling by 1 /
+// EXP_num_p Precise value of the constant is not needed.
+static constexpr double LN2_INV = 1.44269504088896340735 * EXP_num_p;
+
+// LN2_HIGH + LN2_LOW = ln(2) with precision higher than double(ln(2))
+static constexpr double LN2_HIGH = 0x1.62e42fefa0000p-1 / EXP_num_p;
+static constexpr double LN2_LOW = 0x1.cf79abc9e3b3ap-40 / EXP_num_p;
+
+struct exe_eval_result_t {
+  // exp(x) = mult_exp * (r + 1.0)
+  // where mult_exp = 2^e and r in range [~-0.3, ~0.41]
+  double mult_exp;
+  double r;
+};
+
+// Value to add to make integer digits always positive.
+// Needed for correct rounding to the nearest value in both positive and
+// negative range.
+static exe_eval_result_t exp_eval(double x) {
+  int64_t ps = fputil::multiply_add(LN2_INV, x, double(ADD_TO_POS) + 0.5);
+  // Negative sign due to multiply_add optimization
+  double mult_e1, ml;
+  {
+    uint64_t ps_shifted = ps + (1 << (EXP_bits_p - 1));
+    int64_t dg = int64_t(ps_shifted >> EXP_bits_p) - (ADD_TO_POS >> EXP_bits_p);
+    int table_index = ps_shifted & (EXP_num_p - 1);
+    fputil::FPBits<double> bs;
+    bs.set_unbiased_exponent(fputil::FPBits<double>::EXPONENT_BIAS + dg);
+    ml = EXP_2_POW[table_index];
+    mult_e1 = bs.get_val();
+  }
+  double psx = double(ADD_TO_POS) - ps;
+  double dx = fputil::multiply_add(psx, LN2_LOW,
+                                   fputil::multiply_add(psx, LN2_HIGH, x));
+
+  // Taylor series coefficients
+  double pe = dx * fputil::polyeval(dx, 1.0, 1.0 / 2.0, 1.0 / 6.0, 1.0 / 24.0,
+                                    1.0 / 120.0, 1.0 / 720.0);
+
+  double r = fputil::multiply_add(ml, pe, pe) + ml;
+  return {mult_e1, r};
+}
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_GENERIC_EXPXF_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1299,6 +1299,18 @@
     libc.src.__support.FPUtil.fputil
 )
 
+add_fp_unittest(
+  expxf_test
+    NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    expxf_test.cpp
+  DEPENDS
+    libc.src.math.generic.common_constants
+    libc.src.__support.FPUtil.fputil
+)
+
 add_subdirectory(generic)
 add_subdirectory(exhaustive)
 add_subdirectory(differential_testing)
diff --git a/libc/test/src/math/expxf_test.cpp b/libc/test/src/math/expxf_test.cpp
new file mode 100644
--- /dev/null
+++ b/libc/test/src/math/expxf_test.cpp
@@ -0,0 +1,37 @@
+//===-- Unittests for expf ------------------------------------------------===//
+//
+// 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/__support/FPUtil/FPBits.h"
+#include "src/math/generic/expxf.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/FPMatcher.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+#include <errno.h>
+#include <stdint.h>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+DECLARE_SPECIAL_CONSTANTS(float)
+
+TEST(LlvmLibcExpfTest, InFloatRange) {
+  constexpr uint32_t COUNT = 1000000;
+  constexpr uint32_t STEP = UINT32_MAX / COUNT;
+  auto fx = [](float x) -> float {
+    auto result = __llvm_libc::exp_eval(x);
+    return static_cast<float>(result.mult_exp * result.r + result.mult_exp);
+  };
+  for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+    float x = float(FPBits(v));
+    if (isnan(x) || isinf(x) || x < -70 || x > 70 || fabsf(x) < 0x1.0p-10)
+      continue;
+
+    EXPECT_MPFR_MATCH_DEFAULT(mpfr::Operation::Exp, x, fx(x), 0.5);
+  }
+}