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

[sivachandra]

Summary

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.

[lntue]

libc/utils/FPUtil/DivisionAndRemainderOperations.h

73

When the remainder is in the subnormal range, I think the final check below is correct, but with an implicit assumption that the rounding from NormalFloat<T> to T is truncation.
Can you add that to comments, in case we change the rounding NormalFloat<T> to T to the default round-to-nearest, tie-to-even? So that at least we will add more tests for subnormal range to make sure that it still works correctly, even when NormalFloat<T> rounding is updated.

[sivachandra]

Add comment as suggested.

[sivachandra]

Added comment as suggested.

libc/config/linux/api.td

def MathAPI : PublicAPI<"math.h"> {
    "logb",
    "logbf",
    "logbl",
    "modf",
    "modff",
    "modfl",
    "expf",
    "exp2f",
+   "remainderf",
+   "remainder",
+   "remainderl",
+   "remquof",
+   "remquo",
+   "remquol",
    "round",
    "roundf",
    "roundl",
    "sincosf",
    "sinf",
    "sqrt",
    "sqrtf",
    "sqrtl",

libc/config/linux/x86_64/entrypoints.txt

set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.frexpf
     libc.src.math.frexpl
     libc.src.math.logb
     libc.src.math.logbf
     libc.src.math.logbl
     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.src.math.sincosf
     libc.src.math.sinf
     libc.src.math.sqrt
     libc.src.math.sqrtf
     libc.src.math.sqrtl
     libc.src.math.trunc
     libc.src.math.truncf
     libc.src.math.truncl
 )

libc/spec/stdc.td

HeaderSpec Math = HeaderSpec<
           FunctionSpec<"modfl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoublePtr>]>,
 
           FunctionSpec<"cosf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"sinf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
 
           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>]>,
 
           FunctionSpec<"sqrt", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"sqrtf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"sqrtl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>]>,
 

libc/src/math/CMakeLists.txt

SRCS
     sqrtl.cpp
   HDRS
     sqrtl.h
   DEPENDS
     libc.utils.FPUtil.fputil
   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.h

+//===-- 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/remainder.cpp

+//===-- 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/remainderf.h

+//===-- 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/remainderf.cpp

+//===-- 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/remainderl.h

+//===-- 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/remainderl.cpp

+//===-- 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/remquo.h

+//===-- 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/remquo.cpp

+//===-- 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/remquof.h

+//===-- 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/remquof.cpp

+//===-- 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/remquol.h

+//===-- 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/src/math/remquol.cpp

+//===-- 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/test/src/math/CMakeLists.txt

SUITE
     libc_math_unittests
   SRCS
     sqrtl_test.cpp
   DEPENDS
     libc.include.math
     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

+//===-- 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

+//===-- 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

+//===-- 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

 if(${LIBC_TARGET_MACHINE} MATCHES "^x86.*")
   set(LONG_DOUBLE_HDR LongDoubleBitsX86.h)
 else()
   set(LONG_DOUBLE_HDR)
 endif()
 
 add_header_library(
   fputil
   HDRS
     ${LONG_DOUBLE_HDR}
     BasicOperations.h
     BitPatterns.h
     ClassificationFunctions.h
+    DivisionAndRemainderOperations.h
     FloatOperations.h
     FloatProperties.h
     FPBits.h
     BasicOperations.h
     ManipulationFunctions.h
     NearestIntegerOperations.h
     NormalFloat.h
   DEPENDS

libc/utils/FPUtil/DivisionAndRemainderOperations.h

+//===-- 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);
+  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

template <typename T> struct __attribute__((packed)) FPBits {
 
   UIntType mantissa : MantissaWidth<T>::value;
   uint16_t exponent : ExponentWidth<T>::value;
   uint8_t sign : 1;
 
   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,
             cpp::EnableIfType<cpp::IsSame<T, XType>::Value, int> = 0>
   explicit FPBits(XType x) {
     *this = *reinterpret_cast<FPBits<T> *>(&x);
   }
 

libc/utils/FPUtil/LongDoubleBitsX86.h

 // x86_64 padding.
 template <> struct Padding<8> { static constexpr unsigned value = 48; };
 
 template <> struct __attribute__((packed)) FPBits<long double> {
   using UIntType = __uint128_t;
 
   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;
   uint16_t exponent : ExponentWidth<long double>::value;
   uint8_t sign : 1;
   uint64_t padding : Padding<sizeof(uintptr_t)>::value;
 
   template <typename XType,