[libc] Add implementations of remquo[f|l] and remainder[f|l].

The implementation is not fully standards compliant in the sense that
errno is not set on error, and floating point exceptions are not raised.

Subnormal range and normal range are tested separately in the tests.

Reviewed By: lntue

Differential Revision: https://reviews.llvm.org/D86666

GitOrigin-RevId: 8514ecb02d4330bc075b9c8fef77c87810088d2f
diff --git a/config/linux/api.td b/config/linux/api.td
index 063fe40..33ae64c 100644
--- a/config/linux/api.td
+++ b/config/linux/api.td
@@ -199,6 +199,12 @@
    "modfl",
    "expf",
    "exp2f",
+   "remainderf",
+   "remainder",
+   "remainderl",
+   "remquof",
+   "remquo",
+   "remquol",
    "round",
    "roundf",
    "roundl",
diff --git a/config/linux/x86_64/entrypoints.txt b/config/linux/x86_64/entrypoints.txt
index c24173b..6aca5e4 100644
--- a/config/linux/x86_64/entrypoints.txt
+++ b/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/spec/stdc.td b/spec/stdc.td
index 15fc12d..77fa971 100644
--- a/spec/stdc.td
+++ b/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/src/math/CMakeLists.txt b/src/math/CMakeLists.txt
index 0c878de..3b4f821 100644
--- a/src/math/CMakeLists.txt
+++ b/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/src/math/remainder.cpp b/src/math/remainder.cpp
new file mode 100644
index 0000000..880e6a6
--- /dev/null
+++ b/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/src/math/remainder.h b/src/math/remainder.h
new file mode 100644
index 0000000..8a720fc
--- /dev/null
+++ b/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/src/math/remainderf.cpp b/src/math/remainderf.cpp
new file mode 100644
index 0000000..bab3201
--- /dev/null
+++ b/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/src/math/remainderf.h b/src/math/remainderf.h
new file mode 100644
index 0000000..19a16d0
--- /dev/null
+++ b/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/src/math/remainderl.cpp b/src/math/remainderl.cpp
new file mode 100644
index 0000000..bd9bc49
--- /dev/null
+++ b/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/src/math/remainderl.h b/src/math/remainderl.h
new file mode 100644
index 0000000..f283763
--- /dev/null
+++ b/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/src/math/remquo.cpp b/src/math/remquo.cpp
new file mode 100644
index 0000000..b61d7d4
--- /dev/null
+++ b/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/src/math/remquo.h b/src/math/remquo.h
new file mode 100644
index 0000000..cb753fe
--- /dev/null
+++ b/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/src/math/remquof.cpp b/src/math/remquof.cpp
new file mode 100644
index 0000000..246bee0
--- /dev/null
+++ b/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/src/math/remquof.h b/src/math/remquof.h
new file mode 100644
index 0000000..feb2e4f
--- /dev/null
+++ b/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/src/math/remquol.cpp b/src/math/remquol.cpp
new file mode 100644
index 0000000..8e02876
--- /dev/null
+++ b/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/src/math/remquol.h b/src/math/remquol.h
new file mode 100644
index 0000000..d1b0e20
--- /dev/null
+++ b/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/test/src/math/CMakeLists.txt b/test/src/math/CMakeLists.txt
index 07b5052..e1bac1a 100644
--- a/test/src/math/CMakeLists.txt
+++ b/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/test/src/math/remquo_test.cpp b/test/src/math/remquo_test.cpp
new file mode 100644
index 0000000..4ea61dd
--- /dev/null
+++ b/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/test/src/math/remquof_test.cpp b/test/src/math/remquof_test.cpp
new file mode 100644
index 0000000..0c51d5f
--- /dev/null
+++ b/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/test/src/math/remquol_test.cpp b/test/src/math/remquol_test.cpp
new file mode 100644
index 0000000..eab3a5f
--- /dev/null
+++ b/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/utils/FPUtil/CMakeLists.txt b/utils/FPUtil/CMakeLists.txt
index 745ede3..8a6cc36 100644
--- a/utils/FPUtil/CMakeLists.txt
+++ b/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/utils/FPUtil/DivisionAndRemainderOperations.h b/utils/FPUtil/DivisionAndRemainderOperations.h
new file mode 100644
index 0000000..ceae538
--- /dev/null
+++ b/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/utils/FPUtil/FPBits.h b/utils/FPUtil/FPBits.h
index 2c630db..89bdd92 100644
--- a/utils/FPUtil/FPBits.h
+++ b/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/utils/FPUtil/LongDoubleBitsX86.h b/utils/FPUtil/LongDoubleBitsX86.h
index 482c4de..fdde1bd 100644
--- a/utils/FPUtil/LongDoubleBitsX86.h
+++ b/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;