Added std::assoc_legendre and std::legendre [sf.cmath]

[thebrandre]

Summary

Added std::assoc_legendre and std::legendre from C++17, 29.9.5 Mathematical special functions [sf.cmath].

Resubmitted parts of Added implementations of five Mathematical Special Functions added in C++17
It is tested with a table of values generated by gcc-8.3.1 and boost-1.66.0 using boost::multiprecision::cpp_bin_float_oct (256bit float).

Notes:

• The float version uses double internally because assoc_legendre will easily overflow for float. But it can be shown that assoc_legendre is bound by a value less than std::numeric_limits<double>::max() for the whole range required by the standard, x in [-1,1], l < 128, and m < 128. (There is a unit test for that and references to the literature in the code.)

So that seemed to be kind of a natural choice to me. Analog to Added std::assoc_laguerre and std::laguerre and boost, the double version is also promoted long double for internal calculations.

• The implementation is inline for the same reason as D60869.

Related revisions:

• Added std::assoc_laguerre and std::laguerre
• Added implementation of expint

[mclow.lists]

Thanks for doing all this. I'm not really ignoring these, it just seems like it. :-(

[thebrandre]

No worries. :-)
I know that this is time-consuming and I am well aware that is not a high-priority feature. And now I had time to reconsider and find my own mistakes. ;-)
The integration should be now much cleaner than my first draft.

Btw .... I should also point out that this draft implicitly checks the flag _LIBCPP_NO_EXCEPTIONS (via __throw_domain_error) in the src folder. I keep wondering which of these #defines may actually be used in cpp files. Should this be avoided so there's no reason to build different versions of the library?

[thebrandre]

To avoid AVX-SSE transition penalties (and for simplicity) the implementation is now inline.
See Added std::assoc_laguerre and std::laguerre.

I am quite certain by now that the benefits of the inline implementation outweigh the drawbacks ... but I am open for discussion/advice.

[thebrandre]

I recently became aware of the fact that gcc and msvc don't aggree on "Exceptions or Error Codes?" (Section III.D in N1542)

N1542 is in favor of signaling domain errors via

#include <cerrno>
errno = EDOM;

for consistency with other functions from cmath such as the log function and for compatibility with C.

VS2019's implementation does exactly that whereas gcc throws an std::domain_error. (As an aside: the exceptions of legendre got removed in gcc-9 but laguerre still throws for x<0 .)

The C++17 standard and N3060 simply refer to it as "domain error". Unfortunately, this term is ambiguous because both EDOM and std::domain_error are known by this name.

As N1542 suggests to use errno for perfectly valid reasons, I am inclined to change my patches accordingly even though it means that we'd diverge from libstdc++ in that matter.
Maybe/hopefully someone can help out with any (official) statements that give clearer directions?

Interestingly, boost is in-between: when you use boost::math::policies::errno_on_error instead of the default policy, you'll get close to MSVC's behavior. The default policy throws std::domain_error.

I'd also like to point out that the situation is even worse than that. For example, log(-1) not only sets errno but also raises a floating-point exception as if

std::feraiseexcept(FE_INVALID);

were called. This is also referred to as a "domain error".

We could mimic the behavior of the log function by introducing something like the following:

#include <cerrno>
#include <cmath>
#include <cfenv>

void __raise_math_domain_error() _NOEXCEPT {
#if math_errhandling & MATH_ERREXCEPT == MATH_ERREXCEPT
  std::feraiseexcept(FE_INVALID);
#endif
#if math_errhandling & MATH_ERRNO == MATH_ERRNO
  errno = EDOM;
#endif
}

Let's call that one "errno + floating-point exception".

However, neither boost nor msvc seem to go beyond "errno".
So when it comes down to the question "std::domain_error, errno, or 'errno + floating-point exception'?", I personally think errno is the way to go ...
but I would be sooooo happy for some feedback on that. ;-)

[thebrandre]

After a vivid discussion with myself, I went with "errno + floating-point exception".
In fact, the C11 Standard has a clear definition of a domain error in 7.12.1.
I added cmath.cpp with a definition of a function that sets the domain error as specified (similar to the comment above, just minus the bugs). It is not inline to avoid an include of cfenv and cerrno.

Note: the tests on invalid arguments only check for errno because clang doesn't support the STDC FENV_ACCESS pragma.

Also some minor changes:

• Moved the tests for consistency with the latest changes in cmath (lerp).
• Simplification of the tests: now using only one set of reference values generated by boost with multiprecision floats (256bit rounded to 128bit). I am pretty optimistic that every digit is correct, so that the comparison to other implementations no longer serves any purpose.
• Rebased against current master branch.

libcxx/include/__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__CMATH
+#define _LIBCPP__CMATH
+
+#include <__config>
+
+#if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
+#pragma GCC system_header
+#endif
+
+_LIBCPP_BEGIN_NAMESPACE_STD
+
+// __libcpp_assoc_legendre
+
+_LIBCPP_FUNC_VIS
+float __libcpp_assoc_legendref(unsigned int __lcpp_l, unsigned int __lcpp_m,
+                               float __lcpp_x);
+
+_LIBCPP_FUNC_VIS
+double __libcpp_assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m,
+                               double __lcpp_x);
+
+_LIBCPP_FUNC_VIS
+long double __libcpp_assoc_legendrel(unsigned int __lcpp_l,
+                                     unsigned int __lcpp_m,
+                                     long double __lcpp_x);
+
+
+// __libcpp_legendre
+
+_LIBCPP_FUNC_VIS
+float __libcpp_legendref(unsigned int __lcpp_n, float __lcpp_x);
+
+_LIBCPP_FUNC_VIS
+double __libcpp_legendre(unsigned int __lcpp_n, double __lcpp_x);
+
+_LIBCPP_FUNC_VIS
+long double __libcpp_legendrel(unsigned int __lcpp_n, long double __lcpp_x);
+
+_LIBCPP_END_NAMESPACE_STD
+
+#endif  // _LIBCPP__CMATH

libcxx/include/cmath

 float          truncf(float x);
 long double    truncl(long double x);
 
+// C++17 Mathematical special functions [sf.cmath]
+
+floating_point assoc_legendre (unsigned l, unsigned m, arithmetic x);
+float          assoc_legendref(unsigned l, unsigned m, float x);
+long double    assoc_legendrel(unsigned l, unsigned m, long double x);
+
+floating_point legendre (unsigned n, arithmetic x);
+float          legendref(unsigned n, float x);
+long double    legendrel(unsigned n, long double x);
+
 }  // std
 
 */

 #include <__config>
 #include <math.h>
 #include <version>
+#include <__cmath>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
 #pragma GCC system_header

 using ::truncl;
 
 #if _LIBCPP_STD_VER > 14
