Added implementations of five Mathematical Special Functions added in C++17.

[thebrandre]

Summary

Added implementations of std::assoc_laguerre, std::assoc_legendre, std::hermite, std::laguerre, and std::legendre as specified in ISO/IEC JTC 1/SC 22/WG 21 N3060 and in the C++17 standard.

The implementations were tested against the implementations of Visual Studio 2017 and libstdc++.
I got the same results up to floating-point precision with one exception:
libstdc++'s implementation of std::assoc_legendre doesn't fully adhere to the standard: it has the (pretty trivial) Condon-Shortley phase term, which is not part of the standard and explicitly excluded on cppreference.com.
Visual Studio 2017 complies to the standard in this respect.
I plan to submit that change to libstdc++ for consistency.

All float versions use double internally! It turns out that float is too small for, e.g, std::assoc_legendre in the range of m, l up to 127. You can pretty easily get NaN with float if you don't want to check for inf in each iteration.
The double implementation doesn't run into this problem in the range specified by the standard.

[mclow.lists]

Thanks for doing this; I'll be commenting on this today. Sorry for the lag.

[thebrandre]

• Fixed compilation error of unit tests with clang (previously only tested with gcc).
• Removed the checks for STDCPP_WANT_MATH_SPEC_FUNCS as specified in ISO/IEC JTC 1/SC 22/WG 21 N3060. The special functions are thus available if and only if C++17 is enabled.
• Adjusted Doxygen style to LLVM Coding Standards (backslash instead of @)

[thebrandre]

Minor changes:

• Applied clang-format and LLVM Coding Standards to all files.
• Added unit tests for std::hermite. We now check the roots of the first 20 Hermite polynomials against values from the literature. And also check if the recurrence relation holds

[thebrandre]

I am going to resubmit the functions separately ... and will integrate it into the project just like https://reviews.llvm.org/D59937

libcxx/include/experimental/__hermite

+//===------------------------ __hermite -------------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the internal implementations of std::hermite.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef _LIBCPP_EXPERIMENTAL___HERMITE
+#define _LIBCPP_EXPERIMENTAL___HERMITE
+
+#include <experimental/__config>
+#include <cmath>
+#include <limits>
+
+/// \return the hermite polynomial \f$ H_{n}(x) \f$
+/// \note The implementation is based on the recurrence formula
+/// \f[
+/// nH_{n+1}(x) = 2x H_{n}(x) - 2 n H_{n-1}
+/// \f]
+/// Press, William H., et al. Numerical recipes 3rd edition: The art of
+/// scientific computing. Cambridge university press, 2007, p. 182.
+template <class _Real>
+_Real __libcpp_hermite_recurrence(unsigned __n, _Real __x) {
+  if (__n == 0u)
+    return _Real(1);
+
+  _Real __t2(1);
+  _Real __t1 = _Real(2) * __x;
+  for (unsigned __i = 1; __i < __n; ++__i) {
+    const _Real __t0 = _Real(2) * (__x * __t1 - _Real(__i) * __t2);
+    __t2 = __t1;
+    __t1 = __t0;
+  }
+  return __t1;
+}
+
+template <class _Real> _Real __libcpp_hermite(unsigned __n, _Real __x) {
+  if (std::isnan(__x))
+    return std::numeric_limits<_Real>::quiet_NaN();
+
+  return __libcpp_hermite_recurrence(__n, __x);
+}
+
+#endif // _LIBCPP_EXPERIMENTAL___HERMITE

libcxx/include/experimental/__laguerre

+//===------------------------ __laguerre ------------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the internal implementations of std::laguerre
+/// and std::assoc_laguerre.
+///
+//===----------------------------------------------------------------------===//
+
+
+#ifndef _LIBCPP_EXPERIMENTAL___LAGUERRE
+#define _LIBCPP_EXPERIMENTAL___LAGUERRE
+
+#include <experimental/__config>
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+
+/// \return the generalized laguerre polynomial \f$ L_{n}^{(\alpha)}(x) \f$
+/// \note The implementation is based on the recurrence formula
+/// \f[
+/// nL_{n}^{(\alpha)}(x) = (-x + 2n + \alpha - 1) L_{n-1}^{(\alpha)}(x) -
+/// (n + \alpha - 1) L_{n-2}^{(\alpha)}(x)
+/// \f]
+/// Press, William H., et al. Numerical recipes 3rd edition: The art of
+/// scientific computing. Cambridge university press, 2007, p. 182.
+template <class _Real>
+_Real __libcpp_generalized_laguerre_recurrence(unsigned __n, _Real __alpha,
+                                               _Real __x) {
+  if (__n == 0u)
+    return _Real(1);
+
+  _Real __delta = __alpha - __x;
+  _Real __li = _Real(1) + __delta;
+  const _Real __alpham1 = __alpha - _Real(1);
+  for (unsigned __i = 2; __i <= __n; ++__i) {
+    __delta = (__delta * (_Real(__i) + __alpham1) - __x * __li) / _Real(__i);
+    __li += __delta;
+  }
+  return __li;
+}
+
+template <class _Real>
+_Real __libcpp_assoc_laguerre(unsigned __n, unsigned __m, _Real __x) {
+  if (std::isnan(__x))
+    return std::numeric_limits<_Real>::quiet_NaN();
+
+  if (__x < _Real(0))
+    _VSTD::__throw_domain_error(
+        "Argument of assoc_laguerre function is out of range");
+
+  return __libcpp_generalized_laguerre_recurrence(__n, _Real(__m), __x);
+}
+
+template <class _Real> _Real __libcpp_laguerre(unsigned __n, _Real __x) {
+  if (std::isnan(__x))
+    return std::numeric_limits<_Real>::quiet_NaN();
+
+  if (__x < _Real(0))
+    _VSTD::__throw_domain_error(
+        "Argument of laguerre function is out of range");
+
+  return __libcpp_generalized_laguerre_recurrence(__n, _Real(0), __x);
+}
+
+#endif // _LIBCPP_EXPERIMENTAL___LAGUERRE

