[libc++][math] Mathematical special functions: Implementing std::legendre

libcxx/include/__math/special_functions.h

@@ -11,11 +11,13 @@
 #define _LIBCPP___MATH_SPECIAL_FUNCTIONS_H
 
 #include <__config>
+#include <__math/abs.h>
 #include <__math/copysign.h>
 #include <__math/traits.h>
 #include <__type_traits/enable_if.h>
 #include <__type_traits/is_integral.h>
 #include <limits>
+#include <stdexcept>
 
 #if !defined(_LIBCPP_HAS_NO_PRAGMA_SYSTEM_HEADER)
 #  pragma GCC system_header
@@ -49,7 +51,7 @@ _LIBCPP_HIDE_FROM_ABI _Real __hermite(unsigned __n, _Real __x) {
   }
 
   if (!__math::isfinite(__H_n)) {
-    // Overflow occured. Two possible cases:
+    // Overflow occurred. Two possible cases:
     //    n is odd:  return infinity of the same sign as x.
     //    n is even: return +Inf
     _Real __inf = std::numeric_limits<_Real>::infinity();
@@ -77,6 +79,51 @@ _LIBCPP_HIDE_FROM_ABI double hermite(unsigned __n, _Integer __x) {
   return std::hermite(__n, static_cast<double>(__x));
 }
 
+template <class _Real>
+_LIBCPP_HIDE_FROM_ABI _Real __legendre(unsigned __n, _Real __x) {
+  // The Legendre polynomial P_n(x).
+  // The implementation is based on the recurrence formula: (n+1) P_{n+1}(x) = (2n+1) x P_n(x) - n P_{n-1}(x).
+  // Press, William H., et al. Numerical recipes 3rd edition: The art of scientific computing.
+  // Cambridge university press, 2007, p. 183.
+
+  // NOLINTBEGIN(readability-identifier-naming)
+  if (__math::isnan(__x))
+    return __x;
+
+  if (__math::fabs(__x) > 1)
+    std::__throw_domain_error("Argument `x` of Legendre function is out of range: `|x| <= 1`.");
+
+  _Real __P_0{1};
+  if (__n == 0)
+    return __P_0;
+
+  _Real __P_n_prev = __P_0;
+  _Real __P_n      = __x;
+  for (unsigned __i = 1; __i < __n; ++__i) {
+    _Real __P_n_next = ((2 * __i + 1) * __x * __P_n - __i * __P_n_prev) / static_cast<_Real>(__i + 1);
+    __P_n_prev       = __P_n;
+    __P_n            = __P_n_next;
+  }
+
+  return __P_n;
+  // NOLINTEND(readability-identifier-naming)
+}
+
+inline _LIBCPP_HIDE_FROM_ABI double legendre(unsigned __n, double __x) { return std::__legendre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI float legendre(unsigned __n, float __x) { return std::__legendre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI long double legendre(unsigned __n, long double __x) { return std::__legendre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI float legendref(unsigned __n, float __x) { return std::legendre(__n, __x); }
+
+inline _LIBCPP_HIDE_FROM_ABI long double legendrel(unsigned __n, long double __x) { return std::legendre(__n, __x); }
+
+template <class _Integer, std::enable_if_t<std::is_integral_v<_Integer>, int> = 0>
+_LIBCPP_HIDE_FROM_ABI double legendre(unsigned __n, _Integer __x) {
+  return std::legendre(__n, static_cast<double>(__x));
+}
+
 #endif // _LIBCPP_STD_VER >= 17
 
 _LIBCPP_END_NAMESPACE_STD

libcxx/include/cmath

@@ -224,6 +224,14 @@ int ilogb (arithmetic x);
 int ilogbf(float x);
 int ilogbl(long double x);
 
+double         legendre(unsigned n, double x);                    // C++17
+float          legendre(unsigned n, float x);                     // C++17
+long double    legendre(unsigned n, long double x);               // C++17
+float          legendref(unsigned n, float x);                    // C++17
+long double    legendrel(unsigned n, long double x);              // C++17
+template <class Integer>
+double         legendre(unsigned n, Integer x);                   // C++17
+
 floating_point lgamma (arithmetic x);
 float          lgammaf(float x);
 long double    lgammal(long double x);

libcxx/modules/std/cmath.inc

@@ -346,12 +346,14 @@ export namespace std {
   using std::laguerre;
   using std::laguerref;
   using std::laguerrel;
+#endif
 
   // [sf.cmath.legendre], Legendre polynomials
   using std::legendre;
   using std::legendref;
   using std::legendrel;
 
+#if 0
   // [sf.cmath.riemann.zeta], Riemann zeta function
   using std::riemann_zeta;
   using std::riemann_zetaf;

libcxx/test/std/numerics/c.math/legendre.pass.cpp

@@ -0,0 +1,179 @@
+//===----------------------------------------------------------------------===//
+//
+// 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++03, c++11, c++14
+
+// <cmath>
+
+// double         legendre(unsigned n, double x);
+// float          legendre(unsigned n, float x);
+// long double    legendre(unsigned n, long double x);
+// float          legendref(unsigned n, float x);
+// long double    legendrel(unsigned n, long double x);
+// template <class Integer>
+// double         legendre(unsigned n, Integer x);
+
+#include <array>
+#include <cassert>
+#include <cmath>
+#include <limits>
+
+#include "type_algorithms.h"
+
+inline constexpr unsigned g_max_n = 128;
+
+template <class T>
+std::array<T, 7> sample_points() {
+  return {-1.0, -0.5, -0.1, 0.0, 0.1, 0.5, 1.0};
+}
+
+template <class Real>
+class CompareFloatingValues {
+private:
+  Real tol;
+
+public:
+  CompareFloatingValues() {
+    if (std::is_same_v<Real, float>)
+      tol = 1e-6f;
+    else if (std::is_same_v<Real, double>)
+      tol = 1e-9;
+    else
+      tol = 1e-12l;
+  }
+
+  bool operator()(Real result, Real expected) const {
+    if (std::isinf(expected) && std::isinf(result))
+      return result == expected;
+
+    if (std::isnan(expected) || std::isnan(result))
+      return false;
+
+    return std::abs(result - expected) < (tol + std::abs(expected) * tol);
+  }
+};
+
+template <class Real>
+void test() {
+  { // checks if NaNs are reported correctly (i.e. output == input for input == NaN)
+    using nl = std::numeric_limits<Real>;
+    for (Real NaN : {nl::quiet_NaN(), nl::signaling_NaN()})
+      for (unsigned n = 0; n < g_max_n; ++n)
+        assert(std::isnan(std::legendre(n, NaN)));
+  }
+
+  { // simple sample points for n=0..127 should not produce NaNs.
+    for (Real x : sample_points<Real>())
+      for (unsigned n = 0; n < g_max_n; ++n)
+        assert(!std::isnan(std::legendre(n, x)));
+  }
+
+  { // For x with abs(x) > 1 an domain_error exception should be thrown.
+#ifndef _LIBCPP_NO_EXCEPTIONS
+    for (double absX : {2.0, 7.77, 42.42, std::numeric_limits<double>::infinity()})
+      for (Real x : {-absX, absX})
+        for (unsigned n = 0; n < g_max_n; ++n) {
+          bool throws = false;
+          try {
+            std::legendre(n, x);
+          } catch (const std::domain_error&) {
+            throws = true;
+          }
+          assert(throws);
+        }
+#endif // _LIBCPP_NO_EXCEPTIONS
+  }
+
+  { // check against analytic polynoms for order n=0..6
+    for (Real x : sample_points<Real>()) {
+      const auto l0 = [](Real) -> Real { return 1; };
+      const auto l1 = [](Real y) -> Real { return y; };
+      const auto l2 = [](Real y) -> Real { return (3 * y * y - 1) / 2; };
+      const auto l3 = [](Real y) -> Real { return (5 * y * y - 3) * y / 2; };
+      const auto l4 = [](Real y) -> Real { return (35 * std::pow(y, 4) - 30 * y * y + 3) / 8; };
+      const auto l5 = [](Real y) -> Real { return (63 * std::pow(y, 4) - 70 * y * y + 15) * y / 8; };
+      const auto l6 = [](Real y) -> Real {
+        return (231 * std::pow(y, 6) - 315 * std::pow(y, 4) + 105 * y * y - 5) / 16;
+      };
+
+      const CompareFloatingValues<Real> compare;
+      assert(compare(std::legendre(0, x), l0(x)));
+      assert(compare(std::legendre(1, x), l1(x)));
+      assert(compare(std::legendre(2, x), l2(x)));
+      assert(compare(std::legendre(3, x), l3(x)));
+      assert(compare(std::legendre(4, x), l4(x)));
+      assert(compare(std::legendre(5, x), l5(x)));
+      assert(compare(std::legendre(6, x), l6(x)));
+    }
+  }
+
+  { // checks std::legendref for bitwise equality with std::legendre(unsigned, float)
+    if constexpr (std::is_same_v<Real, float>)
+      for (unsigned n = 0; n < g_max_n; ++n)
+        for (float x : sample_points<float>())
+          assert(std::legendre(n, x) == std::legendref(n, x));
+  }
+
+  { // checks std::legendrel for bitwise equality with std::legendre(unsigned, long double)
+    if constexpr (std::is_same_v<Real, long double>)
+      for (unsigned n = 0; n < g_max_n; ++n)
+        for (long double x : sample_points<long double>()) {
+          assert(std::legendre(n, x) == std::legendrel(n, x));
+          assert(std::legendre(n, x) == std::legendrel(n, x));
+        }
+  }
+
+  { // evaluation at x=1: P_n(1) = 1
+    const CompareFloatingValues<Real> compare;
+    for (unsigned n = 0; n < g_max_n; ++n)
+      assert(compare(std::legendre(n, Real{1}), 1));
+  }
+
+  { // evaluation at x=-1: P_n(-1) = (-1)^n
+    const CompareFloatingValues<Real> compare;
+    for (unsigned n = 0; n < g_max_n; ++n)
+      assert(compare(std::legendre(n, Real{-1}), std::pow(-1, n)));
+  }
+
+  { // evaluation at x=0:
+    //    P_{2n  }(0) = (-1)^n (2n-1)!! / (2n)!!
+    //    P_{2n+1}(0) = 0
+    const CompareFloatingValues<Real> compare;
+    Real doubleFactorials{1};
+    for (unsigned n = 0; n < g_max_n; ++n) {
+      if (n & 1) // odd
+        doubleFactorials *= n;
+      else if (n != 0) // even and not zero
+        doubleFactorials /= n;
+
+      assert(compare(std::legendre(n, Real{0}), (n & 1) ? Real{0} : std::pow(-1, n / 2) * doubleFactorials));
+    }
+  }
+}
+
+struct TestFloat {
+  template <class Real>
+  void operator()() {
+    test<Real>();
+  }
+};
+
+struct TestInt {
+  template <class Integer>
+  void operator()() {
+    // checks that std::legendre(unsigned, Integer) actually wraps std::legendre(unsigned, double)
+    for (unsigned n = 0; n < g_max_n; ++n)
+      for (Integer x : {-1, 0, 1})
+        assert(std::legendre(n, x) == std::legendre(n, static_cast<double>(x)));
+  }
+};
+
+int main() {
+  types::for_each(types::floating_point_types(), TestFloat());
+  types::for_each(types::type_list<short, int, long, long long>(), TestInt());
+}