[libc] Add implementations of nextafter[f|l] functions.

A differential fuzzer for these functions has also been added.
Along the way, a small correction has been done to the normal/subnormal
limits of x86 long double values.

Reviewed By: lntue

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

21 files changed, 618 insertions, 1 deletions

diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 67e2b72..b904262 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -88,6 +88,9 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.modf
     libc.src.math.modff
     libc.src.math.modfl
+    libc.src.math.nextafter
+    libc.src.math.nextafterf
+    libc.src.math.nextafterl
     libc.src.math.remainderf
     libc.src.math.remainder
     libc.src.math.remainderl
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 45571ee..a34c5964 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -144,6 +144,9 @@ set(TARGET_LIBM_ENTRYPOINTS
     libc.src.math.nearbyint
     libc.src.math.nearbyintf
     libc.src.math.nearbyintl
+    libc.src.math.nextafter
+    libc.src.math.nextafterf
+    libc.src.math.nextafterl
     libc.src.math.remainderf
     libc.src.math.remainder
     libc.src.math.remainderl
diff --git a/libc/fuzzing/math/CMakeLists.txt b/libc/fuzzing/math/CMakeLists.txt
index 7c710fc..af9a5d3 100644
--- a/libc/fuzzing/math/CMakeLists.txt
+++ b/libc/fuzzing/math/CMakeLists.txt
@@ -48,3 +48,15 @@ add_libc_fuzzer(
     libc.utils.FPUtil.fputil
     libc.utils.CPP.standalone_cpp
 )
+
+add_libc_fuzzer(
+  nextafter_differential_fuzz
+  SRCS
+    nextafter_differential_fuzz.cpp
+  HDRS
+    TwoInputSingleOutputDiff.h
+  DEPENDS
+    libc.src.math.nextafter
+    libc.src.math.nextafterf
+    libc.src.math.nextafterl
+)
diff --git a/libc/fuzzing/math/Compare.h b/libc/fuzzing/math/Compare.h
index 62ab6c7..11d5465 100644
--- a/libc/fuzzing/math/Compare.h
+++ b/libc/fuzzing/math/Compare.h
@@ -10,6 +10,7 @@
 #define LLVM_LIBC_FUZZING_MATH_COMPARE_H
 
 #include "utils/CPP/TypeTraits.h"
+#include "utils/FPUtil/FPBits.h"
 
 template <typename T>
 __llvm_libc::cpp::EnableIfType<__llvm_libc::cpp::IsFloatingPointType<T>::Value,
diff --git a/libc/fuzzing/math/nextafter_differential_fuzz.cpp b/libc/fuzzing/math/nextafter_differential_fuzz.cpp
new file mode 100644
index 0000000..f4a7891
--- /dev/null
+++ b/libc/fuzzing/math/nextafter_differential_fuzz.cpp
@@ -0,0 +1,26 @@
+//===-- nextafter_differential_fuzz.cpp
+//---------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// Differential fuzz test for llvm-libc nextafter implementation.
+///
+//===----------------------------------------------------------------------===//
+
+#include "fuzzing/math/TwoInputSingleOutputDiff.h"
+
+#include "src/math/nextafter.h"
+#include "src/math/nextafterf.h"
+#include "src/math/nextafterl.h"
+
+extern "C" int LLVMFuzzerTestOneInput(const uint8_t *data, size_t size) {
+  TwoInputSingleOutputDiff<float, float>(&__llvm_libc::nextafterf,
+                                         &::nextafterf, data, size);
+  TwoInputSingleOutputDiff<double, double>(&__llvm_libc::nextafter,
+                                           &::nextafter, data, size);
+  return 0;
+}
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index 355321c..41f6083 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -394,6 +394,10 @@ def StdC : StandardSpec<"stdc"> {
           FunctionSpec<"nearbyint", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
           FunctionSpec<"nearbyintf", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
           FunctionSpec<"nearbyintl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>]>,
+
+          FunctionSpec<"nextafterf", RetValSpec<FloatType>, [ArgSpec<FloatType>, ArgSpec<FloatType>]>,
+          FunctionSpec<"nextafter", RetValSpec<DoubleType>, [ArgSpec<DoubleType>, ArgSpec<DoubleType>]>,
+          FunctionSpec<"nextafterl", RetValSpec<LongDoubleType>, [ArgSpec<LongDoubleType>, ArgSpec<LongDoubleType>]>,
       ]
   >;
 
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index aa4b4a5..15a2e46 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -941,3 +941,40 @@ add_entrypoint_object(
   COMPILE_OPTIONS
     -O2
 )
+
+add_entrypoint_object(
+  nextafter
+  SRCS
+    nextafter.cpp
+  HDRS
+    nextafter.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
+
+add_entrypoint_object(
+  nextafterf
+  SRCS
+    nextafterf.cpp
+  HDRS
+    nextafterf.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
+
+add_entrypoint_object(
+  nextafterl
+  SRCS
+    nextafterl.cpp
+  HDRS
+    nextafterl.h
+  DEPENDS
+    libc.utils.FPUtil.fputil
+  COMPILE_OPTIONS
+    -O2
+)
+
diff --git a/libc/src/math/nextafter.cpp b/libc/src/math/nextafter.cpp
new file mode 100644
index 0000000..8ba1553
--- /dev/null
+++ b/libc/src/math/nextafter.cpp
@@ -0,0 +1,18 @@
+//===-- Implementation of nextafter 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/ManipulationFunctions.h"
+
+namespace __llvm_libc {
+
+double LLVM_LIBC_ENTRYPOINT(nextafter)(double x, double y) {
+  return fputil::nextafter(x, y);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/nextafter.h b/libc/src/math/nextafter.h
new file mode 100644
index 0000000..1f6d24d
--- /dev/null
+++ b/libc/src/math/nextafter.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for nextafter ---------------------*- 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_NEXTAFTER_H
+#define LLVM_LIBC_SRC_MATH_NEXTAFTER_H
+
+namespace __llvm_libc {
+
+double nextafter(double x, double y);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_NEXTAFTER_H
diff --git a/libc/src/math/nextafterf.cpp b/libc/src/math/nextafterf.cpp
new file mode 100644
index 0000000..dab591e
--- /dev/null
+++ b/libc/src/math/nextafterf.cpp
@@ -0,0 +1,18 @@
+//===-- Implementation of nextafterf 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/ManipulationFunctions.h"
+
+namespace __llvm_libc {
+
+float LLVM_LIBC_ENTRYPOINT(nextafterf)(float x, float y) {
+  return fputil::nextafter(x, y);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/nextafterf.h b/libc/src/math/nextafterf.h
new file mode 100644
index 0000000..8afac0e
--- /dev/null
+++ b/libc/src/math/nextafterf.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for nextafterf --------------------*- 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_NEXTAFTERF_H
+#define LLVM_LIBC_SRC_MATH_NEXTAFTERF_H
+
+namespace __llvm_libc {
+
+float nextafterf(float x, float y);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_NEXTAFTERF_H
diff --git a/libc/src/math/nextafterl.cpp b/libc/src/math/nextafterl.cpp
new file mode 100644
index 0000000..4a9b64c
--- /dev/null
+++ b/libc/src/math/nextafterl.cpp
@@ -0,0 +1,18 @@
+//===-- Implementation of nextafterl 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/ManipulationFunctions.h"
+
+namespace __llvm_libc {
+
+long double LLVM_LIBC_ENTRYPOINT(nextafterl)(long double x, long double y) {
+  return fputil::nextafter(x, y);
+}
+
+} // namespace __llvm_libc
diff --git a/libc/src/math/nextafterl.h b/libc/src/math/nextafterl.h
new file mode 100644
index 0000000..9128778
--- /dev/null
+++ b/libc/src/math/nextafterl.h
@@ -0,0 +1,18 @@
+//===-- Implementation header for nextafterl --------------------*- 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_NEXTAFTERL_H
+#define LLVM_LIBC_SRC_MATH_NEXTAFTERL_H
+
+namespace __llvm_libc {
+
+long double nextafterl(long double x, long double y);
+
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_SRC_MATH_NEXTAFTERL_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index f07843a..3f4d8c7 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1007,3 +1007,45 @@ add_fp_unittest(
     libc.src.math.hypot
     libc.utils.FPUtil.fputil
 )
+
+add_fp_unittest(
+  nextafter_test
+  SUITE
+    libc_math_unittests
+  SRCS
+    nextafter_test.cpp
+  HDRS
+    NextAfterTest.h
+  DEPENDS
+    libc.include.math
+    libc.src.math.nextafter
+    libc.utils.FPUtil.fputil
+)
+
+add_fp_unittest(
+  nextafterf_test
+  SUITE
+    libc_math_unittests
+  SRCS
+    nextafterf_test.cpp
+  HDRS
+    NextAfterTest.h
+  DEPENDS
+    libc.include.math
+    libc.src.math.nextafterf
+    libc.utils.FPUtil.fputil
+)
+
+add_fp_unittest(
+  nextafterl_test
+  SUITE
+    libc_math_unittests
+  SRCS
+    nextafterl_test.cpp
+  HDRS
+    NextAfterTest.h
+  DEPENDS
+    libc.include.math
+    libc.src.math.nextafterl
+    libc.utils.FPUtil.fputil
+)
diff --git a/libc/test/src/math/NextAfterTest.h b/libc/test/src/math/NextAfterTest.h
new file mode 100644
index 0000000..6286aa9
--- /dev/null
+++ b/libc/test/src/math/NextAfterTest.h
@@ -0,0 +1,193 @@
+//===-- Utility class to test different flavors of nextafter ----*- 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_TEST_SRC_MATH_NEXTAFTERTEST_H
+#define LLVM_LIBC_TEST_SRC_MATH_NEXTAFTERTEST_H
+
+#include "utils/CPP/TypeTraits.h"
+#include "utils/FPUtil/BasicOperations.h"
+#include "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/TestHelpers.h"
+#include "utils/UnitTest/Test.h"
+#include <math.h>
+
+template <typename T>
+class NextAfterTestTemplate : public __llvm_libc::testing::Test {
+  using FPBits = __llvm_libc::fputil::FPBits<T>;
+  using MantissaWidth = __llvm_libc::fputil::MantissaWidth<T>;
+  using UIntType = typename FPBits::UIntType;
+
+#if (defined(__x86_64__) || defined(__i386__))
+  static constexpr int bitWidthOfType =
+      __llvm_libc::cpp::IsSame<T, long double>::Value ? 80 : (sizeof(T) * 8);
+#else
+  static constexpr int bitWidthOfType = sizeof(T) * 8;
+#endif
+
+  const T zero = FPBits::zero();
+  const T negZero = FPBits::negZero();
+  const T inf = FPBits::inf();
+  const T negInf = FPBits::negInf();
+  const T nan = FPBits::buildNaN(1);
+  const UIntType minSubnormal = FPBits::minSubnormal;
+  const UIntType maxSubnormal = FPBits::maxSubnormal;
+  const UIntType minNormal = FPBits::minNormal;
+  const UIntType maxNormal = FPBits::maxNormal;
+
+public:
+  typedef T (*NextAfterFunc)(T, T);
+
+  void testNaN(NextAfterFunc func) {
+    ASSERT_FP_EQ(func(nan, 0), nan);
+    ASSERT_FP_EQ(func(0, nan), nan);
+  }
+
+  void testBoundaries(NextAfterFunc func) {
+    ASSERT_FP_EQ(func(zero, negZero), negZero);
+    ASSERT_FP_EQ(func(negZero, zero), zero);
+
+    // 'from' is zero|negZero.
+    T x = zero;
+    T result = func(x, T(1));
+    UIntType expectedBits = 1;
+    T expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    result = func(x, T(-1));
+    expectedBits = (UIntType(1) << (bitWidthOfType - 1)) + 1;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    x = negZero;
+    result = func(x, 1);
+    expectedBits = 1;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    result = func(x, -1);
+    expectedBits = (UIntType(1) << (bitWidthOfType - 1)) + 1;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    // 'from' is max subnormal value.
+    x = *reinterpret_cast<const T *>(&maxSubnormal);
+    result = func(x, 1);
+    expected = *reinterpret_cast<const T *>(&minNormal);
+    ASSERT_FP_EQ(result, expected);
+
+    result = func(x, 0);
+    expectedBits = maxSubnormal - 1;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    x = -x;
+
+    result = func(x, -1);
+    expectedBits = (UIntType(1) << (bitWidthOfType - 1)) + minNormal;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    result = func(x, 0);
+    expectedBits = (UIntType(1) << (bitWidthOfType - 1)) + maxSubnormal - 1;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    // 'from' is min subnormal value.
+    x = *reinterpret_cast<const T *>(&minSubnormal);
+    result = func(x, 1);
+    expectedBits = minSubnormal + 1;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+    ASSERT_FP_EQ(func(x, 0), 0);
+
+    x = -x;
+    result = func(x, -1);
+    expectedBits = (UIntType(1) << (bitWidthOfType - 1)) + minSubnormal + 1;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+    ASSERT_FP_EQ(func(x, 0), T(-0.0));
+
+    // 'from' is min normal.
+    x = *reinterpret_cast<const T *>(&minNormal);
+    result = func(x, 0);
+    expectedBits = maxSubnormal;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    result = func(x, inf);
+    expectedBits = minNormal + 1;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    x = -x;
+    result = func(x, 0);
+    expectedBits = (UIntType(1) << (bitWidthOfType - 1)) + maxSubnormal;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    result = func(x, -inf);
+    expectedBits = (UIntType(1) << (bitWidthOfType - 1)) + minNormal + 1;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+
+    // 'from' is max normal and 'to' is infinity.
+    x = *reinterpret_cast<const T *>(&maxNormal);
+    result = func(x, inf);
+    ASSERT_FP_EQ(result, inf);
+
+    result = func(-x, -inf);
+    ASSERT_FP_EQ(result, -inf);
+
+    // 'from' is infinity.
+    x = inf;
+    result = func(x, 0);
+    expectedBits = maxNormal;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+    ASSERT_FP_EQ(func(x, inf), inf);
+
+    x = negInf;
+    result = func(x, 0);
+    expectedBits = (UIntType(1) << (bitWidthOfType - 1)) + maxNormal;
+    expected = *reinterpret_cast<T *>(&expectedBits);
+    ASSERT_FP_EQ(result, expected);
+    ASSERT_FP_EQ(func(x, negInf), negInf);
+
+    // 'from' is a power of 2.
+    x = T(32.0);
+    result = func(x, 0);
+    FPBits xBits = FPBits(x);
+    FPBits resultBits = FPBits(result);
+    ASSERT_EQ(resultBits.exponent, uint16_t(xBits.exponent - 1));
+    ASSERT_EQ(resultBits.mantissa, (UIntType(1) << MantissaWidth::value) - 1);
+
+    result = func(x, T(33.0));
+    resultBits = FPBits(result);
+    ASSERT_EQ(resultBits.exponent, xBits.exponent);
+    ASSERT_EQ(resultBits.mantissa, xBits.mantissa + UIntType(1));
+
+    x = -x;
+
+    result = func(x, 0);
+    resultBits = FPBits(result);
+    ASSERT_EQ(resultBits.exponent, uint16_t(xBits.exponent - 1));
+    ASSERT_EQ(resultBits.mantissa, (UIntType(1) << MantissaWidth::value) - 1);
+
+    result = func(x, T(-33.0));
+    resultBits = FPBits(result);
+    ASSERT_EQ(resultBits.exponent, xBits.exponent);
+    ASSERT_EQ(resultBits.mantissa, xBits.mantissa + UIntType(1));
+  }
+};
+
+#define LIST_NEXTAFTER_TESTS(T, func)                                          \
+  using NextAfterTest = NextAfterTestTemplate<T>;                              \
+  TEST_F(NextAfterTest, TestNaN) { testNaN(&func); }                           \
+  TEST_F(NextAfterTest, TestBoundaries) { testBoundaries(&func); }
+
+#endif // LLVM_LIBC_TEST_SRC_MATH_NEXTAFTERTEST_H
diff --git a/libc/test/src/math/nextafter_test.cpp b/libc/test/src/math/nextafter_test.cpp
new file mode 100644
index 0000000..2dde48c
--- /dev/null
+++ b/libc/test/src/math/nextafter_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for nextafter -------------------------------------------===//
+//
+// 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 "NextAfterTest.h"
+
+#include "src/math/nextafter.h"
+
+LIST_NEXTAFTER_TESTS(double, __llvm_libc::nextafter)
diff --git a/libc/test/src/math/nextafterf_test.cpp b/libc/test/src/math/nextafterf_test.cpp
new file mode 100644
index 0000000..14234c7
--- /dev/null
+++ b/libc/test/src/math/nextafterf_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for nextafterf ------------------------------------------===//
+//
+// 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 "NextAfterTest.h"
+
+#include "src/math/nextafterf.h"
+
+LIST_NEXTAFTER_TESTS(float, __llvm_libc::nextafterf)
diff --git a/libc/test/src/math/nextafterl_test.cpp b/libc/test/src/math/nextafterl_test.cpp
new file mode 100644
index 0000000..db85d83
--- /dev/null
+++ b/libc/test/src/math/nextafterl_test.cpp
@@ -0,0 +1,13 @@
+//===-- Unittests for nextafterl ------------------------------------------===//
+//
+// 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 "NextAfterTest.h"
+
+#include "src/math/nextafterl.h"
+
+LIST_NEXTAFTER_TESTS(long double, __llvm_libc::nextafterl)
diff --git a/libc/utils/FPUtil/LongDoubleBitsX86.h b/libc/utils/FPUtil/LongDoubleBitsX86.h
index fdde1bd..c051660 100644
--- a/libc/utils/FPUtil/LongDoubleBitsX86.h
+++ b/libc/utils/FPUtil/LongDoubleBitsX86.h
@@ -36,7 +36,7 @@ template <> struct __attribute__((packed)) FPBits<long double> {
   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;
+      (UIntType(1) << (MantissaWidth<long double>::value)) - 1;
   static constexpr UIntType minNormal =
       (UIntType(3) << MantissaWidth<long double>::value);
   static constexpr UIntType maxNormal =
diff --git a/libc/utils/FPUtil/ManipulationFunctions.h b/libc/utils/FPUtil/ManipulationFunctions.h
index 79dc741..f230a46 100644
--- a/libc/utils/FPUtil/ManipulationFunctions.h
+++ b/libc/utils/FPUtil/ManipulationFunctions.h
@@ -143,7 +143,42 @@ static inline T ldexp(T x, int exp) {
   return normal;
 }
 
+template <typename T,
+          cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, int> = 0>
+static inline T nextafter(T from, T to) {
+  FPBits<T> fromBits(from);
+  if (fromBits.isNaN())
+    return from;
+
+  FPBits<T> toBits(to);
+  if (toBits.isNaN())
+    return to;
+
+  if (from == to)
+    return to;
+
+  using UIntType = typename FPBits<T>::UIntType;
+  auto intVal = fromBits.bitsAsUInt();
+  UIntType signMask = (UIntType(1) << (sizeof(T) * 8 - 1));
+  if (from != T(0.0)) {
+    if ((from < to) == (from > T(0.0))) {
+      ++intVal;
+    } else {
+      --intVal;
+    }
+  } else {
+    intVal = (toBits.bitsAsUInt() & signMask) + UIntType(1);
+  }
+
+  return *reinterpret_cast<T *>(&intVal);
+  // TODO: Raise floating point exceptions as required by the standard.
+}
+
 } // namespace fputil
 } // namespace __llvm_libc
 
+#if (defined(__x86_64__) || defined(__i386__))
+#include "NextAfterLongDoubleX86.h"
+#endif // defined(__x86_64__) || defined(__i386__)
+
 #endif // LLVM_LIBC_UTILS_FPUTIL_MANIPULATION_FUNCTIONS_H
diff --git a/libc/utils/FPUtil/NextAfterLongDoubleX86.h b/libc/utils/FPUtil/NextAfterLongDoubleX86.h
new file mode 100644
index 0000000..e1555b7
--- /dev/null
+++ b/libc/utils/FPUtil/NextAfterLongDoubleX86.h
@@ -0,0 +1,114 @@
+//===-- nextafter implementation for x86 long double numbers ----*- 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_NEXT_AFTER_LONG_DOUBLE_X86_H
+#define LLVM_LIBC_UTILS_FPUTIL_NEXT_AFTER_LONG_DOUBLE_X86_H
+
+#include "FPBits.h"
+
+#include <stdint.h>
+
+namespace __llvm_libc {
+namespace fputil {
+
+static inline long double nextafter(long double from, long double to) {
+  using FPBits = FPBits<long double>;
+  FPBits fromBits(from);
+  if (fromBits.isNaN())
+    return from;
+
+  FPBits toBits(to);
+  if (toBits.isNaN())
+    return to;
+
+  if (from == to)
+    return to;
+
+  // Convert pseudo subnormal number to normal number.
+  if (fromBits.implicitBit == 1 && fromBits.exponent == 0) {
+    fromBits.exponent = 1;
+  }
+
+  using UIntType = FPBits::UIntType;
+  constexpr UIntType signVal = (UIntType(1) << 79);
+  constexpr UIntType mantissaMask =
+      (UIntType(1) << MantissaWidth<long double>::value) - 1;
+  auto intVal = fromBits.bitsAsUInt();
+  if (from < 0.0l) {
+    if (from > to) {
+      if (intVal == (signVal + FPBits::maxSubnormal)) {
+        // We deal with normal/subnormal boundary separately to avoid
+        // dealing with the implicit bit.
+        intVal = signVal + FPBits::minNormal;
+      } else if ((intVal & mantissaMask) == mantissaMask) {
+        fromBits.mantissa = 0;
+        // Incrementing exponent might overflow the value to infinity,
+        // which is what is expected. Since NaNs are handling separately,
+        // it will never overflow "beyond" infinity.
+        ++fromBits.exponent;
+        return fromBits;
+      } else {
+        ++intVal;
+      }
+    } else {
+      if (intVal == (signVal + FPBits::minNormal)) {
+        // We deal with normal/subnormal boundary separately to avoid
+        // dealing with the implicit bit.
+        intVal = signVal + FPBits::maxSubnormal;
+      } else if ((intVal & mantissaMask) == 0) {
+        fromBits.mantissa = mantissaMask;
+        // from == 0 is handled separately so decrementing the exponent will not
+        // lead to underflow.
+        --fromBits.exponent;
+        return fromBits;
+      } else {
+        --intVal;
+      }
+    }
+  } else if (from == 0.0l) {
+    if (from > to)
+      intVal = signVal + 1;
+    else
+      intVal = 1;
+  } else {
+    if (from > to) {
+      if (intVal == FPBits::minNormal) {
+        intVal = FPBits::maxSubnormal;
+      } else if ((intVal & mantissaMask) == 0) {
+        fromBits.mantissa = mantissaMask;
+        // from == 0 is handled separately so decrementing the exponent will not
+        // lead to underflow.
+        --fromBits.exponent;
+        return fromBits;
+      } else {
+        --intVal;
+      }
+    } else {
+      if (intVal == FPBits::maxSubnormal) {
+        intVal = FPBits::minNormal;
+      } else if ((intVal & mantissaMask) == mantissaMask) {
+        fromBits.mantissa = 0;
+        // Incrementing exponent might overflow the value to infinity,
+        // which is what is expected. Since NaNs are handling separately,
+        // it will never overflow "beyond" infinity.
+        ++fromBits.exponent;
+        return fromBits;
+      } else {
+        ++intVal;
+      }
+    }
+  }
+
+  return *reinterpret_cast<long double *>(&intVal);
+  // TODO: Raise floating point exceptions as required by the standard.
+}
+
+} // namespace fputil
+} // namespace __llvm_libc
+
+#endif // LLVM_LIBC_UTILS_FPUTIL_NEXT_AFTER_LONG_DOUBLE_X86_H