libcxx/include/experimental/__legendre

+//===------------------------ __legendre ------------------------*- 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the internal implementations of std::legendre
+/// and std::assoc_legendre.
+///
+//===----------------------------------------------------------------------===//
+
+
+#ifndef _LIBCPP_EXPERIMENTAL___LEGENDRE
+#define _LIBCPP_EXPERIMENTAL___LEGENDRE
+
+#include <experimental/__config>
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+
+/// \return the Legendre polynomial \f$ P_{n}(x) \f$
+/// \note The implementation is based on the recurrence formula
+/// \f[
+/// (n+1)P_{n+1}(x) = (2n+1)xP_{n}(x) - nP_{n-1}(x)
+/// \f]
+/// Press, William H., et al. Numerical recipes 3rd edition: The art of
+/// scientific computing. Cambridge university press, 2007, p. 182.
+template <class _Real>
+_Real __libcpp_legendre_recurrence(unsigned __n, _Real __x) {
+  if (__n == 0u)
+    return _Real(1);
+
+  _Real __t2(1);
+  _Real __t1 = __x;
+  for (unsigned __i = 1; __i < __n; ++__i) {
+    const _Real __k = _Real(__i);
+    _Real __t0 = ((_Real(2) * __k + _Real(1)) * __x * __t1 - __k * __t2) /
+                 (__k + _Real(1));
+    __t2 = __t1;
+    __t1 = __t0;
+  }
+  return __t1;
+}
+
+template <class _Real> _Real __libcpp_legendre(unsigned __n, _Real __x) {
+  if (std::isnan(__x))
+    return std::numeric_limits<_Real>::quiet_NaN();
+
+  if (std::abs(__x) > _Real(1))
+    _VSTD::__throw_domain_error(
+        "Argument of legendre function is out of range");
+
+  return __libcpp_legendre_recurrence(__n, __x);
+}
+
+/// \return \f$ s^{-m} P_{l}^{m}(x) \f$ with an additonal scaling factor to
+/// prevent overflow. \note The implementation is based on the recurrence
+/// formula \f[ (l-m+1)P_{l+1}^{m}(x) = (2l+1)xP_{l}^{m}(x) -
+/// (l+m)P_{l-1}^{m}(x) \f] with \f[ P_{m}^{m}(x) = \sqrt{1 -
+/// x^2}^{m}\frac{(2m)!}{2^m m!} \f] and \f[ P_{m-1}^{m}(x) = 0 \f] \attention
+/// The starting point of the recursion grows exponentially with __m! For large
+/// m, we have the following relation: \f[ P_{m}^{m}(x) \approx \sqrt{1 -
+/// x^2}^{m}\sqrt{2} 2^{n} \exp( n(\ln n - 1 )) \f] For example, for \f$ m = 40
+/// \f$, we already have \f$ P_{40}^{40}(0) \approx 8 \cdot 10^{58}  \f$
+/// \attention The so-called Condon-Shortley phase term is omitted in the C++17
+/// standard's definition of std::assoc_laguerre.
+template <class _Real>
+_Real __libcpp_assoc_legendre_recurrence(unsigned __l, unsigned __m, _Real __x,
+                                         _Real __scale = _Real(1)) {
+  if (__m == 0u)
+    return __libcpp_legendre_recurrence(__l, __x);
+
+  if (__l < __m)
+    return _Real(0);
+
+  if (__l == 0u)
+    return _Real(1);
+
+  _Real __pmm = _Real(1);
+  // Note: (1-x)*(1+x) is more accurate than (1-x*x)
+  // "What Every Computer Scientist Should Know About Floating-Point
+  // Arithmetic", David Goldberg, p. 38
+  const _Real __t =
+      std::sqrt((_Real(1) - __x) * (_Real(1) + __x)) / (_Real(2) * __scale);
+  for (unsigned __i = 2u * __m; __i > __m; --__i)
+    __pmm *= __t * __i;
+
+  if (__l == __m)
+    return __pmm;
+
+  // Actually, we'd start with _pmm but it grows exponentially with __m.
+  // Luckily, the recursion scales. So we can start with 1 and multiply
+  // afterwards.
+  _Real __t2 = _Real(1);
+  _Real __t1 = _Real(2u * __m + 1u) * __x; // first iteration unfolded
+  for (unsigned __i = __m + 1u; __i < __l; ++__i) {
+    // As soon as one of the terms becomes inf, this will quickly lead to NaNs.
+    // float just doesn't do it for the whole range up to l==127.
+    const _Real __t0 =
+        (_Real(2u * __i + 1u) * __x * __t1 - _Real(__i + __m) * __t2) /
+        _Real(__i - __m + 1u);
+    __t2 = __t1;
+    __t1 = __t0;
+  }
+  return __t1 * __pmm;
+}
+
+template <class _Real>
+_Real __libcpp_assoc_legendre(unsigned __n, unsigned __m, _Real __x) {
+  if (std::isnan(__x))
+    return std::numeric_limits<_Real>::quiet_NaN();
+
+  if (std::abs(__x) > _Real(1))
+    _VSTD::__throw_domain_error(
+        "Argument of assoc_legendre function is out of range");
+
+  return __libcpp_assoc_legendre_recurrence(__n, __m, __x);
+}
+
+#endif // _LIBCPP_EXPERIMENTAL___LEGENDRE

libcxx/include/experimental/cmath

