[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

libc/config/linux/api.td

@@ -199,6 +199,12 @@ def MathAPI : PublicAPI<"math.h"> {
    "modfl",
    "expf",
    "exp2f",
+   "remainderf",
+   "remainder",
+   "remainderl",
+   "remquof",
+   "remquo",
+   "remquol",
    "round",
    "roundf",
    "roundl",

libc/config/linux/x86_64/entrypoints.txt

@@ -103,6 +103,12 @@ set(TARGET_LIBM_ENTRYPOINTS
     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

libc/spec/stdc.td

@@ -310,6 +310,14 @@ def StdC : StandardSpec<"stdc"> {
           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>]>,

libc/src/math/CMakeLists.txt

@@ -521,3 +521,75 @@ add_entrypoint_object(
   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
+)

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

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

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

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

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

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

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

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

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

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

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

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

libc/test/src/math/CMakeLists.txt

@@ -552,3 +552,42 @@ add_fp_unittest(
     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
+)

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);
+  }
+}

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);
+  }
+}

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);
+    }
+  }
+}

libc/utils/FPUtil/CMakeLists.txt

@@ -11,6 +11,7 @@ add_header_library(
     BasicOperations.h
     BitPatterns.h
     ClassificationFunctions.h
+    DivisionAndRemainderOperations.h
     FloatOperations.h
     FloatProperties.h
     FPBits.h

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

libc/utils/FPUtil/FPBits.h

@@ -73,6 +73,14 @@ template <typename T> struct __attribute__((packed)) FPBits {
   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,

libc/utils/FPUtil/LongDoubleBitsX86.h

@@ -33,6 +33,15 @@ template <> struct __attribute__((packed)) FPBits<long double> {
 
   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;