+
+// assoc_legendre
+
+inline _LIBCPP_INLINE_VISIBILITY
+float
+assoc_legendref(unsigned int __lcpp_l, unsigned int __lcpp_m, float __lcpp_x) {
+    return _VSTD::__libcpp_assoc_legendref(__lcpp_l, __lcpp_m, __lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+double
+assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, double __lcpp_x) {
+    return _VSTD::__libcpp_assoc_legendre(__lcpp_l, __lcpp_m, __lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+long double
+assoc_legendrel(unsigned int __lcpp_l, unsigned int __lcpp_m, long double __lcpp_x) {
+    return _VSTD::__libcpp_assoc_legendrel(__lcpp_l, __lcpp_m, __lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+float
+assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, float __lcpp_x) {
+    return _VSTD::assoc_legendref(__lcpp_l, __lcpp_m, __lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+long double
+assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, long double __lcpp_x) {
+    return _VSTD::assoc_legendrel(__lcpp_l, __lcpp_m, __lcpp_x);
+}
+
+template <class _A1>
+inline _LIBCPP_INLINE_VISIBILITY
+typename std::enable_if<std::is_integral<_A1>::value, double>::type
+assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m, _A1 __lcpp_x) {
+  return _VSTD::assoc_legendre(__lcpp_l, __lcpp_m, (double)__lcpp_x);
+}
+
+// legendre
+
+inline _LIBCPP_INLINE_VISIBILITY
+float
+legendref(unsigned int __lcpp_n, float __lcpp_x) {
+    return _VSTD::__libcpp_legendref(__lcpp_n, __lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+double
+legendre(unsigned int __lcpp_n, double __lcpp_x) {
+    return _VSTD::__libcpp_legendre(__lcpp_n, __lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+long double
+legendrel(unsigned int __lcpp_n, long double __lcpp_x) {
+    return _VSTD::__libcpp_legendrel(__lcpp_n, __lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+float
+legendre(unsigned int __lcpp_n, float __lcpp_x) {
+    return _VSTD::legendref(__lcpp_n, __lcpp_x);
+}
+
+inline _LIBCPP_INLINE_VISIBILITY
+long double
+legendre(unsigned int __lcpp_n, long double __lcpp_x) {
+    return _VSTD::legendrel(__lcpp_n, __lcpp_x);
+}
+
+template <class _A1>
+inline _LIBCPP_INLINE_VISIBILITY
+typename std::enable_if<std::is_integral<_A1>::value, double>::type
+legendre(unsigned int __lcpp_n, _A1 __lcpp_x) {
+  return _VSTD::legendre(__lcpp_n, (double)__lcpp_x);
+}
+
+// hypot
+
 inline _LIBCPP_INLINE_VISIBILITY float       hypot(       float x,       float y,       float z ) { return sqrt(x*x + y*y + z*z); }
 inline _LIBCPP_INLINE_VISIBILITY double      hypot(      double x,      double y,      double z ) { return sqrt(x*x + y*y + z*z); }
 inline _LIBCPP_INLINE_VISIBILITY long double hypot( long double x, long double y, long double z ) { return sqrt(x*x + y*y + z*z); }

libcxx/lib/CMakeLists.txt

   endif()
 endif()
 
+list(APPEND LIBCXX_SOURCES ../src/cmath/legendre.cpp)
+
 # Add all the headers to the project for IDEs.
 if (LIBCXX_CONFIGURE_IDE)
   file(GLOB_RECURSE LIBCXX_HEADERS ${CMAKE_CURRENT_SOURCE_DIR}/../include/*)

libcxx/src/cmath/legendre.cpp

+//===----------------------- legendre.cpp -----------------------*- 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
+/// and std::assoc_legendre  (as of C++17)
+/// std::assoc_legendre, std::assoc_legendref, and std::assoc_legendrel
+/// correspond to __libcpp_assoc_legendre, __libcpp_assoc_legendref, and
+/// __libcpp_assoc_legendrel.
+/// std::legendre, std::legendref, and std::legendrel correspond to
+/// __libcpp_legendre, __libcpp_legendref,  and __libcpp_legendrel.
+///
+//===----------------------------------------------------------------------===//
+
+#include <cmath>
+#include <limits>
+#include <stdexcept>
+
+_LIBCPP_BEGIN_NAMESPACE_STD
+
+/// \brief Evaluates the Legendre polynomial \f$ P_{n}(x) \f$.
+///
+/// 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]
+/// where \f$ P_{0}(x) = 1 \f$ and \f$ P_{1}(x) = x \f$.
+/// Press, William H., et al. Numerical recipes 3rd edition: The art of
+/// scientific computing. Cambridge university press, 2007, p. 182.
+/// \tparam T one of float, double, or long double
+/// \param [in] __n corresponds to the degree \f$ n \f$.
+/// \param [in] __x corresponds to the argument \f$ x \f$.
+/// \return the Legendre polynomial \f$ P_{n}(x) \f$
+template <class _Real>
+inline _LIBCPP_INLINE_VISIBILITY _Real
+__libcpp_legendre_recurrence(unsigned __n, _Real __x) _NOEXCEPT {
+  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;
+}
+
+/// \brief Evaluates the Legendre polynomial \f$ P_{n}(x) \f$.
+///
+/// \tparam T one of float, double, or long double
+/// \param [in] __n corresponds to the degree \f$ n \f$.
+/// \param [in] __x corresponds to the argument \f$ x \f$.
+/// \return the Legendre polynomial \f$ P_{n}(x) \f$
+/// \throw std::domain_error if __x is not inside the interval [-1,1].
+template <class _Real>
+inline _LIBCPP_INLINE_VISIBILITY _Real __libcpp_legendre_impl(unsigned __n,
+                                                              _Real __x) {
+  if (std::isnan(__x))
+    return std::numeric_limits<_Real>::quiet_NaN();
+
+  if (std::abs(__x) > _Real(1))
+    __throw_domain_error("Argument of legendre function is out of range");
+
+  return __libcpp_legendre_recurrence(__n, __x);
+}
+
+/// \brief Evaluates the associated Legendre polynomial \f$ P_{l}^{m} \f$.
+///
+/// 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]
+///
+/// \note The starting point of the recursion grows exponentially with
+/// \f$ m \f$. For large \f$ m \f$, 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_legendre.
+/// \param [in] __l corresponds to \f$ l \f$.
+/// \param [in] __m corresponds to \f$ m \f$.
+/// \param [in] __x corresponds to the argument \f$ x \f$.
+/// \return \f$ P_{l}^{m}(x) \f$
+template <class _Real>
+inline _LIBCPP_INLINE_VISIBILITY _Real __libcpp_assoc_legendre_recurrence(
+    unsigned __l, unsigned __m, _Real __x) _NOEXCEPT {
+  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));
+  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;
+}
+
+/// \brief Evaluates the associated Legendre polynomial \f$ P_{l}^{m} \f$.
+///
+/// \param [in] __l corresponds to \f$ l \f$.
+/// \param [in] __m corresponds to \f$ m \f$.
+/// \param [in] __x corresponds to the argument \f$ x \f$.
+/// \return \f$ P_{l}^{m}(x) \f$
+/// \throw std::domain_error if __x is not inside the interval [-1,1].
+template <class _Real>
+inline _LIBCPP_INLINE_VISIBILITY _Real
+__libcpp_assoc_legendre_impl(unsigned __l, unsigned __m, _Real __x) {
+  if (std::isnan(__x))
+    return std::numeric_limits<_Real>::quiet_NaN();
+
+  if (std::abs(__x) > _Real(1))
+    __throw_domain_error("Argument of assoc_legendre function is out of range");
+
+  return __libcpp_assoc_legendre_recurrence(__l, __m, __x);
+}
+
+/// \brief Evaluates the associated Legendre polynomial \f$ P_{l}^{m} \f$.
+///
+/// We use the double version internally because assoc_legendre
+/// is prone to overflow. It can be shown that assoc_legendre is bound
+/// by a value less than std::numeric_limits<double>::max() for the whole
+/// range required by the standard, __lcpp_x in [-1,1}, __lcpp_l < 128,
+/// and __lcpp_m < 128.
+/// This is elaborated in more detail in the unit tests
+/// See Eq. (22) on p. 183 of Lohoefer, G. (1998).
+/// Inequalities for the associated Legendre functions.
+float __libcpp_assoc_legendref(unsigned int __lcpp_l, unsigned int __lcpp_m,
+                               float __lcpp_x) {
+  return static_cast<float>(__libcpp_assoc_legendre_impl<double>(__lcpp_l, __lcpp_m, static_cast<double>(__lcpp_x)));
+}
+
+double __libcpp_assoc_legendre(unsigned int __lcpp_l, unsigned int __lcpp_m,
+                               double __lcpp_x) {
+  return __libcpp_assoc_legendre_impl<double>(__lcpp_l, __lcpp_m, __lcpp_x);
+}
+
+long double __libcpp_assoc_legendrel(unsigned int __lcpp_l,
+                                     unsigned int __lcpp_m,
+                                     long double __lcpp_x) {
+  return __libcpp_assoc_legendre_impl<long double>(__lcpp_l, __lcpp_m,
+                                                   __lcpp_x);
+}
+
+/// \brief Evaluates the Legendre polynomial \f$ P_{n}(x) \f$.
+///
+/// In contrast to assoc_legendre, legendre won't overflow in common scenarios.
+/// Nevertheless, we also use the double implementation to make sure that
+/// legendre(n,x) returns assoc_legendre(n,0,x).
+float __libcpp_legendref(unsigned int __lcpp_n, float __lcpp_x) {
+  return static_cast<float>(__libcpp_legendre_impl<double>(__lcpp_n, static_cast<double>(__lcpp_x)));
+}
+
+double __libcpp_legendre(unsigned int __lcpp_n, double __lcpp_x) {
+  return __libcpp_legendre_impl<double>(__lcpp_n, __lcpp_x);
+}
+
+long double __libcpp_legendrel(unsigned int __lcpp_n, long double __lcpp_x) {
+  return __libcpp_legendre_impl<long double>(__lcpp_n, __lcpp_x);
+}
+
+_LIBCPP_END_NAMESPACE_STD

libcxx/test/std/numerics/legendre/assoc_legendre_basic.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
+//
+//===----------------------------------------------------------------------===//
+
+// UNSUPPORTED: c++98, c++03, c++11, c++14
+
+// <cmath>
+
+// floating_point assoc_legendre (unsigned l, unsigned m, arithmetic x);
+// float          assoc_legendref(unsigned l, unsigned m, float x);
+// long double    assoc_legendrel(unsigned l, unsigned m, long double x);
+
+#include <cmath>
+#include <cassert>
+#include <limits>
+
+/// \brief std::assoc_legendre(n,m,x) must silently return NaN if x is NaN.
+///
+/// \tparam T one of float, double, or long double
+template <class T>
+void testNaNPropagation() {
+  for (unsigned l = 0; l < 128; ++l) {
+    for (unsigned m = 0; m < 128; ++m) {
+      const T ExpectNaN =
+          std::assoc_legendre(l, m, std::numeric_limits<T>::quiet_NaN());
+      assert(std::isnan(ExpectNaN));
+    }
+  }
+}
+
+/// \brief Tests if the associated Legendre polynomial behaves
+/// correctly under the transformation x -> -x
+///
+/// The expected property is
+/// \f$ P_{l} ^{m} (-x) = (-1)^{l + m} P_{l}^{m}(x) \f$.
+/// For odd \f$ l + m \f$, this test checks for a root at \f$ x = 0\f$.
+/// \tparam T one of float, double, or long double
+/// \param [in] l first argument of assoc_legendre
+/// \param [in] m second argument of assoc_legendre
+/// \param [in] x argument of the Legendre polynomial
+/// \param [in] RelTolerance relative tolerance for the value of legendre
+/// \param [in] AbsTolerance absolute tolerance for the value of legendre
+template <class T>
+void testParity(unsigned l, unsigned m, T x, T RelTolerance, T AbsTolerance) {
+  const T Parity = (l + m) % 2 == 0 ? T(1) : T(-1);
+  const T ExpectedValue = Parity * std::assoc_legendre(l, m, x);
+  const T Tolerance = RelTolerance * std::abs(ExpectedValue) + AbsTolerance;
+  const T Error = std::abs(std::assoc_legendre(l, m, -x) - ExpectedValue);
+  assert(Error <= Tolerance);
+}
+
+template <class T>
+void testParity(T RelTolerance, T AbsTolerance) {
+  const T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)};
+  for (unsigned l = 0; l < 128; ++l) {
+    for (unsigned m = 0; m < 128; ++m) {
+      for (T x : Samples)
+        testParity<T>(l, m, x, RelTolerance, AbsTolerance);
+    }
+  }
+}
+
+/// \brief Tests if the associated Legendre polynomial adheres
+/// to an upper bound from the literature.
+///
+/// For $ 1 \leq m \leq l $, we have:
+/// \f[
+/// \frac{1}{\sqrt{2.22(m+1)}} <
+/// \sqrt{\frac{(l-m)!}{(l+m)!}} \cdot
+/// \max\limits_{x\in\left[-1,1\right]}\left| P_{l}^{m} (x) \right|
+/// < \sqrt[4]{\frac{2^5}{\pi^3 m}}
+/// \f]
+/// See Eq. (22) on p. 183 of Lohoefer, G. (1998).
+/// Inequalities for the associated Legendre functions.
+/// Journal of Approximation Theory, 95(2), 178-193.
+/// \note For \f$ l \leq 127, the maximum upper bound is
+/// 3.44177e+250, which is well below std::numeric_limits<double>::max()
+/// of 1.79769e+308 but larger than the maximum value of float (3.40282e+38).
+/// The lower bound for the maximum is 6.80027e+249. So the float version
+/// will have to deal with overflows.
+/// \tparam T one of float, double, or long double
+/// \param [in] l first argument of assoc_legendre
+/// \param [in] m second argument of assoc_legendre
+/// \param [in] x argument of the Legendre polynomial
+/// \param [in] RelTolerance relative tolerance for the value of legendre
+/// \param [in] AbsTolerance absolute tolerance for the value of legendre
+template <class T>
+void testUpperBound(unsigned l, unsigned m, T x, T RelTolerance,
+                    T AbsTolerance) {
+  const T Pi = T(3.14159265358979323846264338327950288L);
+  const T LogUpperBound =
+      std::lgamma(T(l + m + 1u)) / T(2) - std::lgamma(T(l - m + 1u)) / T(2) -
+      std::log(m) / T(4) + T(1.25) * std::log(2) - T(0.75) * std::log(Pi);
+  const T UpperBound = std::exp(LogUpperBound);
+  if (std::isinf(UpperBound))
+    return;
+
+  const T UpperBoundWithTolerance =
+      UpperBound + UpperBound * RelTolerance + AbsTolerance;
+  assert(std::assoc_legendre(l, m, x) <= UpperBoundWithTolerance);
+}
+
+template <class T>
+void testUpperBound(T RelTolerance, T AbsTolerance) {
+  const T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)};
+  for (unsigned l = 1; l < 128; ++l) {
+    for (unsigned m = 1; m <= l; ++m) {
+      for (T x : Samples)
+        testUpperBound<T>(l, m, x, RelTolerance, AbsTolerance);
+    }
+  }
+}
+
+int main(int, char**) {
+  testNaNPropagation<float>();
+  testNaNPropagation<double>();
+  testNaNPropagation<long double>();
+
+  // We only test double here because the float version of std::assoc_legendre
+  // often returns inf for large l and m.
+  testParity<double>(1e-12, 0.0);
+  testUpperBound<double>(1e-12, 0.0);
+}

libcxx/test/std/numerics/legendre/assoc_legendre_comparisons.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
+//
+//===----------------------------------------------------------------------===//
+
+// UNSUPPORTED: c++98, c++03, c++11, c++14
+
+// <cmath>
+
+// floating_point assoc_legendre (unsigned l, unsigned m, arithmetic x);
+// float          assoc_legendref(unsigned l, unsigned m, float x);
+// long double    assoc_legendrel(unsigned l, unsigned m, long double x);
+
+#include <cmath>
+#include <map>
+#include <vector>
+#include <cassert>
+
+template <class Real, class Iterator>
+Real evaluatePolynomial(Iterator CoefficientsBegin, Iterator CoefficientsEnd,
+                        Real x) {
+  Real Result = Real(0);
+  for (Iterator Iter = CoefficientsBegin; Iter != CoefficientsEnd; ++Iter) {
+    Result = Result * x + *Iter;
+  }
+  return Result;
+}
+
+template <class Real, class Container>
+Real evaluatePolynomial(const Container& Coefficients, Real x) {
+  return evaluatePolynomial(std::begin(Coefficients), std::end(Coefficients),
+                            x);
+}
+
+template <class T, class TestFunction>
+void compareWithAnalyticExpressions(const T x, TestFunction&& ValueSink) {
+  std::vector<std::vector<std::vector<T> > > TableOfCoefficients;
+
+  // m = 1 and l = 1..4
+  TableOfCoefficients.emplace_back();
+  TableOfCoefficients.back().push_back({T(1)});
+  TableOfCoefficients.back().push_back({T(3), T(0)});
+  TableOfCoefficients.back().push_back({T(15) / T(2), T(0), -T(3) / T(2)});
+  TableOfCoefficients.back().push_back(
+      {T(35) / T(2), T(0), -T(15) / T(2), T(0)});
+
+  // m = 2 and l = 1..4 without the factor (1-x^2)^{m/2}
+  TableOfCoefficients.emplace_back();
+  TableOfCoefficients.back().emplace_back();
+  TableOfCoefficients.back().push_back({T(3)});
+  TableOfCoefficients.back().push_back({T(15), T(0)});
+  TableOfCoefficients.back().push_back({T(105) / T(2), T(0), -T(15) / T(2)});
+
+  // m = 3 and l = 1..4 without the factor (1-x^2)^{m/2}
+  TableOfCoefficients.emplace_back();
+  TableOfCoefficients.back().emplace_back();
+  TableOfCoefficients.back().emplace_back();
+  TableOfCoefficients.back().push_back({T(15)});
+  TableOfCoefficients.back().push_back({T(105), T(0)});
+
+  // m = 4 and l = 1..4 without the factor (1-x^2)^{m/2}
+  TableOfCoefficients.emplace_back();
+  TableOfCoefficients.back().emplace_back();
+  TableOfCoefficients.back().emplace_back();
+  TableOfCoefficients.back().emplace_back();
+  TableOfCoefficients.back().push_back({T(105)});
+
+  for (size_t i = 0; i < TableOfCoefficients.size(); ++i) {
+    const unsigned m = static_cast<unsigned int>(i) + 1;
+    const T Factor = std::pow((T(1) - x) * (T(1) + x), T(m) / T(2));
+    for (size_t j = 0; j < TableOfCoefficients[i].size(); ++j) {
+      const unsigned int l = static_cast<unsigned int>(j) + 1;
+      const T ExpectedValue =
+          Factor * evaluatePolynomial(TableOfCoefficients[i][j], x);
+      ValueSink(l, m, x, ExpectedValue);
+    }
+  }
+}
+
+template <class T>
+void compareWithAnalyticExpressions(const T AbsTolerance,
+                                    const T RelTolerance) {
+  const auto Samples = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)};
+  auto testFunction = [RelTolerance, AbsTolerance](unsigned l, unsigned m, T x,
+                                                   T ExpectedValue) {
+    const T Result = std::assoc_legendre(l, m, x);
+    const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance;
+    const T Error = std::abs(Result - ExpectedValue);
+    assert(Error <= Tolerance);
+  };
+  for (T x : Samples)
+    compareWithAnalyticExpressions(x, testFunction);
+}
+
+template <class T>
+void compareWithTable(
+    std::map<std::pair<unsigned, unsigned>, std::map<long double, long double> >
+        Table,
+    T RelTolerance, T AbsTolerance) {
+  for (const auto& Entry : Table) {
+    const unsigned l = Entry.first.first;
+    const unsigned m = Entry.first.second;
+    for (const auto& xy : Entry.second) {
+      const T ExpectedValue = xy.second;
+      const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance;
+      const T Error =
+          std::abs(std::assoc_legendre(l, m, xy.first) - ExpectedValue);
+      assert(Error <= Tolerance);
+    }
+  }
+}
+
+template <class T>
+void compareWithMSVC(T RelTolerance, T AbsTolerance) {
+  // Generated with msvc 191627027
+  std::map<std::pair<unsigned, unsigned>, std::map<long double, long double> >
+      Table;
+  Table[{10, 10}] = {
+      {0.0L, 6.54729075000000000e+08L}, {0.1L, 6.22640835705230832e+08L},
+      {0.2L, 5.33848212079902649e+08L}, {0.3L, 4.08571989131588876e+08L},
+      {0.4L, 2.73815518201505244e+08L}, {0.5L, 1.55370278540039063e+08L},
+      {0.6L, 7.03009991216332912e+07L}, {0.7L, 2.25898063438872844e+07L},
+      {0.8L, 3.95889634812671319e+06L}, {0.9L, 1.62117400787842256e+05L}};
+  Table[{20, 10}] = {
+      {0.0L, -1.61205295667431641e+12L}, {0.1L, 3.56041756390717590e+11L},
+      {0.2L, 1.46660926765988135e+12L},  {0.3L, -1.09242045429771387e+12L},
+      {0.4L, -8.99973487310552368e+11L}, {0.5L, 1.74450229594620776e+12L},
+      {0.6L, -4.24471915195055054e+11L}, {0.7L, -1.54836452494957959e+12L},
+      {0.8L, 2.39275793328105029e+12L},  {0.9L, 9.90017476694990845e+11L}};
+  Table[{20, 20}] = {
+      {0.0L, 3.19830986772877752e+23L}, {0.1L, 2.89249411469768665e+23L},
+      {0.2L, 2.12634078007975140e+23L}, {0.3L, 1.24547341322977873e+23L},
+      {0.4L, 5.59388325840124914e+22L}, {0.5L, 1.80108069781796049e+22L},
+      {0.6L, 3.68740022490078904e+21L}, {0.7L, 3.80734558806206054e+20L},
+      {0.8L, 1.16935276168332288e+19L}, {0.9L, 1.96090497120238360e+16L}};
+  Table[{50, 10}] = {
+      {0.0L, 1.14522382492463460e+16L},  {0.1L, 2.79436816274998850e+15L},
+      {0.2L, -9.92643944667355200e+15L}, {0.3L, -9.45138403285157200e+15L},
+      {0.4L, 8.89201701866905625e+14L},  {0.5L, 9.18003538846575400e+15L},
+      {0.6L, 1.21766907555422340e+16L},  {0.7L, 1.21535350115069900e+16L},
+      {0.8L, 9.32311375306574375e+14L},  {0.9L, -8.99466171009315500e+15L}};
+  Table[{50, 20}] = {
+      {0.0L, -7.87359356973322099e+32L}, {0.1L, 5.41513582908981306e+31L},
+      {0.2L, 7.93092660564343958e+32L},  {0.3L, -4.00540672362456500e+31L},
+      {0.4L, -8.23243252797740372e+32L}, {0.5L, -4.01698608142744523e+32L},
+      {0.6L, 4.08601287568084080e+32L},  {0.7L, 8.60225211647169535e+32L},
+      {0.8L, 1.06226766578920635e+33L},  {0.9L, 1.66309251586455094e+33L}};
+  Table[{50, 50}] = {
+      {0.0L, 2.72539213975072988e+78L}, {0.1L, 2.11986821903666169e+78L},
+      {0.2L, 9.82222379316811538e+77L}, {0.3L, 2.57907402241500578e+77L},
+      {0.4L, 3.48667245334433154e+76L}, {0.5L, 2.05097602570372398e+75L},
+      {0.6L, 3.88980964317819002e+73L}, {0.7L, 1.33254965595665771e+71L},
+      {0.8L, 2.20288344031012791e+67L}, {0.9L, 2.53669949744314049e+60L}};
+  Table[{100, 10}] = {
+      {0.0L, -8.24945971816703795e+18L}, {0.1L, 6.86278552250344038e+18L},
+      {0.2L, -2.34688103580907827e+18L}, {0.3L, -4.90893583880523878e+18L},
+      {0.4L, 8.23708346184265728e+18L},  {0.5L, 4.23871230219633664e+18L},
+      {0.6L, -9.59361116591237120e+17L}, {0.7L, 4.43550484874968269e+18L},
+      {0.8L, 1.57786735456734694e+18L},  {0.9L, 2.56433959576586016e+17L}};
+  Table[{100, 20}] = {
+      {0.0L, 7.77274398361595813e+38L},  {0.1L, -7.05032664517025152e+38L},
+      {0.2L, 4.39662513074438558e+38L},  {0.3L, 1.14744968919013248e+38L},
+      {0.4L, -7.61205860438438888e+38L}, {0.5L, 3.07504970399995784e+38L},
+      {0.6L, 8.26800055741207868e+38L},  {0.7L, 7.90222368531099936e+38L},
+      {0.8L, -8.39638951169623065e+38L}, {0.9L, 7.16043448787788025e+37L}};
+  Table[{100, 50}] = {
+      {0.0L, -1.17121450315783826e+98L}, {0.1L, 9.01193385501805844e+97L},
+      {0.2L, -2.77911495881219591e+97L}, {0.3L, -3.35206020674796850e+97L},
+      {0.4L, 7.14313850937400452e+97L},  {0.5L, -8.22526203397273908e+97L},
+      {0.6L, 5.93640459256693667e+97L},  {0.7L, 1.80893339034652566e+97L},
+      {0.8L, -1.10297979774543506e+98L}, {0.9L, 9.32312785168941904e+96L}};
+  Table[{100, 100}] = {
+      {0.0L, 6.66630867007295431e+186L}, {0.1L, 4.03315719080569865e+186L},
+      {0.2L, 8.65858791475270039e+185L}, {0.3L, 5.96973475268218189e+184L},
+      {0.4L, 1.09106273304589198e+183L}, {0.5L, 3.77527496828897112e+180L},
+      {0.6L, 1.35795105902891242e+177L}, {0.7L, 1.59365468505949483e+172L},
+      {0.8L, 4.35522360415857506e+164L}, {0.9L, 5.77517922557585923e+150L}};
+  compareWithTable(Table, RelTolerance, AbsTolerance);
+}
+
+template <class T>
+void compareWithGCC(T RelTolerance, T AbsTolerance) {
+  // Generated with gcc 8.3.1
+  std::map<std::pair<unsigned, unsigned>, std::map<long double, long double> >
+      Table;
+  Table[{10, 10}] = {{0.0L, 6.547290750000000000000e+08L},
+                     {0.1L, 6.226408357052308421698e+08L},
+                     {0.2L, 5.338482120799027200555e+08L},
+                     {0.3L, 4.085719891315887827659e+08L},
+                     {0.4L, 2.738155182015052800416e+08L},
+                     {0.5L, 1.553702785400390624272e+08L},
+                     {0.6L, 7.030099912163327995950e+07L},
+                     {0.7L, 2.258980634388728250269e+07L},
+                     {0.8L, 3.958896348126719998163e+06L},
+                     {0.9L, 1.621174007878425001223e+05L}};
+  Table[{20, 10}] = {{0.0L, -1.612052956674316406250e+12L},
+                     {0.1L, 3.560417563907174735665e+11L},
+                     {0.2L, 1.466609267659881798506e+12L},
+                     {0.3L, -1.092420454297714176655e+12L},
+                     {0.4L, -8.999734873105529630780e+11L},
+                     {0.5L, 1.744502295946207409143e+12L},
+                     {0.6L, -4.244719151950563148558e+11L},
+                     {0.7L, -1.548364524949579275370e+12L},
+                     {0.8L, 2.392757933281052241802e+12L},
+                     {0.9L, 9.900174766949915345311e+11L}};
+  Table[{20, 20}] = {{0.0L, 3.198309867728777708175e+23L},
+                     {0.1L, 2.892494114697687005594e+23L},
+                     {0.2L, 2.126340780079752328315e+23L},
+                     {0.3L, 1.245473413229778328863e+23L},
+                     {0.4L, 5.593883258401251200614e+22L},
+                     {0.5L, 1.801080697817960689050e+22L},
+                     {0.6L, 3.687400224900788126976e+21L},
+                     {0.7L, 3.807345588062060330880e+20L},
+                     {0.8L, 1.169352761683327100000e+19L},
+                     {0.9L, 1.960904971202388847266e+16L}};
+  Table[{50, 10}] = {{0.0L, 1.145223824924634340332e+16L},
+                     {0.1L, 2.794368162749986069824e+15L},
+                     {0.2L, -9.926439446673560716797e+15L},
+                     {0.3L, -9.451384032851574123047e+15L},
+                     {0.4L, 8.892017018669203554688e+14L},
+                     {0.5L, 9.180035388465756383789e+15L},
+                     {0.6L, 1.217669075554222969727e+16L},
+                     {0.7L, 1.215353501150698340625e+16L},
+                     {0.8L, 9.323113753066150476685e+14L},
+                     {0.9L, -8.994661710093201598633e+15L}};
+  Table[{50, 20}] = {{0.0L, -7.873593569733218616572e+32L},
+                     {0.1L, 5.415135829089829040302e+31L},
+                     {0.2L, 7.930926605643441888714e+32L},
+                     {0.3L, -4.005406723624505754186e+31L},
+                     {0.4L, -8.232432527977407811105e+32L},
+                     {0.5L, -4.016986081427444523960e+32L},
+                     {0.6L, 4.086012875680830742676e+32L},
+                     {0.7L, 8.602252116471695076665e+32L},
+                     {0.8L, 1.062267665789207148447e+33L},
+                     {0.9L, 1.663092515864551144796e+33L}};
+  Table[{50, 50}] = {{0.0L, 2.725392139750729502972e+78L},
+                     {0.1L, 2.119868219036661812206e+78L},
+                     {0.2L, 9.822223793168123441126e+77L},
+                     {0.3L, 2.579074022415003126022e+77L},
+                     {0.4L, 3.486672453344334274446e+76L},
+                     {0.5L, 2.050976025703723878091e+75L},
+                     {0.6L, 3.889809643178187639203e+73L},
+                     {0.7L, 1.332549655956656855416e+71L},
+                     {0.8L, 2.202883440310146416308e+67L},
+                     {0.9L, 2.536699497443158134445e+60L}};
+  Table[{100, 10}] = {{0.0L, -8.249459718167038282000e+18L},
+                      {0.1L, 6.862785522503444437500e+18L},
+                      {0.2L, -2.346881035809095210750e+18L},
+                      {0.3L, -4.908935838805249369500e+18L},
+                      {0.4L, 8.237083461842661207000e+18L},
+                      {0.5L, 4.238712302196339605250e+18L},
+                      {0.6L, -9.593611165912185080625e+17L},
+                      {0.7L, 4.435504848749747270500e+18L},
+                      {0.8L, 1.577867354567400466500e+18L},
+                      {0.9L, 2.564339595766638741406e+17L}};
+  Table[{100, 20}] = {{0.0L, 7.772743983615954055061e+38L},
+                      {0.1L, -7.050326645170251933206e+38L},
+                      {0.2L, 4.396625130744397503342e+38L},
+                      {0.3L, 1.147449689190144464586e+38L},
+                      {0.4L, -7.612058604384389981878e+38L},
+                      {0.5L, 3.075049703999956739724e+38L},
+                      {0.6L, 8.268000557412077024506e+38L},
+                      {0.7L, 7.902223685310973476875e+38L},
+                      {0.8L, -8.396389511696198700493e+38L},
+                      {0.9L, 7.160434487878827530346e+37L}};
+  Table[{100, 50}] = {{0.0L, -1.171214503157837848611e+98L},
+                      {0.1L, 9.011933855018056106092e+97L},
+                      {0.2L, -2.779114958812187282420e+97L},
+                      {0.3L, -3.352060206747950450847e+97L},
+                      {0.4L, 7.143138509374017832078e+97L},
+                      {0.5L, -8.225262033972737461155e+97L},
+                      {0.6L, 5.936404592566909963016e+97L},
+                      {0.7L, 1.808933390346613713076e+97L},
+                      {0.8L, -1.102979797745434149687e+98L},
+                      {0.9L, 9.323127851689438561522e+96L}};
+  Table[{100, 100}] = {{0.0L, 6.666308670072953743181e+186L},
+                       {0.1L, 4.033157190805700092363e+186L},
+                       {0.2L, 8.658587914752715633363e+185L},
+                       {0.3L, 5.969734752682171415182e+184L},
+                       {0.4L, 1.091062733045894109826e+183L},
+                       {0.5L, 3.775274968288970779608e+180L},
+                       {0.6L, 1.357951059028910874245e+177L},
+                       {0.7L, 1.593654685059493405648e+172L},
+                       {0.8L, 4.355223604158650443018e+164L},
+                       {0.9L, 5.775179225575939991155e+150L}};
+  compareWithTable(Table, RelTolerance, AbsTolerance);
+}
+
+template <class T>
+void compareWithBoost(T RelTolerance, T AbsTolerance) {
+  // Generated with gcc 8.3.1 and boost 1.66.0
+  std::map<std::pair<unsigned, unsigned>, std::map<long double, long double> >
+      Table;
+  Table[{10, 10}] = {{0.0L, 6.547290750000000000000e+08L},
+                     {0.1L, 6.226408357052308425191e+08L},
+                     {0.2L, 5.338482120799027199391e+08L},
+                     {0.3L, 4.085719891315887824458e+08L},
+                     {0.4L, 2.738155182015052799543e+08L},
+                     {0.5L, 1.553702785400390625000e+08L},
+                     {0.6L, 7.030099912163327999588e+07L},
+                     {0.7L, 2.258980634388728249723e+07L},
+                     {0.8L, 3.958896348126719997708e+06L},
+                     {0.9L, 1.621174007878425002218e+05L}};
+  Table[{20, 10}] = {{0.0L, -1.612052956674316406250e+12L},
+                     {0.1L, 3.560417563907174737155e+11L},
+                     {0.2L, 1.466609267659881798267e+12L},
+                     {0.3L, -1.092420454297714175999e+12L},
+                     {0.4L, -8.999734873105529626012e+11L},
+                     {0.5L, 1.744502295946207410216e+12L},
+                     {0.6L, -4.244719151950563156009e+11L},
+                     {0.7L, -1.548364524949579275131e+12L},
+                     {0.8L, 2.392757933281052242041e+12L},
+                     {0.9L, 9.900174766949915348291e+11L}};
+  Table[{20, 20}] = {{0.0L, 3.198309867728777708175e+23L},
+                     {0.1L, 2.892494114697687008379e+23L},
+                     {0.2L, 2.126340780079752327168e+23L},
+                     {0.3L, 1.245473413229778327224e+23L},
+                     {0.4L, 5.593883258401251196928e+22L},
+                     {0.5L, 1.801080697817960690381e+22L},
+                     {0.6L, 3.687400224900788130304e+21L},
+                     {0.7L, 3.807345588062060328960e+20L},
+                     {0.8L, 1.169352761683327099800e+19L},
+                     {0.9L, 1.960904971202388849219e+16L}};
+  Table[{50, 10}] = {{0.0L, 1.145223824924634340332e+16L},
+                     {0.1L, 2.794368162749986070557e+15L},
+                     {0.2L, -9.926439446673560716797e+15L},
+                     {0.3L, -9.451384032851574119141e+15L},
+                     {0.4L, 8.892017018669203523560e+14L},
+                     {0.5L, 9.180035388465756391602e+15L},
+                     {0.6L, 1.217669075554222970703e+16L},
+                     {0.7L, 1.215353501150698341113e+16L},
+                     {0.8L, 9.323113753066150515747e+14L},
+                     {0.9L, -8.994661710093201604004e+15L}};
+  Table[{50, 20}] = {{0.0L, -7.873593569733218616572e+32L},
+                     {0.1L, 5.415135829089829046900e+31L},
+                     {0.2L, 7.930926605643441882380e+32L},
+                     {0.3L, -4.005406723624505746709e+31L},
+                     {0.4L, -8.232432527977407807587e+32L},
+                     {0.5L, -4.016986081427444519738e+32L},
+                     {0.6L, 4.086012875680830746546e+32L},
+                     {0.7L, 8.602252116471695071036e+32L},
+                     {0.8L, 1.062267665789207148447e+33L},
+                     {0.9L, 1.663092515864551146625e+33L}};
+  Table[{50, 50}] = {{0.0L, 2.725392139750729502972e+78L},
+                     {0.1L, 2.119868219036661818834e+78L},
+                     {0.2L, 9.822223793168123432087e+77L},
+                     {0.3L, 2.579074022415003118238e+77L},
+                     {0.4L, 3.486672453344334267855e+76L},
+                     {0.5L, 2.050976025703723882014e+75L},
+                     {0.6L, 3.889809643178187647171e+73L},
+                     {0.7L, 1.332549655956656853859e+71L},
+                     {0.8L, 2.202883440310146415724e+67L},
+                     {0.9L, 2.536699497443158141239e+60L}};
+  Table[{100, 10}] = {{0.0L, -8.249459718167038282000e+18L},
+                      {0.1L, 6.862785522503444443500e+18L},
+                      {0.2L, -2.346881035809095210750e+18L},
+                      {0.3L, -4.908935838805249365500e+18L},
+                      {0.4L, 8.237083461842661209000e+18L},
+                      {0.5L, 4.238712302196339609250e+18L},
+                      {0.6L, -9.593611165912185088125e+17L},
+                      {0.7L, 4.435504848749747265500e+18L},
+                      {0.8L, 1.577867354567400459375e+18L},
+                      {0.9L, 2.564339595766638627500e+17L}};
+  Table[{100, 20}] = {{0.0L, 7.772743983615954055061e+38L},
+                      {0.1L, -7.050326645170251938371e+38L},
+                      {0.2L, 4.396625130744397503342e+38L},
+                      {0.3L, 1.147449689190144463479e+38L},
+                      {0.4L, -7.612058604384389980402e+38L},
+                      {0.5L, 3.075049703999956748025e+38L},
+                      {0.6L, 8.268000557412077031146e+38L},
+                      {0.7L, 7.902223685310973468759e+38L},
+                      {0.8L, -8.396389511696198700493e+38L},
+                      {0.9L, 7.160434487878827394302e+37L}};
+  Table[{100, 50}] = {{0.0L, -1.171214503157837848611e+98L},
+                      {0.1L, 9.011933855018056135735e+97L},
+                      {0.2L, -2.779114958812187282420e+97L},
+                      {0.3L, -3.352060206747950445289e+97L},
+                      {0.4L, 7.143138509374017815033e+97L},
+                      {0.5L, -8.225262033972737473753e+97L},
+                      {0.6L, 5.936404592566909988953e+97L},
+                      {0.7L, 1.808933390346613708259e+97L},
+                      {0.8L, -1.102979797745434147612e+98L},
+                      {0.9L, 9.323127851689438586533e+96L}};
+  Table[{100, 100}] = {{0.0L, 6.666308670072953744596e+186L},
+                       {0.1L, 4.033157190805700115242e+186L},
+                       {0.2L, 8.658587914752715622749e+185L},
+                       {0.3L, 5.969734752682171379064e+184L},
+                       {0.4L, 1.091062733045894106141e+183L},
+                       {0.5L, 3.775274968288970794454e+180L},
+                       {0.6L, 1.357951059028910880506e+177L},
+                       {0.7L, 1.593654685059493402464e+172L},
+                       {0.8L, 4.355223604158650442269e+164L},
+                       {0.9L, 5.775179225575940023096e+150L}};
+  compareWithTable(Table, RelTolerance, AbsTolerance);
+}
+
+template <class T>
+void testReferenceValues(T RelTolerance, T AbsTolerance) {
+  compareWithMSVC<T>(RelTolerance, AbsTolerance);
+  compareWithGCC<T>(RelTolerance, AbsTolerance);
+  compareWithBoost<T>(RelTolerance, AbsTolerance);
+}
+
+int main(int, char**) {
+  compareWithAnalyticExpressions<float>(1e-6f, 1e-6f);
+  compareWithAnalyticExpressions<double>(1e-12, 1e-12);
+
+  testReferenceValues<float>(1e-6f, 0.0f);
+  testReferenceValues<double>(1e-12, 0.0);
+}

libcxx/test/std/numerics/legendre/assoc_legendre_exceptions.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
+//
+//===----------------------------------------------------------------------===//
+
+// UNSUPPORTED: c++98, c++03, c++11, c++14, libcpp-no-exceptions
+
+// <cmath>
+
+// floating_point assoc_legendre (unsigned l, unsigned m, arithmetic x);
+// float          assoc_legendref(unsigned l, unsigned m, float x);
+// long double    assoc_legendrel(unsigned l, unsigned m, long double x);
+
+#include <cmath>
+#include <stdexcept>
+#include <cassert>
+
+/// \brief std::legendre(n,x) throw if x is not in the range [-1,1] and not NaN.
+///
+/// \tparam T one of float, double, or long double
+template <class T>
+void testInvalidArguments() {
+  const T InvalidArguments[] = {T(-5), T(2)};
+  for (unsigned l = 0; l < 128; ++l) {
+    for (unsigned m = 0; m < 128; ++m) {
+      for (T x : InvalidArguments) {
+        bool ThrowsDomainError = false;
+        try {
+          const T Value = std::assoc_legendre(l, m, x);
+          (void)Value;
+        } catch (const std::domain_error&) {
+          ThrowsDomainError = true;
+        } catch (...) {
+        }
+        assert(ThrowsDomainError);
+      }
+    }
+  }
+}
+
+int main(int, char**) {
+  testInvalidArguments<float>();
+  testInvalidArguments<double>();
+  testInvalidArguments<long double>();
+}

libcxx/test/std/numerics/legendre/legendre_basic.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
+//
+//===----------------------------------------------------------------------===//
+
+// UNSUPPORTED: c++98, c++03, c++11, c++14
+
+// <cmath>
+
+// floating_point legendre (unsigned n, arithmetic x);
+// float          legendref(unsigned n, float x);
+// long double    legendrel(unsigned n, long double x);
+
+#include <cmath>
+#include <cassert>
+#include <limits>
+
+/// \brief Tests if the Legendre polynomial behaves correctly
+/// under the transformation x -> -x
+///
+/// The expected property is
+/// \f$ P_{n}(-x) = (-1)^{n} P_{n}(x) \f$.
+/// For odd \f$ n \f$, this test checks for a root at \f$ x = 0\f$.
+/// \tparam T one of float, double, or long double
+/// \param [in] n degree of the Legendre polynomial
+/// \param [in] x argument of the Legendre polynomial
+/// \param [in] RelTolerance relative tolerance for the value of legendre
+/// \param [in] AbsTolerance absolute tolerance for the value of legendre
+template <class T>
+void testParity(unsigned n, T x, T RelTolerance, T AbsTolerance) {
+  const T Parity = n % 2 == 0 ? T(1) : T(-1);
+  const T ExpectedValue = Parity * std::legendre(n, x);
+  const T Tolerance = RelTolerance * std::abs(ExpectedValue) + AbsTolerance;
+  const T Error = std::abs(std::legendre(n, -x) - ExpectedValue);
+  assert(Error <= Tolerance);
+}
+
+template <class T>
+void testParity(const T RelTolerance, const T AbsTolerance) {
+  const T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)};
+  for (unsigned n = 0; n < 128; ++n) {
+    for (T x : Samples)
+      testParity<T>(n, x, RelTolerance, AbsTolerance);
+  }
+}
+
+/// \brief Tests if the values of the Legendre polynomial at x = -1
+/// and x = -1 are correct.
+///
+/// The expected values are:
+/// \f$ P_{n}(1) = 1 \f$ and \f$ P_{n}(1) = (-1)^{n} \f$.
+/// These properties can be easily proven by induction.
+/// \tparam T one of float, double, or long double
+/// \param [in] n degree of the legendre polynomial.
+/// \param [in] RelTolerance relative tolerance for the value of legendre(x)
+/// \param [in] AbsTolerance absolute tolerance for the value of legendre(x)
+template <class T>
+void testBoundaryValues(unsigned n, T RelTolerance, T AbsTolerance) {
+  const T Parity = n % 2 == 0 ? T(1) : T(-1);
+  const T Tolerance = RelTolerance + AbsTolerance;
+  const T ErrorRight = std::abs(std::legendre(n, T(1)) - T(1));
+  const T ErrorLeft = std::abs(std::legendre(n, T(-1)) - Parity);
+  assert(ErrorRight <= Tolerance);
+  assert(ErrorLeft <= Tolerance);
+}
+
+template <class T>
+void testBoundaryValues(const T RelTolerance, const T AbsTolerance) {
+  for (unsigned n = 0; n < 128; ++n)
+    testBoundaryValues<T>(n, RelTolerance, AbsTolerance);
+}
+
+/// \brief Tests if an analytic upper bound holds
+/// for the implementation of legendre.
+///
+/// The upper bound is given by
+/// \f[
+/// \left| P_n (x) \right| \leq
+/// \left( 1 + n(n+1) (1-x^2)\right)^{-\frac{1}{4}}.
+/// \f]
+/// In particular, we check for \f$ \left| P_n (x) \right| \leq 1 \f$.
+/// See Eq. A1, p. 1435, in Martin, A. (1963).
+/// Unitarity and high-energy behavior of scattering amplitudes.
+/// Physical Review, 129(3), 1432.
+/// \tparam T one of float, double, or long double
+/// \param [in] n degree of the Legendre polynomial
+/// \param [in] x argument of the Legendre polynomial
+/// \param [in] RelTolerance relative tolerance for the value of legendre
+/// \param [in] AbsTolerance absolute tolerance for the value of legendre
+template <class T>
+void testUpperBounds(unsigned n, T x, T RelTolerance, T AbsTolerance) {
+  const T UpperBound =
+      T(1) /
+      std::sqrt(std::sqrt(T(1) + T(n) * T(n + 1u) * (T(1) - x) * (T(1) + x)));
+  const T LegendreAbs = std::abs(std::legendre(n, x));
+  const T UpperBoundWithTolerance =
+      UpperBound + UpperBound * RelTolerance + AbsTolerance;
+  assert(LegendreAbs <= UpperBoundWithTolerance);
+}
+
+template <class T>
+void testUpperBounds(T RelTolerance, T AbsTolerance) {
+  const T Samples[] = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)};
+  for (unsigned n = 0; n < 128; ++n) {
+    for (T x : Samples)
+      testUpperBounds<T>(n, x, RelTolerance, AbsTolerance);
+  }
+}
+
+/// \brief std::legendre(n,x) must silently return NaN if x is NaN.
+///
+/// \tparam T one of float, double, or long double
+template <class T>
+void testNaNPropagation() {
+  for (unsigned n = 0; n < 128; ++n) {
+    const T ShouldBeNaN = std::legendre(n, std::numeric_limits<T>::quiet_NaN());
+    assert(std::isnan(ShouldBeNaN));
+  }
+}
+
+int main(int, char**) {
+  testParity<float>(1e-6f, 0.0f);
+  testParity<double>(1e-12, 0.0);
+
+  testBoundaryValues<float>(1e-6f, 0.0f);
+  testBoundaryValues<double>(1e-12, 0.0);
+
+  testUpperBounds<float>(1e-6f, 0.0f);
+  testUpperBounds<double>(1e-12, 0.0);
+
+  testNaNPropagation<float>();
+  testNaNPropagation<double>();
+}

libcxx/test/std/numerics/legendre/legendre_comparisons.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
+//
+//===----------------------------------------------------------------------===//
+
+// UNSUPPORTED: c++98, c++03, c++11, c++14
+
+// <cmath>
+
+// floating_point legendre (unsigned n, arithmetic x);
+// float          legendref(unsigned n, float x);
+// long double    legendrel(unsigned n, long double x);
+
+#include <cmath>
+#include <map>
+#include <vector>
+#include <cassert>
+
+template <class Real, class Iterator>
+Real evaluatePolynomial(Iterator CoefficientsBegin, Iterator CoefficientsEnd,
+                        Real x) {
+  Real Result = Real(0);
+  for (Iterator Iter = CoefficientsBegin; Iter != CoefficientsEnd; ++Iter) {
+    Result = Result * x + *Iter;
+  }
+  return Result;
+}
+
+template <class Real, class Container>
+Real evaluatePolynomial(const Container& Coefficients, Real x) {
+  return evaluatePolynomial(std::begin(Coefficients), std::end(Coefficients),
+                            x);
+}
+
+/// \brief Checks if the Legendre polynomials up to n=16
+/// have roots at values from the literature.
+///
+/// The test checks for a change of sign in the interval that is
+/// centered around the expected root and has the width of the
+/// provided tolerance.
+/// The expected values were taken from:
+/// Lowan, A. N., Davids, N., & Levenson, A. (1942).
+/// Table of the zeros of the Legendre polynomials of order 1-16 and
+/// the weight coefficients for Gauss� mechanical quadrature formula.
+/// Bulletin of the American Mathematical Society, 48(10), 739-743.
+/// For n=11 and n=12, I had to correct some misprints (in accordance
+/// with wolframalpha).
+/// \tparam T one of float, double, or long double
+/// \param [in] RelTolerance Tolerance for x; contributes to the
+/// width of the interval where legendre changes its sign.
+/// \param [in] AbsTolerance Tolerance for x; contributes
+/// to the
+/// width of the interval where legendre changes its sign.
+template <class T>
+void compareZerosWithLowan1942(T RelTolerance, T AbsTolerance) {
+  std::vector<std::vector<T> > LegendreZeros = {
+      {},
+      {T(0)},
+      {T(0.577350269189626L)},
+      {T(0), T(0.774596669241483)},
+      {T(0.339981043584856), T(0.86113631159405)},
+      {T(0), T(0.538469310105683), T(0.90617984593866)},
+      {T(0.238619186083197), T(0.661209386466265), T(0.93246951420315)},
+      {T(0), T(0.405845151377397), T(0.741531185599394), T(0.949107912342759)},
+      {T(0.183434642495650), T(0.525532409916329), T(0.796666477413627),
+       T(0.960289856497536)},
+      {T(0), T(0.324253423403809), T(0.613371432700590), T(0.836031107326636),
+       T(0.968160239507626)},
+      {T(0.148874338981631), T(0.433395394129247), T(0.679409568299024),
+       T(0.865063366688985), T(0.973906528517172)},
+      {T(0), T(0.269543155952345), T(0.51909612920681), T(0.730152005574049),
+       T(0.887062599768095), T(0.978228658146057)},
+      {T(0.125233408511469), T(0.367831498998180), T(0.587317954286617),
+       T(0.769902674194305), T(0.904117256370475), T(0.981560634246719)},
+      {T(0), T(0.230458315955135), T(0.448492751036447), T(0.642349339440340),
+       T(0.801578090733310), T(0.917598399222978), T(0.984183054718588)},
+      {T(0.108054948707344), T(0.319112368927890), T(0.515248636358154),
+       T(0.687292904811685), T(0.827201315069765), T(0.928434883663574),
+       T(0.986283808696812)},
+      {T(0), T(0.201194093997435), T(0.394151347077563), T(0.570972172608539),
+       T(0.724417731360170), T(0.848206583410427), T(0.937273392400706),
+       T(0.987992518020485)},
+      {T(0.095012509837637), T(0.281603550779259), T(0.458016777657227),
+       T(0.617876244402644), T(0.755404408355003), T(0.865631202387832),
+       T(0.944575023073233), T(0.989400934991650)}};
+  for (size_t n = 0; n < LegendreZeros.size(); ++n) {
+    for (T Root : LegendreZeros[n]) {
+      const T HalfWidthInterval = T(0.5) * (RelTolerance * Root + AbsTolerance);
+      const T LeftOfRoot = Root - HalfWidthInterval;
+      const T RightOfRoot = Root + HalfWidthInterval;
+      assert(LeftOfRoot != Root);
+      assert(RightOfRoot != Root);
+      const T LegendreAtLeft =
+          std::legendre(static_cast<unsigned int>(n), LeftOfRoot);
+      const T LegendreAtRight =
+          std::legendre(static_cast<unsigned int>(n), RightOfRoot);
+      assert(std::signbit(LegendreAtLeft) ^ std::signbit(LegendreAtRight));
+    }
+  }
+}
+
+/// \brief Compares std::legendre with analytic expressions
+/// for n = 0..6.
+///
+/// \tparam T one of float, double, or long double
+/// \tparam TestFunction a callable with arguments (unsigned,T,T)
+/// \param [in] x argument of legendre
+/// \param [in] ValueSink will be invoked with the arguments of
+/// legendre and its expected outcome, i.e., (n,x,expected).
+template <class T, class TestFunction>
+void compareWithAnalyticExpressions(T x, TestFunction&& ValueSink) {
+  std::vector<std::vector<T> > VectorOfCoeffients;
+  VectorOfCoeffients.push_back({T(1)});
+  VectorOfCoeffients.push_back({T(1), T(0)});
+  VectorOfCoeffients.push_back({T(3) / T(2), T(0), -T(1) / T(2)});
+  VectorOfCoeffients.push_back({T(5) / T(2), T(0), -T(3) / T(2), T(0)});
+  VectorOfCoeffients.push_back(
+      {T(35) / T(8), T(0), -T(30) / T(8), T(0), T(3) / T(8)});
+  VectorOfCoeffients.push_back(
+      {T(63) / T(8), T(0), -T(70) / T(8), T(0), T(15) / T(8), T(0)});
+  VectorOfCoeffients.push_back({T(231) / T(16), T(0), -T(315) / T(16), T(0),
+                                T(105) / T(16), T(0), -T(5) / T(16)});
+
+  for (const auto& Coefficients : VectorOfCoeffients) {
+    const unsigned int n = static_cast<unsigned int>(Coefficients.size()) - 1;
+    const T ExpectedValue = ::evaluatePolynomial(Coefficients, x);
+    ValueSink(n, x, ExpectedValue);
+  }
+}
+
+/// \brief Compares std::legendre with analytic expressions
+/// for n = 0..6.
+///
+/// \tparam T one of float, double, or long double
+/// \param [in] RelTolerance relative tolerance for the value of legendre
+/// \param [in] AbsTolerance absolute tolerance for the value of legendre
+template <class T>
+void compareWithAnalyticExpressions(T RelTolerance, T AbsTolerance) {
+  auto testFunction = [RelTolerance, AbsTolerance](unsigned n, T x,
+                                                   T ExpectedValue) {
+    const T Value = std::legendre(n, x);
+    const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance;
+    const T Error = std::abs(Value - ExpectedValue);
+    assert(Error <= Tolerance);
+  };
+
+  const auto Samples = {T(0), T(0.1), T(0.2), T(0.5), T(0.9), T(1.0)};
+  for (T x : Samples)
+    compareWithAnalyticExpressions(x, testFunction);
+}
+
+template <class T>
+void compareWithTable(
+    const std::map<unsigned, std::map<long double, long double> >& Table,
+    T RelTolerance, T AbsTolerance) {
+  for (const auto& Entry : Table) {
+    const unsigned n = Entry.first;
+    for (const auto& xy : Entry.second) {
+      const T Value = std::legendre(n, xy.first);
+      const T ExpectedValue = xy.second;
+      const T Tolerance = std::abs(ExpectedValue) * RelTolerance + AbsTolerance;
+      const T Error = std::abs(Value - ExpectedValue);
+      assert(Error <= Tolerance);
+    }
+  }
+}
+
+template <class T>
+void compareWithMSVC(T RelTolerance, T AbsTolerance) {
+  // Generated with msvc 191627027
+  std::map<unsigned, std::map<long double, long double> > Table;
+  Table[10u] = {
+      {0.0L, -2.46093750000000000e-01L}, {0.1L, -1.22124997387109330e-01L},
+      {0.2L, 1.29072025600000029e-01L},  {0.3L, 2.51476349516015663e-01L},
+      {0.4L, 9.68390643999999245e-02L},  {0.5L, -1.88228607177734375e-01L},
+      {0.6L, -2.43662745600000008e-01L}, {0.7L, 8.58057955316403331e-02L},
+      {0.8L, 3.00529795600000038e-01L},  {0.9L, -2.63145617855859604e-01L}};
+  Table[20u] = {
+      {0.0L, 1.76197052001953125e-01L},  {0.1L, -8.20771309445277736e-02L},
+      {0.2L, -9.80421943445946159e-02L}, {0.3L, 1.80287159479980419e-01L},
+      {0.4L, -1.01592615586281562e-01L}, {0.5L, -4.83583810673735570e-02L},
+      {0.6L, 1.59167529100981087e-01L},  {0.7L, -2.04573944638341720e-01L},
+      {0.8L, 2.24204605417413466e-01L},  {0.9L, -1.49308235309848353e-01L}};
+  Table[50u] = {
+      {0.0L, -1.12275172659217048e-01L}, {0.1L, -3.82058126613136831e-02L},
+      {0.2L, 8.34322722041973136e-02L},  {0.3L, 1.09110515747147960e-01L},
+      {0.4L, 4.15690333818253752e-02L},  {0.5L, -3.10590992396098109e-02L},
+      {0.6L, -5.88607988440020963e-02L}, {0.7L, -1.45727316458928516e-02L},
+      {0.8L, 1.38797373450931127e-01L},  {0.9L, -1.70037659943836711e-01L}};
+  Table[100u] = {
+      {0.0L, 7.95892373871787268e-02L},  {0.1L, -6.38950984347501499e-02L},
+      {0.2L, 1.46818353556592614e-02L},  {0.3L, 5.71273922028013650e-02L},
+      {0.4L, -7.22582021256844703e-02L}, {0.5L, -6.05180259618611979e-02L},
+      {0.6L, -2.37470239051331411e-02L}, {0.7L, -7.71325071997786410e-02L},
+      {0.8L, 5.08611679135842279e-02L},  {0.9L, 1.02265820558718926e-01L}};
+
+  compareWithTable<T>(Table, RelTolerance, AbsTolerance);
+}
+
+template <class T>
+void compareWithGCC(T RelTolerance, T AbsTolerance) {
+  // Generated with gcc 8.3.1
+  std::map<unsigned, std::map<long double, long double> > Table;
+  Table[10u] = {{0.0L, -2.460937500000000000000e-01L},
+                {0.1L, -1.221249973871093750063e-01L},
+                {0.2L, 1.290720256000000000328e-01L},
+                {0.3L, 2.514763495160156249876e-01L},
+                {0.4L, 9.683906439999999996141e-02L},
+                {0.5L, -1.882286071777343750000e-01L},
+                {0.6L, -2.436627456000000000027e-01L},
+                {0.7L, 8.580579553164062492864e-02L},
+                {0.8L, 3.005297956000000000312e-01L},
+                {0.9L, -2.631456178558593748918e-01L}};
+  Table[20u] = {{0.0L, 1.761970520019531250000e-01L},
+                {0.1L, -8.207713094452772853985e-02L},
+                {0.2L, -9.804219434459463679125e-02L},
+                {0.3L, 1.802871594799804655558e-01L},
+                {0.4L, -1.015926155862814718285e-01L},
+                {0.5L, -4.835838106737355701625e-02L},
+                {0.6L, 1.591675291009810431662e-01L},
+                {0.7L, -2.045739446383416052931e-01L},
+                {0.8L, 2.242046054174133682042e-01L},
+                {0.9L, -1.493082353098486828345e-01L}};
+  Table[50u] = {{0.0L, -1.122751726592170484764e-01L},
+                {0.1L, -3.820581266131368189813e-02L},
+                {0.2L, 8.343227220419731034374e-02L},
+                {0.3L, 1.091105157471479767954e-01L},
+                {0.4L, 4.156903338182522686674e-02L},
+                {0.5L, -3.105909923960982286266e-02L},
+                {0.6L, -5.886079884400199673851e-02L},
+                {0.7L, -1.457273164589237150986e-02L},
+                {0.8L, 1.387973734509308074478e-01L},
+                {0.9L, -1.700376599438367846763e-01L}};
+  Table[100u] = {{0.0L, 7.958923738717876149887e-02L},
+                 {0.1L, -6.389509843475016413679e-02L},
+                 {0.2L, 1.468183535565932807048e-02L},
+                 {0.3L, 5.712739220280141857569e-02L},
+                 {0.4L, -7.225820212568462422759e-02L},
+                 {0.5L, -6.051802596186118686668e-02L},
+                 {0.6L, -2.374702390513330860656e-02L},
+                 {0.7L, -7.713250719977911342573e-02L},
+                 {0.8L, 5.086116791358337655272e-02L},
+                 {0.9L, 1.022658205587185505935e-01L}};
+
+  compareWithTable<T>(Table, RelTolerance, AbsTolerance);
+}
+
+template <class T>
+void compareWithBoost(T RelTolerance, T AbsTolerance) {
+  // Generated with gcc 8.3.1 and boost 1.66.0
+  std::map<unsigned, std::map<long double, long double> > Table;
+  Table[10u] = {{0.0L, -2.460937500000000000000e-01L},
+                {0.1L, -1.221249973871093750063e-01L},
+                {0.2L, 1.290720256000000000057e-01L},
+                {0.3L, 2.514763495160156249876e-01L},
+                {0.4L, 9.683906439999999996818e-02L},
+                {0.5L, -1.882286071777343750000e-01L},
+                {0.6L, -2.436627455999999999756e-01L},
+                {0.7L, 8.580579553164062490831e-02L},
+                {0.8L, 3.005297955999999999499e-01L},
+                {0.9L, -2.631456178558593748647e-01L}};
+  Table[20u] = {{0.0L, 1.761970520019531250000e-01L},
+                {0.1L, -8.207713094452772851275e-02L},
+                {0.2L, -9.804219434459463677092e-02L},
+                {0.3L, 1.802871594799804655558e-01L},
+                {0.4L, -1.015926155862814718082e-01L},
+                {0.5L, -4.835838106737355701625e-02L},
+                {0.6L, 1.591675291009810431391e-01L},
+                {0.7L, -2.045739446383416053744e-01L},
+                {0.8L, 2.242046054174133682177e-01L},
+                {0.9L, -1.493082353098486828480e-01L}};
+  Table[50u] = {{0.0L, -1.122751726592170484764e-01L},
+                {0.1L, -3.820581266131368188457e-02L},
+                {0.2L, 8.343227220419731035052e-02L},
+                {0.3L, 1.091105157471479767683e-01L},
+                {0.4L, 4.156903338182522687013e-02L},
+                {0.5L, -3.105909923960982287282e-02L},
+                {0.6L, -5.886079884400199674528e-02L},
+                {0.7L, -1.457273164589237145396e-02L},
+                {0.8L, 1.387973734509308073665e-01L},
+                {0.9L, -1.700376599438367847305e-01L}};
+  Table[100u] = {{0.0L, 7.958923738717876149209e-02L},
+                 {0.1L, -6.389509843475016413679e-02L},
+                 {0.2L, 1.468183535565932810267e-02L},
+                 {0.3L, 5.712739220280141854181e-02L},
+                 {0.4L, -7.225820212568462420048e-02L},
+                 {0.5L, -6.051802596186118688023e-02L},
+                 {0.6L, -2.374702390513330857607e-02L},
+                 {0.7L, -7.713250719977911350027e-02L},
+                 {0.8L, 5.086116791358337655272e-02L},
+                 {0.9L, 1.022658205587185507019e-01L}};
+
+  compareWithTable<T>(Table, RelTolerance, AbsTolerance);
+}
+
+template <class T>
+void testReferenceValues(T RelTolerance, T AbsTolerance) {
+  compareWithMSVC<T>(RelTolerance, AbsTolerance);
+  compareWithGCC<T>(RelTolerance, AbsTolerance);
+  compareWithBoost<T>(RelTolerance, AbsTolerance);
+}
+
+int main(int, char**) {
+  compareZerosWithLowan1942<float>(0.0f, 1e-6f);
+  compareZerosWithLowan1942<double>(0.0f, 1e-12);
+
+  compareWithAnalyticExpressions<float>(1e-6f, 1e-6f);
+  compareWithAnalyticExpressions<double>(1e-12, 1e-12);
+
+  testReferenceValues<float>(1e-6f, 0.0f);
+  testReferenceValues<double>(1e-12, 0.0);
+}

libcxx/test/std/numerics/legendre/legendre_exceptions.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
+//
+//===----------------------------------------------------------------------===//
+
+// UNSUPPORTED: c++98, c++03, c++11, c++14, libcpp-no-exceptions
+
+// <cmath>
+
+// floating_point legendre (unsigned n, arithmetic x);
+// float          legendref(unsigned n, float x);
+// long double    legendrel(unsigned n, long double x);
+
+#include <cmath>
+#include <stdexcept>
+#include <cassert>
+
+/// \brief std::legendre(n,x) must throw if x is not in the range [-1,1] and not NaN.
+///
+/// \tparam T one of float, double, or long double
+template <class T>
+void testInvalidArguments() {
+  const T InvalidArguments[] = {T(-5), T(2)};
+  for (unsigned n = 0; n < 128; ++n) {
+    for (T x : InvalidArguments) {
+      bool ThrowsDomainError = false;
+      try {
+        const T Value = std::legendre(n, x);
+        (void)Value;
+      } catch (const std::domain_error&) {
+        ThrowsDomainError = true;
+      } catch (...) {
+      }
+      assert(ThrowsDomainError);
+    }
+  }
+}
+
+int main(int, char**) {
+  testInvalidArguments<float>();
+  testInvalidArguments<double>();
+  testInvalidArguments<long double>();
+}