+// -*- C++ -*-
+//===---------------------------- cmath ----------------------------------===//
+//
+// 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 _LIBCPP_EXPERIMENTAL_CMATH
+#define _LIBCPP_EXPERIMENTAL_CMATH
+
+// The following macro name shall be conditionally defined by the implementation
+// to indicate conformance to the International Standard ISO/IEC JTC 1/SC 22/WG 21 N3060
+//#define __STDCPP_MATH_SPEC_FUNCS__ 201003L
+
+#include <experimental/__config>
+#include <cmath>
+
+#if _LIBCPP_STD_VER > 14
+
+#include <experimental/__hermite>
+#include <experimental/__laguerre>
+#include <experimental/__legendre>
+
+_LIBCPP_BEGIN_NAMESPACE_EXPERIMENTAL
+
+inline _LIBCPP_INLINE_VISIBILITY double assoc_laguerre(unsigned __lcpp_n, unsigned __lcpp_m, double __lcpp_x)
+{
+    return __libcpp_assoc_laguerre<double>(__lcpp_n , __lcpp_m, __lcpp_x);
+}
+inline _LIBCPP_INLINE_VISIBILITY float assoc_laguerref(unsigned __lcpp_n, unsigned __lcpp_m, float __lcpp_x)
+{
+    return static_cast<float>(__libcpp_assoc_laguerre<double>(__lcpp_n , __lcpp_m, static_cast<double>(__lcpp_x)));
+}
+inline _LIBCPP_INLINE_VISIBILITY long double assoc_laguerrel(unsigned __lcpp_n, unsigned __lcpp_m, long double __lcpp_x)
+{
+    return __libcpp_assoc_laguerre<long double>(__lcpp_n , __lcpp_m, __lcpp_x);
+}
+
+
+inline _LIBCPP_INLINE_VISIBILITY double assoc_legendre(unsigned __lcpp_n, unsigned __lcpp_m, double __lcpp_x)
+{
+    return __libcpp_assoc_legendre<double>(__lcpp_n , __lcpp_m, __lcpp_x);
+}
+inline _LIBCPP_INLINE_VISIBILITY float assoc_legendref(unsigned __lcpp_n, unsigned __lcpp_m, float __lcpp_x)
+{
+    // use double internally -- float is too prone to overflow!
+    return static_cast<float>(__libcpp_assoc_legendre<double>(__lcpp_n , __lcpp_m, static_cast<double>(__lcpp_x)));
+}
+inline _LIBCPP_INLINE_VISIBILITY long double assoc_legendrel(unsigned __lcpp_n, unsigned __lcpp_m, long double __lcpp_x)
+{
+    return __libcpp_assoc_legendre<long double>(__lcpp_n , __lcpp_m, __lcpp_x);
+}
+
+
+inline _LIBCPP_INLINE_VISIBILITY double laguerre(unsigned __lcpp_n, double __lcpp_x)
+{
+    return __libcpp_laguerre<double>(__lcpp_n ,__lcpp_x);
+}
+inline _LIBCPP_INLINE_VISIBILITY float laguerref(unsigned __lcpp_n, float __lcpp_x)
+{
+    return static_cast<float>(__libcpp_laguerre<double>(__lcpp_n, static_cast<double>(__lcpp_x)));
+}
+inline _LIBCPP_INLINE_VISIBILITY long double laguerrel(unsigned __lcpp_n, long double __lcpp_x)
+{
+    return __libcpp_laguerre<long double>(__lcpp_n ,__lcpp_x);
+}
+
+
+inline _LIBCPP_INLINE_VISIBILITY double legendre(unsigned __lcpp_n, double __lcpp_x)
+{
+    return __libcpp_legendre<double>(__lcpp_n ,__lcpp_x);
+}
+inline _LIBCPP_INLINE_VISIBILITY float legendref(unsigned __lcpp_n, float __lcpp_x)
+{
+    return static_cast<float>(__libcpp_legendre<double>(__lcpp_n, static_cast<double>(__lcpp_x)));
+}
+inline _LIBCPP_INLINE_VISIBILITY long double legendrel(unsigned __lcpp_n, long double __lcpp_x)
+{
+    return __libcpp_legendre<long double>(__lcpp_n ,__lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY double hermite(unsigned __lcpp_n, double __lcpp_x)
+{
+    return __libcpp_hermite<double>(__lcpp_n ,__lcpp_x);
+}
+inline _LIBCPP_INLINE_VISIBILITY float hermitef(unsigned __lcpp_n, float __lcpp_x)
+{
+    // use double internally -- float is too prone to overflow!
+    return static_cast<float>(__libcpp_hermite(__lcpp_n , static_cast<double>(__lcpp_x)));
+}
+inline _LIBCPP_INLINE_VISIBILITY long double hermitel(unsigned __lcpp_n, long double __lcpp_x)
+{
+    return __libcpp_hermite<long double>(__lcpp_n ,__lcpp_x);
+}
+
+
+_LIBCPP_END_NAMESPACE_EXPERIMENTAL
+
+#endif // _LIBCPP_STD_VER > 14
+
+#endif // _LIBCPP_EXPERIMENTAL_CMATH

libcxx/test/libcxx/experimental/numerics/c.math/assoc_laguerre.pass.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
+//
+//===----------------------------------------------------------------------===//
+
+// <experimental/cmath>
+
+#include <cassert>
+#include <experimental/cmath>
+#include <limits>
+
+#if _LIBCPP_STD_VER > 14
+
+template <class T> void testAssocLaguerreNaNPropagation() {
+  const unsigned MaxN = 127;
+  const T x = std::numeric_limits<T>::quiet_NaN();
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    for (unsigned m = 0; m <= MaxN; ++m) {
+      assert(std::isnan(std::experimental::assoc_laguerre(n, m, x)));
+    }
+  }
+}
+
+template <class T> void testAssocLaguerreNotNaN(const T x) {
+  assert(!std::isnan(x));
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    for (unsigned m = 0; m <= MaxN; ++m) {
+      assert(!std::isnan(std::experimental::assoc_laguerre(n, m, x)));
+    }
+  }
+}
+
+template <class T> void testAssocLaguerreThrows(const T x) {
+#ifndef _LIBCPP_NO_EXCEPTIONS
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    for (unsigned m = 0; m <= MaxN; ++m) {
+      bool Throws = false;
+      try {
+        std::experimental::assoc_laguerre(n, m, x);
+      } catch (const std::domain_error &) {
+        Throws = true;
+      }
+      assert(Throws);
+    }
+  }
+#endif // _LIBCPP_NO_EXCEPTIONS
+}
+
+template <class T>
+void testAssocLaguerreVsLaguerre(const T x, const T AbsTolerance,
+                                 const T RelTolerance) {
+  assert(!std::isnan(x));
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    for (unsigned m = 0; m <= MaxN; ++m) {
+      const T Result = std::experimental::assoc_laguerre(n, 0, x);
+      const T ExpectedResult = std::experimental::laguerre(n, x);
+      const T Tolerance =
+          AbsTolerance + std::abs(ExpectedResult) * RelTolerance;
+      const T Difference = std::abs(Result - ExpectedResult);
+      assert(Difference <= Tolerance);
+    }
+  }
+}
+
+template <class T>
+void testAssocLaguerreAnalytic(const T x, const T AbsTolerance,
+                               const T RelTolerance) {
+  assert(!std::isnan(x));
+  const auto compareFloatingPoint =
+      [AbsTolerance, RelTolerance](const T Result, const T ExpectedResult) {
+        if (std::isinf(ExpectedResult) && std::isinf(Result))
+          return true;
+
+        if (std::isnan(ExpectedResult) || std::isnan(Result))
+          return false;
+
+        const T Tolerance =
+            AbsTolerance + std::abs(ExpectedResult) * RelTolerance;
+        return std::abs(Result - ExpectedResult) < Tolerance;
+      };
+
+  const auto l0 = [](T, unsigned) { return T(1); };
+  const auto l1 = [](T x, unsigned m) { return -x + T(m + 1); };
+  const auto l2 = [](T x, unsigned m) {
+    return x * x / T(2) - T(m + 2) * x + T(m + 1) * T(m + 2) / T(2);
+  };
+  const auto l3 = [](T x, unsigned m) {
+    return -x * x * x / T(6) + T(m + 3) * x * x / T(2) -
+           T(m + 2) * T(m + 3) * x / T(2) +
+           T(m + 1) * T(m + 2) * T(m + 3) / T(6);
+  };
+
+  for (unsigned m = 0; m < 128; ++m) {
+    assert(compareFloatingPoint(std::experimental::assoc_laguerre(0, m, x),
+                                l0(x, m)));
+    assert(compareFloatingPoint(std::experimental::assoc_laguerre(1, m, x),
+                                l1(x, m)));
+    assert(compareFloatingPoint(std::experimental::assoc_laguerre(2, m, x),
+                                l2(x, m)));
+    assert(compareFloatingPoint(std::experimental::assoc_laguerre(3, m, x),
+                                l3(x, m)));
+  }
+}
+
+template <class T>
+void testAssocLaguerre(const T AbsTolerance, const T RelTolerance) {
+  testAssocLaguerreNaNPropagation<T>();
+  testAssocLaguerreThrows<T>(T(-5));
+
+  const T Samples[] = {T(0.0), T(0.1), T(0.5), T(1.0), T(10.0)};
+
+  for (T x : Samples) {
+    testAssocLaguerreNotNaN(x);
+    testAssocLaguerreAnalytic(x, AbsTolerance, RelTolerance);
+    testAssocLaguerreVsLaguerre(x, AbsTolerance, RelTolerance);
+  }
+}
+
+#endif
+
+int main(int, char **) {
+#if _LIBCPP_STD_VER > 14
+  testAssocLaguerre<float>(1e-5f, 1e-5f);
+  testAssocLaguerre<double>(1e-9, 1e-9);
+  testAssocLaguerre<long double>(1e-12, 1e-12);
+#endif
+  return 0;
+}

libcxx/test/libcxx/experimental/numerics/c.math/assoc_legendre.pass.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
+//
+//===----------------------------------------------------------------------===//
+
+// <experimental/cmath>
+
+#include <cassert>
+#include <experimental/cmath>
+#include <limits>
+
+#if _LIBCPP_STD_VER > 14
+
+template <class T> void testAssocLegendreNaNPropagation() {
+  const unsigned MaxN = 127;
+  const T x = std::numeric_limits<T>::quiet_NaN();
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    for (unsigned m = 0; m <= MaxN; ++m) {
+      assert(std::isnan(std::experimental::assoc_legendre(n, m, x)));
+    }
+  }
+}
+
+template <class T> void testAssocLegendreNotNaN(const T x) {
+  assert(!std::isnan(x));
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    for (unsigned m = 0; m <= MaxN; ++m) {
+      assert(!std::isnan(std::experimental::assoc_legendre(n, m, x)));
+    }
+  }
+}
+
+template <class T> void testAssocLegendreThrows(const T x) {
+#ifndef _LIBCPP_NO_EXCEPTIONS
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    for (unsigned m = 0; m <= MaxN; ++m) {
+      bool Throws = false;
+      try {
+        std::experimental::assoc_legendre(n, m, x);
+      } catch (const std::domain_error &) {
+        Throws = true;
+      }
+      assert(Throws);
+    }
+  }
+#endif // _LIBCPP_NO_EXCEPTIONS
+}
+
+template <class T>
+void testAssocLegendreVsLegendre(const T x, const T AbsTolerance,
+                                 const T RelTolerance) {
+  assert(!std::isnan(x));
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    for (unsigned m = 0; m <= MaxN; ++m) {
+      const T Result = std::experimental::assoc_legendre(n, 0, x);
+      const T ExpectedResult = std::experimental::legendre(n, x);
+      const T Tolerance =
+          AbsTolerance + std::abs(ExpectedResult) * RelTolerance;
+      const T Difference = std::abs(Result - ExpectedResult);
+      assert(Difference <= Tolerance);
+    }
+  }
+}
+
+template <class T>
+void testAssocLegendreAnalytic(const T x, const T AbsTolerance,
+                               const T RelTolerance) {
+  assert(!std::isnan(x));
+  const auto compareFloatingPoint =
+      [AbsTolerance, RelTolerance](const T Result, const T ExpectedResult) {
+        if (std::isinf(ExpectedResult) && std::isinf(Result))
+          return true;
+
+        if (std::isnan(ExpectedResult) || std::isnan(Result))
+          return false;
+
+        const T Tolerance =
+            AbsTolerance + std::abs(ExpectedResult) * RelTolerance;
+        return std::abs(Result - ExpectedResult) < Tolerance;
+      };
+
+  const auto l00 = [](T) { return T(1); };
+
+  const auto l10 = [](T x) { return x; };
+  const auto l11 = [](T x) { return std::sqrt((T(1) - x) * (T(1) + x)); };
+
+  const auto l20 = [](T x) { return (T(3) * x * x - T(1)) / T(2); };
+  const auto l21 = [](T x) {
+    return T(3) * x * std::sqrt((T(1) - x) * (T(1) + x));
+  };
+  const auto l22 = [](T x) { return T(3) * (T(1) - x) * (T(1) + x); };
+
+  const auto l30 = [](T x) { return (T(5) * x * x - T(3)) * x / T(2); };
+  const auto l31 = [](T x) {
+    return T(3) / T(2) * (T(5) * x * x - T(1)) *
+           std::sqrt((T(1) - x) * (T(1) + x));
+  };
+  const auto l32 = [](T x) { return T(15) * x * (T(1) - x) * (T(1) + x); };
+  const auto l33 = [](T x) {
+    const T temp = (T(1) - x) * (T(1) + x);
+    return T(15) * temp * std::sqrt(temp);
+  };
+
+  const auto l40 = [](T x) {
+    return (T(35) * x * x * x * x - T(30) * x * x + T(3)) / T(8);
+  };
+  const auto l41 = [](T x) {
+    return T(5) / T(2) * x * (T(7) * x * x - T(3)) *
+           std::sqrt((T(1) - x) * (T(1) + x));
+  };
+  const auto l42 = [](T x) {
+    return T(15) / T(2) * (T(7) * x * x - 1) * (T(1) - x) * (T(1) + x);
+  };
+  const auto l43 = [](T x) {
+    const T temp = (T(1) - x) * (T(1) + x);
+    return T(105) * x * temp * std::sqrt(temp);
+  };
+  const auto l44 = [](T x) {
+    const T temp = (T(1) - x) * (T(1) + x);
+    return T(105) * temp * temp;
+  };
+
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(0, 0, x), l00(x)));
+
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(1, 0, x), l10(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(1, 1, x), l11(x)));
+
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(2, 0, x), l20(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(2, 1, x), l21(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(2, 2, x), l22(x)));
+
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(3, 0, x), l30(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(3, 1, x), l31(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(3, 2, x), l32(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(3, 3, x), l33(x)));
+
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(4, 0, x), l40(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(4, 1, x), l41(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(4, 2, x), l42(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(4, 3, x), l43(x)));
+  assert(
+      compareFloatingPoint(std::experimental::assoc_legendre(4, 4, x), l44(x)));
+
+  try {
+    const unsigned MaxN = 127;
+    for (unsigned n = 0; n <= MaxN; ++n) {
+      for (unsigned m = n + 1; m <= MaxN; ++m) {
+        assert(std::experimental::assoc_legendre(n, m, x) <= AbsTolerance);
+      }
+    }
+  } catch (const std::domain_error &) {
+    // Should not throw! The expression given in
+    // ISO/IEC JTC 1/SC 22/WG 21 N3060 is actually well-defined for m > n!
+    assert(false);
+  }
+}
+
+template <class T>
+void testAssocLegendre(const T AbsTolerance, const T RelTolerance) {
+  testAssocLegendreNaNPropagation<T>();
+  testAssocLegendreThrows<T>(T(-5));
+  testAssocLegendreThrows<T>(T(5));
+
+  const T Samples[] = {T(-1.0), T(-0.5), T(-0.1), T(0.0),
+                       T(0.1),  T(0.5),  T(1.0)};
+
+  for (T x : Samples) {
+    testAssocLegendreNotNaN(x);
+    testAssocLegendreVsLegendre(x, AbsTolerance, RelTolerance);
+    testAssocLegendreAnalytic(x, AbsTolerance, RelTolerance);
+  }
+}
+
+#endif
+
+int main(int, char **) {
+#if _LIBCPP_STD_VER > 14
+  testAssocLegendre<float>(1e-6f, 1e-6f);
+  testAssocLegendre<double>(1e-9, 1e-9);
+  testAssocLegendre<long double>(1e-12, 1e-12);
+#endif
+  return 0;
+}

libcxx/test/libcxx/experimental/numerics/c.math/hermite.pass.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
+//
+//===----------------------------------------------------------------------===//
+
+// <experimental/cmath>
+
+#include <cassert>
+#include <experimental/cmath>
+#include <iostream>
+#include <limits>
+#include <vector>
+
+#if _LIBCPP_STD_VER > 14
+
+template <class T> void testHermiteNaNPropagation() {
+  const unsigned MaxN = 127;
+  const T x = std::numeric_limits<T>::quiet_NaN();
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    assert(std::isnan(std::experimental::hermite(n, x)));
+  }
+}
+
+template <class T> void testHermiteNotNaN(const T x) {
+  assert(!std::isnan(x));
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    assert(!std::isnan(std::experimental::hermite(n, x)));
+  }
+}
+
+template <class T>
+void testHermiteAnalytic(const T x, const T AbsTolerance,
+                         const T RelTolerance) {
+  assert(!std::isnan(x));
+  const auto compareFloatingPoint =
+      [AbsTolerance, RelTolerance](const T Result, const T Expected) {
+        if (std::isinf(Expected) && std::isinf(Result))
+          return true;
+
+        if (std::isnan(Expected) || std::isnan(Result))
+          return false;
+
+        const T Tolerance = AbsTolerance + std::abs(Expected) * RelTolerance;
+        return std::abs(Result - Expected) < Tolerance;
+      };
+
+  const auto h0 = [](T) { return T(1); };
+  const auto h1 = [](T x) { return T(2) * x; };
+  const auto h2 = [](T x) { return T(4) * x * x - T(2); };
+  const auto h3 = [](T x) { return x * (T(8) * x * x - T(12)); };
+  const auto h4 = [](T x) {
+    return (T(16) * x * x * x * x - T(48) * x * x + T(12));
+  };
+  const auto h5 = [](T x) {
+    return x * (T(32) * x * x * x * x - T(160) * x * x + T(120));
+  };
+
+  assert(compareFloatingPoint(std::experimental::hermite(0, x), h0(x)));
+  assert(compareFloatingPoint(std::experimental::hermite(1, x), h1(x)));
+  assert(compareFloatingPoint(std::experimental::hermite(2, x), h2(x)));
+  assert(compareFloatingPoint(std::experimental::hermite(3, x), h3(x)));
+  assert(compareFloatingPoint(std::experimental::hermite(4, x), h4(x)));
+  assert(compareFloatingPoint(std::experimental::hermite(5, x), h5(x)));
+}
+
+/// \details This method checks if the following recurrence relation holds:
+/// \f[
+/// H_{n+1}(x) = 2x H_{n}(x) - 2n H_{n-1}(x)
+/// \f]
+template <class T>
+void testRecurrenceRelation(T x, T RelTolerance, T AbsTolerance) {
+  const unsigned MaxN = 127;
+  for (unsigned n = 1; n < MaxN; ++n) {
+    const T HermiteNext = std::experimental::hermite(n + 1, x);
+    const T HermiteNextRecurrence =
+        T(2) * x * std::experimental::hermite(n, x) -
+        T(2) * T(n) * std::experimental::hermite(n - 1, x);
+    const T Tolerance = AbsTolerance + std::abs(HermiteNext) * RelTolerance;
+    const T Error = std::abs(HermiteNextRecurrence - HermiteNext);
+
+    if (std::isinf(HermiteNext))
+      break;
+    assert(Error < Tolerance);
+  }
+}
+
+template <class T> void testRecurrenceRelation(T RelTolerance, T AbsTolerance) {
+  const T Samples[] = {T(-1.0), T(-0.5), T(-0.1), T(0.0),
+                       T(0.1),  T(0.5),  T(1.0)};
+  for (T x : Samples)
+    testRecurrenceRelation(x, RelTolerance, AbsTolerance);
+}
+
+/// \note Roots are taken from
+/// Salzer, Herbert E., Ruth Zucker, and Ruth Capuano.
+/// Table of the zeros and weight factors of the first twenty Hermite
+/// polynomials. US Government Printing Office, 1952.
+template <class T> std::vector<T> getHermiteRoots(unsigned n) {
+  if (n == 0u)
+    return {};
+  if (n == 1u)
+    return {T(0)};
+  if (n == 2u)
+    return {T(0.707106781186548)};
+  if (n == 3u)
+    return {T(0),
+            T(1.224744871391589)};
+  if (n == 4u)
+    return {T(0.524647623275290),
+            T(1.650680123885785)};
+  if (n == 5u)
+    return {T(0), T(0.958572464613819),
+            T(2.020182870456086)};
+  if (n == 6u)
+    return {T(0.436077411927617),
+            T(1.335849074013697),
+            T(2.350604973674492)};
+  if (n == 7u)
+    return {T(0),
+            T(0.816287882858965),
+            T(1.673551628767471),
+            T(2.651961356835233)};
+  if (n == 8u)
+    return {T(0.381186990207322),
+            T(1.157193712446780),
+            T(1.981656756695843),
+            T(2.930637420257244)};
+  if (n == 9u)
+    return {T(0),
+            T(0.723551018752838),
+            T(1.468553289216668),
+            T(2.266580584531843),
+            T(3.190993201781528)};
+  if (n == 10u)
+    return {T(0.342901327223705),
+            T(1.036610829789514),
+            T(1.756683649299882),
+            T(2.532731674232790),
+            T(3.436159118837738)};
+  if (n == 11u)
+    return {T(0),
+            T(0.65680956682100),
+            T(1.326557084494933),
+            T(2.025948015825755),
+            T(2.783290099781652),
+            T(3.668470846559583)};
+
+  if (n == 12u)
+    return {T(0.314240376254359),
+            T(0.947788391240164),
+            T(1.597682635152605),
+            T(2.279507080501060),
+            T(3.020637025120890),
+            T(3.889724897869782)};
+
+  if (n == 13u)
+    return {T(0),
+            T(0.605763879171060),
+            T(1.220055036590748),
+            T(1.853107651601512),
+            T(2.519735685678238),
+            T(3.246608978372410),
+            T(4.101337596178640)};
+
+  if (n == 14u)
+    return {T(0.29174551067256),
+            T(0.87871378732940),
+            T(1.47668273114114),
+            T(2.09518325850772),
+            T(2.74847072498540),
+            T(3.46265693360227),
+            T(4.30444857047363)};
+
+  if (n == 15u)
+    return {T(0.00000000000000),
+            T(0.56506958325558),
+            T(1.13611558521092),
+            T(1.71999257518649),
+            T(2.32573248617386),
+            T(2.96716692790560),
+            T(3.66995037340445),
+            T(4.49999070730939)};
+
+  if (n == 16u)
+    return {T(0.27348104613815),
+            T(0.82295144914466),
+            T(1.38025853919888),
+            T(1.95178799091625),
+            T(2.54620215784748),
+            T(3.17699916197996),
+            T(3.86944790486012),
+            T(4.68873893930582)};
+
+  if (n == 17u)
+    return {T(0),
+            T(0.5316330013427),
+            T(1.0676487257435),
+            T(1.6129243142212),
+            T(2.1735028266666),
+            T(2.7577629157039),
+            T(3.3789320911415),
+            T(4.0619466758755),
+            T(4.8713451936744)};
+  if (n == 18u)
+    return {T(0.2582677505191),
+            T(0.7766829192674),
+            T(1.3009208583896),
+            T(1.8355316042616),
+            T(2.3862990891667),
+            T(2.9613775055316),
+            T(3.5737690684863),
+            T(4.2481178735681),
+            T(5.0483640088745)};
+  if (n == 19u)
+    return {T(0),
+            T(0.5035201634239),
+            T(1.0103683871343),
+            T(1.5241706193935),
+            T(2.0492317098506),
+            T(2.5911337897945),
+            T(3.1578488183476),
+            T(3.7621873519640),
+            T(4.4285328066038),
+            T(5.2202716905375)};
+  if (n == 20u)
+    return {T(0.2453407083009),
+            T(0.7374737285454),
+            T(1.2340762153953),
+            T(1.7385377121166),
+            T(2.2549740020893),
+            T(2.7888060584281),
+            T(3.347854567332),
+            T(3.9447640401156),
+            T(4.6036824495507),
+            T(5.3874808900112)};
+
+  return {};
+}
+
+/// \param [in] Tolerance of the root. This value must be smaller than
+/// the smallest difference between adjacent roots in the given range
+/// with n <= 20.
+template <class T> void testHermiteRoots(T Tolerance) {
+  for (unsigned n = 0; n <= 20u; ++n) {
+    const auto Roots = getHermiteRoots<T>(n);
+    for (T x : Roots) {
+      // the roots are symmetric: if x is a root, so is -x
+      if (x > T(0))
+        assert(std::signbit(std::experimental::hermite(n, -x + Tolerance)) !=
+               std::signbit(std::experimental::hermite(n, -x - Tolerance)));
+      assert(std::signbit(std::experimental::hermite(n, x + Tolerance)) !=
+             std::signbit(std::experimental::hermite(n, x - Tolerance)));
+    }
+  }
+}
+
+template <class T>
+void testHermite(const T AbsTolerance, const T RelTolerance) {
+  testHermiteNaNPropagation<T>();
+  const T Samples[] = {T(-1.0), T(-0.5), T(-0.1), T(0.0),
+                       T(0.1),  T(0.5),  T(1.0)};
+
+  for (T x : Samples) {
+    testHermiteNotNaN(x);
+    testHermiteAnalytic(x, AbsTolerance, RelTolerance);
+  }
+}
+
+#endif
+
+int main(int, char **) {
+#if _LIBCPP_STD_VER > 14
+  testHermite<float>(1e-6f, 1e-6f);
+  testHermite<double>(1e-9, 1e-9);
+  testHermite<long double>(1e-12l, 1e-12l);
+
+  testRecurrenceRelation<float>(1e-6f, 1e-6f);
+  testRecurrenceRelation<double>(1e-9, 1e-9);
+  testRecurrenceRelation<long double>(1e-12l, 1e-12l);
+
+  testHermiteRoots<float>(1e-6f);
+  testHermiteRoots<double>(1e-9);
+  testHermiteRoots<long double>(1e-10l);
+#endif
+  return 0;
+}

libcxx/test/libcxx/experimental/numerics/c.math/laguerre.pass.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
+//
+//===----------------------------------------------------------------------===//
+
+// <experimental/cmath>
+
+#include <cassert>
+#include <experimental/cmath>
+#include <limits>
+
+#if _LIBCPP_STD_VER > 14
+
+template <class T> void testLaguerreNaNPropagation() {
+  const unsigned MaxN = 127;
+  const T x = std::numeric_limits<T>::quiet_NaN();
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    assert(std::isnan(std::experimental::laguerre(n, x)));
+  }
+}
+
+template <class T> void testLaguerreNotNaN(const T x) {
+  assert(!std::isnan(x));
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    assert(!std::isnan(std::experimental::laguerre(n, x)));
+  }
+}
+
+template <class T> void testLaguerreThrows(const T x) {
+#ifndef _LIBCPP_NO_EXCEPTIONS
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    bool Throws = false;
+    try {
+      std::experimental::laguerre(n, x);
+    } catch (const std::domain_error &) {
+      Throws = true;
+    }
+    assert(Throws);
+  }
+#endif // _LIBCPP_NO_EXCEPTIONS
+}
+
+template <class T>
+void testLaguerreAnalytic(const T x, const T AbsTolerance,
+                          const T RelTolerance) {
+  assert(!std::isnan(x));
+  const auto compareFloatingPoint =
+      [AbsTolerance, RelTolerance](const T Result, const T ExpectedResult) {
+        if (std::isinf(ExpectedResult) && std::isinf(Result))
+          return true;
+
+        if (std::isnan(ExpectedResult) || std::isnan(Result))
+          return false;
+
+        const T Tolerance =
+            AbsTolerance + std::abs(ExpectedResult) * RelTolerance;
+        return std::abs(Result - ExpectedResult) < Tolerance;
+      };
+
+  const auto l0 = [](T) { return T(1); };
+  const auto l1 = [](T x) { return -x + 1; };
+  const auto l2 = [](T x) { return (x * x - T(4) * x + T(2)) / T(2); };
+  const auto l3 = [](T x) {
+    return (-x * x * x + T(9) * x * x - T(18) * x + T(6)) / T(6);
+  };
+  const auto l4 = [](T x) {
+    return (x * x * x * x - T(16) * x * x * x + T(72) * x * x - T(96) * x +
+            T(24)) /
+           T(24);
+  };
+  const auto l5 = [](T x) {
+    return (-x * x * x * x * x + T(25) * x * x * x * x - T(200) * x * x * x +
+            T(600) * x * x - T(600) * x + T(120)) /
+           T(120);
+  };
+  const auto l6 = [](T x) {
+    return (x * x * x * x * x * x - T(36) * x * x * x * x * x +
+            T(450) * x * x * x * x - T(2400) * x * x * x + T(5400) * x * x -
+            T(4320) * x + T(720)) /
+           T(720);
+  };
+
+  assert(compareFloatingPoint(std::experimental::laguerre(0, x), l0(x)));
+  assert(compareFloatingPoint(std::experimental::laguerre(1, x), l1(x)));
+  assert(compareFloatingPoint(std::experimental::laguerre(2, x), l2(x)));
+  assert(compareFloatingPoint(std::experimental::laguerre(3, x), l3(x)));
+  assert(compareFloatingPoint(std::experimental::laguerre(4, x), l4(x)));
+  assert(compareFloatingPoint(std::experimental::laguerre(5, x), l5(x)));
+  assert(compareFloatingPoint(std::experimental::laguerre(6, x), l6(x)));
+}
+
+template <class T>
+void testLaguerre(const T AbsTolerance, const T RelTolerance) {
+  testLaguerreNaNPropagation<T>();
+  testLaguerreThrows<T>(T(-5));
+
+  const T Samples[] = {T(0.0), T(0.1), T(0.5), T(1.0), T(10.0)};
+
+  for (T x : Samples) {
+    testLaguerreNotNaN(x);
+    testLaguerreAnalytic(x, AbsTolerance, RelTolerance);
+  }
+}
+
+#endif
+
+int main(int, char **) {
+#if _LIBCPP_STD_VER > 14
+  testLaguerre<float>(1e-6f, 1e-6f);
+  testLaguerre<double>(1e-9, 1e-9);
+  testLaguerre<long double>(1e-12, 1e-12);
+#endif
+  return 0;
+}

