diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -65,6 +65,7 @@
     libc.src.math.floor
     libc.src.math.floorf
     libc.src.math.floorl
+    libc.src.math.fmaf
     libc.src.math.fmax
     libc.src.math.fmaxf
     libc.src.math.fmaxl
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -106,6 +106,7 @@
     libc.src.math.floor
     libc.src.math.floorf
     libc.src.math.floorl
+    libc.src.math.fmaf
     libc.src.math.fmin
     libc.src.math.fminf
     libc.src.math.fminl
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -322,6 +322,8 @@
           FunctionSpec<"fmaxf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
           FunctionSpec<"fmaxl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
 
+          FunctionSpec<"fmaf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>, ArgSpec<FloatType>]>,
+
           FunctionSpec<"frexp", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<IntPtr>]>,
           FunctionSpec<"frexpf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<IntPtr>]>,
           FunctionSpec<"frexpl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<IntPtr>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -905,3 +905,15 @@
   COMPILE_OPTIONS
     -O2
 )
+
+add_entrypoint_object(
+  fmaf
+  SRCS
+    fmaf.cpp
+  HDRS
+    fmaf.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
diff --git a/libc/src/math/fmaf.h b/libc/src/math/fmaf.h
new file mode 100644
--- /dev/null
+++ b/libc/src/math/fmaf.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for fmaf --------------------------*- 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_FMAF_H
+#define LLVM_LIBC_SRC_MATH_FMAF_H
+
+namespace __llvm_libc {
+
+float fmaf(float x, float y, float z);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_FMAF_H
diff --git a/libc/src/math/fmaf.cpp b/libc/src/math/fmaf.cpp
new file mode 100644
--- /dev/null
+++ b/libc/src/math/fmaf.cpp
@@ -0,0 +1,47 @@
+//===-- Implementation of fmaf 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/__support/common.h"
+
+#include "utils/FPUtil/FEnv.h"
+#include "utils/FPUtil/FPBits.h"
+
+namespace __llvm_libc {
+
+float LLVM_LIBC_ENTRYPOINT(fmaf)(float x, float y, float z) {
+  fputil::FPBits<float> bitx(x), bity(y), bitz(z);
+
+  // Product is exact.
+  fputil::FPBits<double> bit_prod(static_cast<double>(x) *
+                                  static_cast<double>(y));
+  fputil::FPBits<double> bitz_d(static_cast<double>(z));
+  fputil::FPBits<double> bit_sum(static_cast<double>(bit_prod) +
+                                 static_cast<double>(bitz_d));
+
+  if (!(bit_sum.isInfOrNaN() || bit_sum.isZero())) {
+    // Dekker's 2Sum: find t such that s - t = prod + z exactly
+    fputil::FPBits<double> t(
+        (bit_prod.exponent >= bitz_d.exponent)
+            ? ((static_cast<double>(bit_sum) - bit_prod) - bitz)
+            : ((static_cast<double>(bit_sum) - bitz) - bit_prod));
+
+    // Update sticky bits if t != 0.0 and the least (52 - 23 - 1 = 28) bits are
+    // zero.
+    if (!t.isZero() && ((bit_sum.mantissa & 0xfff'ffffULL) == 0)) {
+      if (bit_sum.sign != t.sign) {
+        ++bit_sum.mantissa;
+      } else if (bit_sum.mantissa) {
+        --bit_sum.mantissa;
+      }
+    }
+  }
+
+  return static_cast<float>(static_cast<double>(bit_sum));
+}
+
+} // namespace __llvm_libc
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
@@ -1007,3 +1007,16 @@
     libc.src.math.hypot
     libc.utils.FPUtil.fputil
 )
+
+add_fp_unittest(
+  fmaf_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    fmaf_test.cpp
+  DEPENDS
+    libc.include.math
+    libc.src.math.fmaf
+    libc.utils.FPUtil.fputil
+)
diff --git a/libc/test/src/math/FmaTest.h b/libc/test/src/math/FmaTest.h
new file mode 100644
--- /dev/null
+++ b/libc/test/src/math/FmaTest.h
@@ -0,0 +1,89 @@
+//===-- Utility class to test different flavors of fma  -------------------===//
+//
+// 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_TEST_SRC_MATH_FMATEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_FMATEST_H
+
+#include "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+#include <random>
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+template <typename T>
+class FmaTestTemplate : public __llvm_libc::testing::Test {
+private:
+  using Func = T (*)(T, T, T);
+  using FPBits = __llvm_libc::fputil::FPBits<T>;
+  using UIntType = typename FPBits::UIntType;
+  const T nan = __llvm_libc::fputil::FPBits<T>::buildNaN(1);
+  const T inf = __llvm_libc::fputil::FPBits<T>::inf();
+  const T negInf = __llvm_libc::fputil::FPBits<T>::negInf();
+  const T zero = __llvm_libc::fputil::FPBits<T>::zero();
+  const T negZero = __llvm_libc::fputil::FPBits<T>::negZero();
+
+  UIntType getRandomBitPattern() {
+    UIntType bits{0};
+    for (size_t i = 0; i < sizeof(UIntType) / 2; ++i) {
+      bits = (bits << 2) + static_cast<uint16_t>(std::rand());
+    }
+    return bits;
+  }
+
+public:
+  void testSpecialNumbers(Func func) {
+    EXPECT_FP_EQ(func(zero, zero, zero), zero);
+    EXPECT_FP_EQ(func(zero, negZero, negZero), negZero);
+    EXPECT_FP_EQ(func(inf, inf, zero), inf);
+    EXPECT_FP_EQ(func(negInf, inf, negInf), negInf);
+    EXPECT_FP_EQ(func(inf, zero, zero), nan);
+    EXPECT_FP_EQ(func(inf, negInf, inf), nan);
+    EXPECT_FP_EQ(func(nan, zero, inf), nan);
+    EXPECT_FP_EQ(func(inf, negInf, nan), nan);
+
+    EXPECT_FP_EQ(func(T(0.5), FPBits(FPBits::minSubnormal),
+                      FPBits(FPBits::minSubnormal)),
+                 FPBits(UIntType(2)));
+    FPBits v(FPBits::minNormal + UIntType(1));
+    EXPECT_FP_EQ(
+        func(T(1) / T(FPBits::minNormal << 1), v, FPBits(FPBits::minNormal)),
+        v);
+  }
+
+  void testSubnormalRange(Func func) {
+    constexpr UIntType count = 1000001;
+    constexpr UIntType step =
+        (FPBits::maxSubnormal - FPBits::minSubnormal) / count;
+    for (UIntType v = FPBits::minSubnormal, w = FPBits::maxSubnormal;
+         v <= FPBits::maxSubnormal && w >= FPBits::minSubnormal;
+         v += step, w -= step) {
+      T x = FPBits(getRandomBitPattern()), y = FPBits(v), z = FPBits(w);
+      T result = func(x, y, z);
+      mpfr::TernaryInput<T> input{x, y, z};
+      ASSERT_MPFR_MATCH(mpfr::Operation::Fma, input, result, 0.5);
+    }
+  }
+
+  void testNormalRange(Func func) {
+    constexpr UIntType count = 1000001;
+    constexpr UIntType step = (FPBits::maxNormal - FPBits::minNormal) / count;
+    for (UIntType v = FPBits::minNormal, w = FPBits::maxNormal;
+         v <= FPBits::maxNormal && w >= FPBits::minNormal;
+         v += step, w -= step) {
+      T x = FPBits(v), y = FPBits(w), z = FPBits(getRandomBitPattern());
+      T result = func(x, y, z);
+      mpfr::TernaryInput<T> input{x, y, z};
+      ASSERT_MPFR_MATCH(mpfr::Operation::Fma, input, result, 0.5);
+    }
+  }
+};
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_FMATEST_H
diff --git a/libc/test/src/math/fmaf_test.cpp b/libc/test/src/math/fmaf_test.cpp
new file mode 100644
--- /dev/null
+++ b/libc/test/src/math/fmaf_test.cpp
@@ -0,0 +1,19 @@
+//===-- Unittests for fmaf  -----------------------------------------------===//
+//
+// 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 "FmaTest.h"
+
+#include "src/math/fmaf.h"
+
+using FmaTest = FmaTestTemplate<float>;
+
+TEST_F(FmaTest, SpecialNumbers) { testSpecialNumbers(&__llvm_libc::fmaf); }
+
+TEST_F(FmaTest, SubnormalRange) { testSubnormalRange(&__llvm_libc::fmaf); }
+
+TEST_F(FmaTest, NormalRange) { testNormalRange(&__llvm_libc::fmaf); }
diff --git a/libc/utils/FPUtil/FPBits.h b/libc/utils/FPUtil/FPBits.h
--- a/libc/utils/FPUtil/FPBits.h
+++ b/libc/utils/FPUtil/FPBits.h
@@ -84,7 +84,10 @@
   // We don't want accidental type promotions/conversions so we require exact
   // type match.
   template <typename XType,
-            cpp::EnableIfType<cpp::IsSame<T, XType>::Value, int> = 0>
+            cpp::EnableIfType<cpp::IsSame<T, XType>::Value ||
+                                  (cpp::IsIntegral<XType>::Value &&
+                                   (sizeof(XType) == sizeof(UIntType))),
+                              int> = 0>
   explicit FPBits(XType x) {
     *this = *reinterpret_cast<FPBits<T> *>(&x);
   }
@@ -106,13 +109,6 @@
   // the potential software implementations of UIntType will not slow real
   // code.
 
-  template <typename XType,
-            cpp::EnableIfType<cpp::IsSame<UIntType, XType>::Value, int> = 0>
-  explicit FPBits<long double>(XType x) {
-    // The last 4 bytes of v are ignored in case of i386.
-    *this = *reinterpret_cast<FPBits<T> *>(&x);
-  }
-
   UIntType bitsAsUInt() const {
     return *reinterpret_cast<const UIntType *>(this);
   }
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -57,8 +57,11 @@
   RemQuo, // The first output, the floating point output, is the remainder.
   EndBinaryOperationsTwoOutputs,
 
+  // Operations wich take three floating point nubmers of the same type as
+  // input and produce a single floating point number of the same type as
+  // output.
   BeginTernaryOperationsSingleOuput,
-  // TODO: Add operations like fma.
+  Fma,
   EndTernaryOperationsSingleOutput,
 };
 
@@ -113,6 +116,11 @@
 bool compareBinaryOperationOneOutput(Operation op, const BinaryInput<T> &input,
                                      T libcOutput, double t);
 
+template <typename T>
+bool compareTernaryOperationOneOutput(Operation op,
+                                      const TernaryInput<T> &input,
+                                      T libcOutput, double t);
+
 template <typename T>
 void explainUnaryOperationSingleOutputError(Operation op, T input, T matchValue,
                                             testutils::StreamWrapper &OS);
@@ -132,6 +140,12 @@
                                           T matchValue,
                                           testutils::StreamWrapper &OS);
 
+template <typename T>
+void explainTernaryOperationOneOutputError(Operation op,
+                                           const TernaryInput<T> &input,
+                                           T matchValue,
+                                           testutils::StreamWrapper &OS);
+
 template <Operation op, typename InputType, typename OutputType>
 class MPFRMatcher : public testing::Matcher<OutputType> {
   InputType input;
@@ -174,7 +188,7 @@
 
   template <typename T>
   static bool match(const TernaryInput<T> &in, T out, double tolerance) {
-    // TODO: Implement the comparision function and error reporter.
+    return compareTernaryOperationOneOutput(op, in, out, tolerance);
   }
 
   template <typename T>
@@ -199,6 +213,12 @@
                            testutils::StreamWrapper &OS) {
     explainBinaryOperationOneOutputError(op, in, out, OS);
   }
+
+  template <typename T>
+  static void explainError(const TernaryInput<T> &in, T out,
+                           testutils::StreamWrapper &OS) {
+    explainTernaryOperationOneOutputError(op, in, out, OS);
+  }
 };
 
 } // namespace internal
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -35,48 +35,71 @@
 namespace testing {
 namespace mpfr {
 
+template <typename T> struct Precision {
+  static constexpr unsigned int value = 128;
+};
+
+template <> struct Precision<float> {
+  static constexpr unsigned int value = 24;
+};
+
+template <> struct Precision<double> {
+  static constexpr unsigned int value = 53;
+};
+
+#if !(defined(__x86_64__) || defined(__i386__))
+template <> struct Precision<long double> {
+  static constexpr unsigned int value = 64;
+};
+#else
+template <> struct Precision<long double> {
+  static constexpr unsigned int value = 113;
+};
+#endif
+
 class MPFRNumber {
   // A precision value which allows sufficiently large additional
   // precision even compared to quad-precision floating point values.
-  static constexpr unsigned int mpfrPrecision = 128;
+  unsigned int mpfrPrecision;
 
   mpfr_t value;
 
 public:
-  MPFRNumber() { mpfr_init2(value, mpfrPrecision); }
+  MPFRNumber() : mpfrPrecision(128) { mpfr_init2(value, mpfrPrecision); }
 
   // We use explicit EnableIf specializations to disallow implicit
   // conversions. Implicit conversions can potentially lead to loss of
   // precision.
   template <typename XType,
             cpp::EnableIfType<cpp::IsSame<float, XType>::Value, int> = 0>
-  explicit MPFRNumber(XType x) {
+  explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
     mpfr_init2(value, mpfrPrecision);
     mpfr_set_flt(value, x, MPFR_RNDN);
   }
 
   template <typename XType,
             cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0>
-  explicit MPFRNumber(XType x) {
+  explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
     mpfr_init2(value, mpfrPrecision);
     mpfr_set_d(value, x, MPFR_RNDN);
   }
 
   template <typename XType,
             cpp::EnableIfType<cpp::IsSame<long double, XType>::Value, int> = 0>
-  explicit MPFRNumber(XType x) {
+  explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
     mpfr_init2(value, mpfrPrecision);
     mpfr_set_ld(value, x, MPFR_RNDN);
   }
 
   template <typename XType,
             cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0>
-  explicit MPFRNumber(XType x) {
+  explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
     mpfr_init2(value, mpfrPrecision);
     mpfr_set_sj(value, x, MPFR_RNDN);
   }
 
-  MPFRNumber(const MPFRNumber &other) {
+  MPFRNumber(const MPFRNumber &other) : mpfrPrecision(other.mpfrPrecision) {
+    mpfr_init2(value, mpfrPrecision);
     mpfr_set(value, other.value, MPFR_RNDN);
   }
 
@@ -85,6 +108,7 @@
   }
 
   MPFRNumber &operator=(const MPFRNumber &rhs) {
+    mpfrPrecision = rhs.mpfrPrecision;
     mpfr_set(value, rhs.value, MPFR_RNDN);
     return *this;
   }
@@ -193,6 +217,12 @@
     return result;
   }
 
+  MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
+    MPFRNumber result(*this);
+    mpfr_fma(result.value, value, b.value, c.value, MPFR_RNDN);
+    return result;
+  }
+
   std::string str() const {
     // 200 bytes should be more than sufficient to hold a 100-digit number
     // plus additional bytes for the decimal point, '-' sign etc.
@@ -328,6 +358,19 @@
   }
 }
 
