diff --git a/libc/config/linux/api.td b/libc/config/linux/api.td
--- a/libc/config/linux/api.td
+++ b/libc/config/linux/api.td
@@ -199,6 +199,12 @@
    "modfl",
    "expf",
    "exp2f",
+   "remainderf",
+   "remainder",
+   "remainderl",
+   "remquof",
+   "remquo",
+   "remquol",
    "round",
    "roundf",
    "roundl",
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
@@ -103,6 +103,12 @@
     libc.src.math.modf
     libc.src.math.modff
     libc.src.math.modfl
+    libc.src.math.remainderf
+    libc.src.math.remainder
+    libc.src.math.remainderl
+    libc.src.math.remquof
+    libc.src.math.remquo
+    libc.src.math.remquol
     libc.src.math.round
     libc.src.math.roundf
     libc.src.math.roundl
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -310,6 +310,14 @@
           FunctionSpec<"expf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"exp2f", 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>]>,
+
           FunctionSpec<"round", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"roundf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"roundl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>]>,
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
@@ -521,3 +521,75 @@
   COMPILE_OPTIONS
     -O2
 )
+
+add_entrypoint_object(
+  remquof
+  SRCS
+    remquof.cpp
+  HDRS
+    remquof.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
+
+add_entrypoint_object(
+  remquo
+  SRCS
+    remquo.cpp
+  HDRS
+    remquo.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
+
+add_entrypoint_object(
+  remquol
+  SRCS
+    remquol.cpp
+  HDRS
+    remquol.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
+
+add_entrypoint_object(
+  remainderf
+  SRCS
+    remainderf.cpp
+  HDRS
+    remainderf.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
+
+add_entrypoint_object(
+  remainder
+  SRCS
+    remainder.cpp
+  HDRS
+    remainder.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
+
+add_entrypoint_object(
+  remainderl
+  SRCS
+    remainderl.cpp
+  HDRS
+    remainderl.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
diff --git a/libc/src/math/remainder.h b/libc/src/math/remainder.h
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remainder.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for remainder ---------------------*- 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_REMAINDER_H
+#define LLVM_LIBC_SRC_MATH_REMAINDER_H
+
+namespace __llvm_libc {
+
+double remainder(double x, double y);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_REMAINDER_H
diff --git a/libc/src/math/remainder.cpp b/libc/src/math/remainder.cpp
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remainder.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of remainder 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/DivisionAndRemainderOperations.h"
+
+namespace __llvm_libc {
+
+double LLVM_LIBC_ENTRYPOINT(remainder)(double x, double y) {
+  int quotient;
+  return fputil::remquo(x, y, quotient);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/remainderf.h b/libc/src/math/remainderf.h
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remainderf.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for remainderf --------------------*- 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_REMAINDERF_H
+#define LLVM_LIBC_SRC_MATH_REMAINDERF_H
+
+namespace __llvm_libc {
+
+float remainderf(float x, float y);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_REMAINDERF_H
diff --git a/libc/src/math/remainderf.cpp b/libc/src/math/remainderf.cpp
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remainderf.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of remainderf 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/DivisionAndRemainderOperations.h"
+
+namespace __llvm_libc {
+
+float LLVM_LIBC_ENTRYPOINT(remainderf)(float x, float y) {
+  int quotient;
+  return fputil::remquo(x, y, quotient);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/remainderl.h b/libc/src/math/remainderl.h
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remainderl.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for remainderl --------------------*- 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_REMAINDERL_H
+#define LLVM_LIBC_SRC_MATH_REMAINDERL_H
+
+namespace __llvm_libc {
+
+long double remainderl(long double x, long double y);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_REMAINDERL_H
diff --git a/libc/src/math/remainderl.cpp b/libc/src/math/remainderl.cpp
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remainderl.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of remainderl 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/DivisionAndRemainderOperations.h"
+
+namespace __llvm_libc {
+
+long double LLVM_LIBC_ENTRYPOINT(remainderl)(long double x, long double y) {
+  int quotient;
+  return fputil::remquo(x, y, quotient);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/remquo.h b/libc/src/math/remquo.h
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remquo.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for remquo ------------------------*- 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_REMQUO_H
+#define LLVM_LIBC_SRC_MATH_REMQUO_H
+
+namespace __llvm_libc {
+
+double remquo(double x, double y, int *exp);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_REMQUO_H
diff --git a/libc/src/math/remquo.cpp b/libc/src/math/remquo.cpp
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remquo.cpp
@@ -0,0 +1,18 @@
+//===-- Implementation of remquo 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/DivisionAndRemainderOperations.h"
+
+namespace __llvm_libc {
+
+double LLVM_LIBC_ENTRYPOINT(remquo)(double x, double y, int *exp) {
+  return fputil::remquo(x, y, *exp);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/remquof.h b/libc/src/math/remquof.h
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remquof.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for remquof -----------------------*- 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_REMQUOF_H
+#define LLVM_LIBC_SRC_MATH_REMQUOF_H
+
+namespace __llvm_libc {
+
+float remquof(float x, float y, int *exp);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_REMQUOF_H
diff --git a/libc/src/math/remquof.cpp b/libc/src/math/remquof.cpp
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remquof.cpp
@@ -0,0 +1,18 @@
+//===-- Implementation of remquof 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/DivisionAndRemainderOperations.h"
+
+namespace __llvm_libc {
+
+float LLVM_LIBC_ENTRYPOINT(remquof)(float x, float y, int *exp) {
+  return fputil::remquo(x, y, *exp);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/remquol.h b/libc/src/math/remquol.h
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remquol.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for remquol -----------------------*- 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_REMQUOL_H
+#define LLVM_LIBC_SRC_MATH_REMQUOL_H
+
+namespace __llvm_libc {
+
+long double remquol(long double x, long double y, int *exp);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_REMQUOL_H
diff --git a/libc/src/math/remquol.cpp b/libc/src/math/remquol.cpp
new file mode 100644
--- /dev/null
+++ b/libc/src/math/remquol.cpp
@@ -0,0 +1,19 @@
+//===-- Implementation of remquol 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/DivisionAndRemainderOperations.h"
+
+namespace __llvm_libc {
+
+long double LLVM_LIBC_ENTRYPOINT(remquol)(long double x, long double y,
+                                          int *exp) {
+  return fputil::remquo(x, y, *exp);
+}
+
+} // 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
@@ -552,3 +552,42 @@
     libc.src.math.sqrtl
     libc.utils.FPUtil.fputil
 )
+
+add_fp_unittest(
+  remquof_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    remquof_test.cpp
+  DEPENDS
+    libc.include.math
+    libc.src.math.remquof
+    libc.utils.FPUtil.fputil
+)
+
+add_fp_unittest(
+  remquo_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    remquo_test.cpp
+  DEPENDS
+    libc.include.math
+    libc.src.math.remquo
+    libc.utils.FPUtil.fputil
+)
+
+add_fp_unittest(
+  remquol_test
+  NEED_MPFR
+  SUITE
+    libc_math_unittests
+  SRCS
+    remquol_test.cpp
+  DEPENDS
+    libc.include.math
+    libc.src.math.remquol
+    libc.utils.FPUtil.fputil
+)
diff --git a/libc/test/src/math/remquo_test.cpp b/libc/test/src/math/remquo_test.cpp
new file mode 100644
--- /dev/null
+++ b/libc/test/src/math/remquo_test.cpp
@@ -0,0 +1,91 @@
+//===-- Unittests for remquo ----------------------------------------------===//
+//
+// 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 "include/math.h"
+#include "src/math/remquo.h"
+#include "utils/FPUtil/BasicOperations.h"
+#include "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+using FPBits = __llvm_libc::fputil::FPBits<double>;
+using UIntType = FPBits::UIntType;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+static const float zero = FPBits::zero();
+static const float negZero = FPBits::negZero();
+static const float nan = FPBits::buildNaN(1);
+static const float inf = FPBits::inf();
+static const float negInf = FPBits::negInf();
+
+TEST(RemquoTest, SpecialNumbers) {
+  int exponent;
+  double x, y;
+
+  y = 1.0;
+  x = inf;
+  EXPECT_NE(isnan(__llvm_libc::remquo(x, y, &exponent)), 0);
+  x = negInf;
+  EXPECT_NE(isnan(__llvm_libc::remquo(x, y, &exponent)), 0);
+
+  x = 1.0;
+  y = zero;
+  EXPECT_NE(isnan(__llvm_libc::remquo(x, y, &exponent)), 0);
+  y = negZero;
+  EXPECT_NE(isnan(__llvm_libc::remquo(x, y, &exponent)), 0);
+
+  y = nan;
+  x = 1.0;
+  EXPECT_NE(isnan(__llvm_libc::remquo(x, y, &exponent)), 0);
+
+  y = 1.0;
+  x = nan;
+  EXPECT_NE(isnan(__llvm_libc::remquo(x, y, &exponent)), 0);
+
+  x = nan;
+  y = nan;
+  EXPECT_NE(isnan(__llvm_libc::remquo(x, y, &exponent)), 0);
+
+  x = zero;
+  y = 1.0;
+  EXPECT_FP_EQ(__llvm_libc::remquo(x, y, &exponent), zero);
+
+  x = negZero;
+  y = 1.0;
+  EXPECT_FP_EQ(__llvm_libc::remquo(x, y, &exponent), negZero);
+}
+
+TEST(RemquoTest, SubnormalRange) {
+  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) {
+    double x = FPBits(v), y = FPBits(w);
+    mpfr::BinaryOutput<double> result;
+    mpfr::BinaryInput<double> input{x, y};
+    result.f = __llvm_libc::remquo(x, y, &result.i);
+    ASSERT_MPFR_MATCH(mpfr::Operation::RemQuo, input, result, 0.0);
+  }
+}
+
+TEST(RemquoTest, NormalRange) {
+  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) {
+    double x = FPBits(v), y = FPBits(w);
+    mpfr::BinaryOutput<double> result;
+    mpfr::BinaryInput<double> input{x, y};
+    result.f = __llvm_libc::remquo(x, y, &result.i);
+    ASSERT_MPFR_MATCH(mpfr::Operation::RemQuo, input, result, 0.0);
+  }
+}
diff --git a/libc/test/src/math/remquof_test.cpp b/libc/test/src/math/remquof_test.cpp
new file mode 100644
--- /dev/null
+++ b/libc/test/src/math/remquof_test.cpp
@@ -0,0 +1,91 @@
+//===-- Unittests for remquof ---------------------------------------------===//
+//
+// 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 "include/math.h"
+#include "src/math/remquof.h"
+#include "utils/FPUtil/BasicOperations.h"
+#include "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+using FPBits = __llvm_libc::fputil::FPBits<float>;
+using UIntType = FPBits::UIntType;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+static const float zero = FPBits::zero();
+static const float negZero = FPBits::negZero();
+static const float nan = FPBits::buildNaN(1);
+static const float inf = FPBits::inf();
+static const float negInf = FPBits::negInf();
+
+TEST(RemquofTest, SpecialNumbers) {
+  int exponent;
+  float x, y;
+
+  y = 1.0f;
+  x = inf;
+  EXPECT_NE(isnan(__llvm_libc::remquof(x, y, &exponent)), 0);
+  x = negInf;
+  EXPECT_NE(isnan(__llvm_libc::remquof(x, y, &exponent)), 0);
+
+  x = 1.0f;
+  y = zero;
+  EXPECT_NE(isnan(__llvm_libc::remquof(x, y, &exponent)), 0);
+  y = negZero;
+  EXPECT_NE(isnan(__llvm_libc::remquof(x, y, &exponent)), 0);
+
+  y = nan;
+  x = 1.0f;
+  EXPECT_NE(isnan(__llvm_libc::remquof(x, y, &exponent)), 0);
+
+  y = 1.0f;
+  x = nan;
+  EXPECT_NE(isnan(__llvm_libc::remquof(x, y, &exponent)), 0);
+
+  x = nan;
+  y = nan;
+  EXPECT_NE(isnan(__llvm_libc::remquof(x, y, &exponent)), 0);
+
+  x = zero;
+  y = 1.0f;
+  EXPECT_FP_EQ(__llvm_libc::remquof(x, y, &exponent), zero);
+
+  x = negZero;
+  y = 1.0f;
+  EXPECT_FP_EQ(__llvm_libc::remquof(x, y, &exponent), negZero);
+}
+
+TEST(RemquofTest, SubnormalRange) {
+  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) {
+    float x = FPBits(v), y = FPBits(w);
+    mpfr::BinaryOutput<float> result;
+    mpfr::BinaryInput<float> input{x, y};
+    result.f = __llvm_libc::remquof(x, y, &result.i);
+    ASSERT_MPFR_MATCH(mpfr::Operation::RemQuo, input, result, 0.0);
+  }
+}
+
+TEST(RemquofTest, NormalRange) {
+  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) {
+    float x = FPBits(v), y = FPBits(w);
+    mpfr::BinaryOutput<float> result;
+    mpfr::BinaryInput<float> input{x, y};
+    result.f = __llvm_libc::remquof(x, y, &result.i);
+    ASSERT_MPFR_MATCH(mpfr::Operation::RemQuo, input, result, 0.0);
+  }
+}
diff --git a/libc/test/src/math/remquol_test.cpp b/libc/test/src/math/remquol_test.cpp
new file mode 100644
--- /dev/null
+++ b/libc/test/src/math/remquol_test.cpp
@@ -0,0 +1,97 @@
+//===-- Unittests for remquol ---------------------------------------------===//
+//
+// 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 "include/math.h"
+#include "src/math/remquol.h"
+#include "utils/FPUtil/BasicOperations.h"
+#include "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include "utils/UnitTest/Test.h"
+
+using FPBits = __llvm_libc::fputil::FPBits<long double>;
+using UIntType = FPBits::UIntType;
+
+namespace mpfr = __llvm_libc::testing::mpfr;
+
+static const long double zero = FPBits::zero();
+static const long double negZero = FPBits::negZero();
+static const long double nan = FPBits::buildNaN(1);
+static const long double inf = FPBits::inf();
+static const long double negInf = FPBits::negInf();
+
+TEST(RemquoTest, SpecialNumbers) {
+  int exponent;
+  long double x, y;
+
+  y = 1.0l;
+  x = inf;
+  EXPECT_NE(isnan(__llvm_libc::remquol(x, y, &exponent)), 0);
+  x = negInf;
+  EXPECT_NE(isnan(__llvm_libc::remquol(x, y, &exponent)), 0);
+
+  x = 1.0l;
+  y = zero;
+  EXPECT_NE(isnan(__llvm_libc::remquol(x, y, &exponent)), 0);
+  y = negZero;
+  EXPECT_NE(isnan(__llvm_libc::remquol(x, y, &exponent)), 0);
+
+  y = nan;
+  x = 1.0l;
+  EXPECT_NE(isnan(__llvm_libc::remquol(x, y, &exponent)), 0);
+
+  y = 1.0l;
+  x = nan;
+  EXPECT_NE(isnan(__llvm_libc::remquol(x, y, &exponent)), 0);
+
+  x = nan;
+  y = nan;
+  EXPECT_NE(isnan(__llvm_libc::remquol(x, y, &exponent)), 0);
+
+  x = zero;
+  y = 1.0l;
+  EXPECT_FP_EQ(__llvm_libc::remquol(x, y, &exponent), zero);
+
+  x = negZero;
+  y = 1.0l;
+  EXPECT_FP_EQ(__llvm_libc::remquol(x, y, &exponent), negZero);
+}
+
+TEST(RemquofTest, SubnormalRange) {
+  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) {
+    long double x = FPBits(v), y = FPBits(w);
+    mpfr::BinaryOutput<long double> result;
+    mpfr::BinaryInput<long double> input{x, y};
+    result.f = __llvm_libc::remquol(x, y, &result.i);
+    ASSERT_MPFR_MATCH(mpfr::Operation::RemQuo, input, result, 0.0);
+  }
+}
+
+TEST(RemquofTest, NormalRange) {
+  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) {
+    long double x = FPBits(v), y = FPBits(w);
+    mpfr::BinaryOutput<long double> result;
+    result.f = __llvm_libc::remquol(x, y, &result.i);
+    // In normal range on x86 platforms, the implicit 1 bit can be zero making
+    // the numbers NaN. Hence we test for them separately.
+    if (isnan(x) || isnan(y)) {
+      ASSERT_NE(isnan(result.f), 0);
+    } else {
+      mpfr::BinaryInput<long double> input{x, y};
+      ASSERT_MPFR_MATCH(mpfr::Operation::RemQuo, input, result, 0.0);
+    }
+  }
+}
diff --git a/libc/utils/FPUtil/CMakeLists.txt b/libc/utils/FPUtil/CMakeLists.txt
--- a/libc/utils/FPUtil/CMakeLists.txt
+++ b/libc/utils/FPUtil/CMakeLists.txt
@@ -11,6 +11,7 @@
     BasicOperations.h
     BitPatterns.h
     ClassificationFunctions.h
+    DivisionAndRemainderOperations.h
     FloatOperations.h
     FloatProperties.h
     FPBits.h
diff --git a/libc/utils/FPUtil/DivisionAndRemainderOperations.h b/libc/utils/FPUtil/DivisionAndRemainderOperations.h
new file mode 100644
--- /dev/null
+++ b/libc/utils/FPUtil/DivisionAndRemainderOperations.h
@@ -0,0 +1,111 @@
+//===-- Floating point divsion and remainder operations ---------*- 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_DIVISION_AND_REMAINDER_OPERATIONS_H
+#define LLVM_LIBC_UTILS_FPUTIL_DIVISION_AND_REMAINDER_OPERATIONS_H
+
+#include "FPBits.h"
+#include "ManipulationFunctions.h"
+#include "NormalFloat.h"
+
+#include "utils/CPP/TypeTraits.h"
+
+namespace __llvm_libc {
+namespace fputil {
+
+static constexpr int quotientLSBBits = 3;
+
+// The implementation is a bit-by-bit algorithm which uses integer division
+// to evaluate the quotient and remainder.
+template <typename T,
+          cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
+static inline T remquo(T x, T y, int &q) {
+  FPBits<T> xbits(x), ybits(y);
+  if (xbits.isNaN())
+    return x;
+  if (ybits.isNaN())
+    return y;
+  if (xbits.isInf() || ybits.isZero())
+    return FPBits<T>::buildNaN(1);
+
+  if (xbits.isZero() || ybits.isInf()) {
+    q = 0;
+    return __llvm_libc::fputil::copysign(T(0.0), x);
+  }
+
+  bool resultSign = (xbits.sign == ybits.sign ? false : true);
+
+  // Once we know the sign of the result, we can just operate on the absolute
+  // values. The correct sign can be applied to the result after the result
+  // is evaluated.
+  xbits.sign = ybits.sign = 0;
+
+  NormalFloat<T> normalx(xbits), normaly(ybits);
+  int exp = normalx.exponent - normaly.exponent;
+  typename NormalFloat<T>::UIntType mx = normalx.mantissa,
+                                    my = normaly.mantissa;
+
+  q = 0;
+  while (exp >= 0) {
+    unsigned shiftCount = 0;
+    typename NormalFloat<T>::UIntType n = mx;
+    for (shiftCount = 0; n < my; n <<= 1, ++shiftCount)
+      ;
+
+    if (static_cast<int>(shiftCount) > exp)
+      break;
+
+    exp -= shiftCount;
+    if (0 <= exp && exp < quotientLSBBits)
+      q |= (1 << exp);
+
+    mx = n - my;
+    if (mx == 0)
+      return __llvm_libc::fputil::copysign(T(0.0), x);
+  }
+
+  NormalFloat<T> remainder(exp + normaly.exponent, mx, 0);
+
+  // Since NormalFloat to native type conversion is a truncation operation
+  // currently, the remainder value in the native type is correct as is.
+  // However, if NormalFloat to native type conversion is updated in future,
+  // then the conversion to native remainder value should be updated
+  // appropriately and some directed tests added.
+  T nativeRemainder(remainder);
+  T absy = T(ybits);
+  int cmp = remainder.mul2(1).cmp(normaly);
+  if (cmp > 0) {
+    q = q + 1;
+    if (x >= T(0.0))
+      nativeRemainder = nativeRemainder - absy;
+    else
+      nativeRemainder = absy - nativeRemainder;
+  } else if (cmp == 0) {
+    if (q & 1) {
+      q += 1;
+      if (x >= T(0.0))
+        nativeRemainder = -nativeRemainder;
+    } else {
+      if (x < T(0.0))
+        nativeRemainder = -nativeRemainder;
+    }
+  } else {
+    if (x < T(0.0))
+      nativeRemainder = -nativeRemainder;
+  }
+
+  q = resultSign ? -q : q;
+  if (nativeRemainder == T(0.0))
+    return __llvm_libc::fputil::copysign(T(0.0), x);
+  return nativeRemainder;
+}
+
+} // namespace fputil
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_UTILS_FPUTIL_DIVISION_AND_REMAINDER_OPERATIONS_H
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
@@ -73,6 +73,14 @@
   static constexpr int exponentBias = (1 << (ExponentWidth<T>::value - 1)) - 1;
   static constexpr int maxExponent = (1 << ExponentWidth<T>::value) - 1;
 
+  static constexpr UIntType minSubnormal = UIntType(1);
+  static constexpr UIntType maxSubnormal =
+      (UIntType(1) << MantissaWidth<T>::value) - 1;
+  static constexpr UIntType minNormal =
+      (UIntType(1) << MantissaWidth<T>::value);
+  static constexpr UIntType maxNormal =
+      ((UIntType(maxExponent) - 1) << MantissaWidth<T>::value) | maxSubnormal;
+
   // We don't want accidental type promotions/conversions so we require exact
   // type match.
   template <typename XType,
diff --git a/libc/utils/FPUtil/LongDoubleBitsX86.h b/libc/utils/FPUtil/LongDoubleBitsX86.h
--- a/libc/utils/FPUtil/LongDoubleBitsX86.h
+++ b/libc/utils/FPUtil/LongDoubleBitsX86.h
@@ -33,6 +33,15 @@
 
   static constexpr int exponentBias = 0x3FFF;
   static constexpr int maxExponent = 0x7FFF;
+  static constexpr UIntType minSubnormal = UIntType(1);
+  // Subnormal numbers include the implicit bit in x86 long double formats.
+  static constexpr UIntType maxSubnormal =
+      (UIntType(1) << (MantissaWidth<long double>::value + 1)) - 1;
+  static constexpr UIntType minNormal =
+      (UIntType(3) << MantissaWidth<long double>::value);
+  static constexpr UIntType maxNormal =
+      ((UIntType(maxExponent) - 1) << (MantissaWidth<long double>::value + 1)) |
+      (UIntType(1) << MantissaWidth<long double>::value) | maxSubnormal;
 
   UIntType mantissa : MantissaWidth<long double>::value;
   uint8_t implicitBit : 1;