libcxx/test/libcxx/experimental/numerics/c.math/legendre.pass.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
+//
+//===----------------------------------------------------------------------===//
+
+// <experimental/cmath>
+
+#include <cassert>
+#include <experimental/cmath>
+#include <limits>
+
+#if _LIBCPP_STD_VER > 14
+
+template <class T> void testLegendreNaNPropagation() {
+  const unsigned MaxN = 127;
+  const T x = std::numeric_limits<T>::quiet_NaN();
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    assert(std::isnan(std::experimental::legendre(n, x)));
+  }
+}
+
+template <class T> void testLegendreNotNaN(const T x) {
+  assert(!std::isnan(x));
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    assert(!std::isnan(std::experimental::legendre(n, x)));
+  }
+}
+
+template <class T> void testLegendreThrows(const T x) {
+#ifndef _LIBCPP_NO_EXCEPTIONS
+  const unsigned MaxN = 127;
+  for (unsigned n = 0; n <= MaxN; ++n) {
+    bool Throws = false;
+    try {
+      std::experimental::legendre(n, x);
+    } catch (const std::domain_error &) {
+      Throws = true;
+    }
+    assert(Throws);
+  }
+#endif // _LIBCPP_NO_EXCEPTIONS
+}
+
+template <class T>
+void testLegendreAnalytic(const T x, const T AbsTolerance,
+                          const T RelTolerance) {
+  assert(!std::isnan(x));
+  const auto compareFloatingPoint =
+      [AbsTolerance, RelTolerance](const T Result, const T ExpectedResult) {
+        if (std::isinf(ExpectedResult) && std::isinf(Result))
+          return true;
+
+        if (std::isnan(ExpectedResult) || std::isnan(Result))
+          return false;
+
+        const T Tolerance =
+            AbsTolerance + std::abs(ExpectedResult) * RelTolerance;
+        return std::abs(Result - ExpectedResult) < Tolerance;
+      };
+
+  const auto l0 = [](T) { return T(1); };
+  const auto l1 = [](T x) { return x; };
+  const auto l2 = [](T x) { return (T(3) * x * x - T(1)) / T(2); };
+  const auto l3 = [](T x) { return (T(5) * x * x - T(3)) * x / T(2); };
+  const auto l4 = [](T x) {
+    return (T(35) * x * x * x * x - T(30) * x * x + T(3)) / T(8);
+  };
+  const auto l5 = [](T x) {
+    return (T(63) * x * x * x * x - T(70) * x * x + T(15)) * x / T(8);
+  };
+  const auto l6 = [](T x) {
+    const T x2 = x * x;
+    return (T(231) * x2 * x2 * x2 - T(315) * x2 * x2 + T(105) * x2 - T(5)) /
+           T(16);
+  };
+
+  assert(compareFloatingPoint(std::experimental::legendre(0, x), l0(x)));
+  assert(compareFloatingPoint(std::experimental::legendre(1, x), l1(x)));
+  assert(compareFloatingPoint(std::experimental::legendre(2, x), l2(x)));
+  assert(compareFloatingPoint(std::experimental::legendre(3, x), l3(x)));
+  assert(compareFloatingPoint(std::experimental::legendre(4, x), l4(x)));
+  assert(compareFloatingPoint(std::experimental::legendre(5, x), l5(x)));
+  assert(compareFloatingPoint(std::experimental::legendre(6, x), l6(x)));
+}
+
+template <class T>
+void testLegendre(const T AbsTolerance, const T RelTolerance) {
+  testLegendreNaNPropagation<T>();
+  testLegendreThrows<T>(T(-5));
+  testLegendreThrows<T>(T(5));
+
+  const T Samples[] = {T(-1.0), T(-0.5), T(-0.1), T(0.0),
+                       T(0.1),  T(0.5),  T(1.0)};
+
+  for (T x : Samples) {
+    testLegendreNotNaN(x);
+    testLegendreAnalytic(x, AbsTolerance, RelTolerance);
+  }
+}
+
+#endif
+
+int main(int, char **) {
+#if _LIBCPP_STD_VER > 14
+  testLegendre<float>(1e-6f, 1e-6f);
+  testLegendre<double>(1e-9, 1e-9);
+  testLegendre<long double>(1e-12, 1e-12);
+#endif
+  return 0;
+}