+template <typename InputType>
+cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
+ternaryOperationOneOutput(Operation op, InputType x, InputType y, InputType z) {
+  constexpr unsigned int prec = Precision<InputType>::value;
+  MPFRNumber inputX(x, prec), inputY(y, prec), inputZ(z, prec);
+  switch (op) {
+  case Operation::Fma:
+    return inputX.fma(inputY, inputZ);
+  default:
+    __builtin_unreachable();
+  }
+}
+
 template <typename T>
 void explainUnaryOperationSingleOutputError(Operation op, T input, T matchValue,
                                             testutils::StreamWrapper &OS) {
@@ -476,6 +519,48 @@
     Operation, const BinaryInput<long double> &, long double,
     testutils::StreamWrapper &);
 
+template <typename T>
+void explainTernaryOperationOneOutputError(Operation op,
+                                           const TernaryInput<T> &input,
+                                           T libcResult,
+                                           testutils::StreamWrapper &OS) {
+  MPFRNumber mpfrX(input.x, Precision<T>::value);
+  MPFRNumber mpfrY(input.y, Precision<T>::value);
+  MPFRNumber mpfrZ(input.z, Precision<T>::value);
+  FPBits<T> xbits(input.x);
+  FPBits<T> ybits(input.y);
+  FPBits<T> zbits(input.z);
+  MPFRNumber mpfrResult =
+      ternaryOperationOneOutput(op, input.x, input.y, input.z);
+  MPFRNumber mpfrMatchValue(libcResult);
+
+  OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str()
+     << " z: " << mpfrZ.str() << '\n';
+  __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
+                                              OS);
+  __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
+                                              OS);
+  __llvm_libc::fputil::testing::describeValue("Third input bits: ", input.z,
+                                              OS);
+
+  OS << "Libc result: " << mpfrMatchValue.str() << '\n'
+     << "MPFR result: " << mpfrResult.str() << '\n';
+  __llvm_libc::fputil::testing::describeValue(
+      "Libc floating point result bits: ", libcResult, OS);
+  __llvm_libc::fputil::testing::describeValue(
+      "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
+  OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
+}
+
+template void explainTernaryOperationOneOutputError<float>(
+    Operation, const TernaryInput<float> &, float, testutils::StreamWrapper &);
+template void explainTernaryOperationOneOutputError<double>(
+    Operation, const TernaryInput<double> &, double,
+    testutils::StreamWrapper &);
+template void explainTernaryOperationOneOutputError<long double>(
+    Operation, const TernaryInput<long double> &, long double,
+    testutils::StreamWrapper &);
+
 template <typename T>
 bool compareUnaryOperationSingleOutput(Operation op, T input, T libcResult,
                                        double ulpError) {
@@ -575,6 +660,27 @@
 template bool compareBinaryOperationOneOutput<long double>(
     Operation, const BinaryInput<long double> &, long double, double);
 
+template <typename T>
+bool compareTernaryOperationOneOutput(Operation op,
+                                      const TernaryInput<T> &input,
+                                      T libcResult, double ulpError) {
+  MPFRNumber mpfrResult =
+      ternaryOperationOneOutput(op, input.x, input.y, input.z);
+  double ulp = mpfrResult.ulp(libcResult);
+
+  bool bitsAreEven = ((FPBits<T>(libcResult).bitsAsUInt() & 1) == 0);
+  return (ulp < ulpError) ||
+         ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
+}
+
+template bool
+compareTernaryOperationOneOutput<float>(Operation, const TernaryInput<float> &,
+                                        float, double);
+template bool compareTernaryOperationOneOutput<double>(
+    Operation, const TernaryInput<double> &, double, double);
+template bool compareTernaryOperationOneOutput<long double>(
+    Operation, const TernaryInput<long double> &, long double, double);
+
 static mpfr_rnd_t getMPFRRoundingMode(RoundingMode mode) {
   switch (mode) {
   case RoundingMode::Upward: