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

libcxx/include/CMakeLists.txt

@@ -516,6 +516,7 @@ set(files
   __math/remainder.h
   __math/roots.h
   __math/rounding_functions.h
+  __math/special_functions.h
   __math/traits.h
   __math/trigonometric_functions.h
   __mbstate_t.h

libcxx/include/__math/special_functions.h

@@ -0,0 +1,105 @@
+// -*- C++ -*-
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef _LIBCPP___MATH_SPECIAL_FUNCTIONS_H
+#define _LIBCPP___MATH_SPECIAL_FUNCTIONS_H
+
+#include <__config>
+#include <__math/abs.h>
+#include <__math/copysign.h>
+#include <__math/roots.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
+#endif
+
+_LIBCPP_BEGIN_NAMESPACE_STD
+
+#if _LIBCPP_STD_VER >= 17
+
+template <class _Real>
+_LIBCPP_HIDE_FROM_ABI _Real __assoc_legendre(unsigned __l, unsigned __m, _Real __x) {
+  // The associated Legendre polynomial P_l^m(x).
+  // The implementation is based on well-known recurrence formulas.
+  // Resources:
+  //    - Book: Press - Numerical Recipes - 3rd Edition - p. 294
+  //    - Excellect write-up: https://justinwillmert.com/articles/2020/calculating-legendre-polynomials/
+
+  // NOLINTBEGIN(readability-identifier-naming)
+  if (__l < __m)
+    return 0;
+  /*
+  // fixme: uncomment once `std::legendre` is integrated and git rebased.
+  if (__m == 0)
+    return std::legendre(__l, __x);
+  */
+  if (__math::isnan(__x))
+    return __x;
+  else if (__math::fabs(__x) > 1)
+    std::__throw_domain_error("Argument `x` of associated Legendre function is out of range: `|x| <= 1`.");
+
+  _Real __Pmm = 1; // init with P_0^0
+  // Compute P_m^m: Explicit loop unrolling to work around computing square root.
+  // Note: (1-x)*(1+x) is more accurate than (1-x*x)
+  for (unsigned __n = 0; __n < __m / 2; ++__n)
+    __Pmm *= (4 * __n + 1) * (4 * __n + 3) * (1 - __x) * (1 + __x);
+  // Odd m case: Cannot combine two iterations. Thus, needs to compute square root.
+  if (__m & 1)
+    __Pmm *= (2 * __m - 1) * __math::sqrt((1 - __x) * (1 + __x));
+
+  if (__l == __m)
+    return __Pmm;
+
+  // Factoring P^m_m out and multiplying it afterwards. Possible as recursion is linear.
+  _Real __Pml      = __x * (2 * __m + 1) /* * __Pmm */; // init with P^m_{m+1}
+  _Real __Pml_prev = 1 /* * __Pmm */;                   // init with P^m_m
+  for (unsigned __n = __m + 2; __n <= __l; ++__n) {
+    _Real __Pml_prev2 = __Pml_prev;                                                                       // P^m_{n-2}
+    __Pml_prev        = __Pml;                                                                            // P^m_{n-1}
+    __Pml             = (__x * (2 * __n - 1) * __Pml_prev - (__n + __m - 1) * __Pml_prev2) / (__n - __m); // P^m_n
+  }
+  return __Pml * __Pmm;
+  // NOLINTEND(readability-identifier-naming)
+}
+
+inline _LIBCPP_HIDE_FROM_ABI double assoc_legendre(unsigned __l, unsigned __m, double __x) {
+  return std::__assoc_legendre(__l, __m, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI float assoc_legendre(unsigned __l, unsigned __m, float __x) {
+  return static_cast<float>(std::__assoc_legendre(__l, __m, static_cast<double>(__x)));
+}
+
+inline _LIBCPP_HIDE_FROM_ABI long double assoc_legendre(unsigned __l, unsigned __m, long double __x) {
+  return std::__assoc_legendre(__l, __m, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI float assoc_legendref(unsigned __l, unsigned __m, float __x) {
+  return std::assoc_legendre(__l, __m, __x);
+}
+
+inline _LIBCPP_HIDE_FROM_ABI long double assoc_legendrel(unsigned __l, unsigned __m, long double __x) {
+  return std::assoc_legendre(__l, __m, __x);
+}
+
+template <class _Integer, std::enable_if_t<std::is_integral_v<_Integer>, int> = 0>
+_LIBCPP_HIDE_FROM_ABI double assoc_legendre(unsigned __l, unsigned __m, _Integer __x) {
+  return std::assoc_legendre(__l, __m, static_cast<double>(__x));
+}
+
+#endif // _LIBCPP_STD_VER >= 17
+
+_LIBCPP_END_NAMESPACE_STD
+
+#endif // _LIBCPP___MATH_SPECIAL_FUNCTIONS_H

libcxx/include/cmath

@@ -160,6 +160,14 @@ floating_point asinh (arithmetic x);
 float          asinhf(float x);
 long double    asinhl(long double x);
 
+double         assoc_legendre(unsigned l, unsigned m, double x);                    // C++17
+float          assoc_legendre(unsigned l, unsigned m, float x);                     // C++17
+long double    assoc_legendre(unsigned l, unsigned m, long double x);               // C++17
+float          assoc_legendref(unsigned l, unsigned m, float x);                    // C++17
+long double    assoc_legendrel(unsigned l, unsigned m, long double x);              // C++17
+template <class Integer>
+double         assoc_legendre(unsigned l, unsigned m, Integer x);                   // C++17
+
 floating_point atanh (arithmetic x);
 float          atanhf(float x);
 long double    atanhl(long double x);
@@ -315,6 +323,7 @@ constexpr long double lerp(long double a, long double b, long double t) noexcept
 #include <limits>
 #include <version>
 
+#include <__math/special_functions.h>
 #include <math.h>
 
 #ifndef _LIBCPP_MATH_H

libcxx/include/module.modulemap

@@ -1488,6 +1488,7 @@ module std_private_math_modulo                          [system] { header "__mat
 module std_private_math_remainder                       [system] { header "__math/remainder.h" }
 module std_private_math_roots                           [system] { header "__math/roots.h" }
 module std_private_math_rounding_functions              [system] { header "__math/rounding_functions.h" }
+module std_private_math_special_functions               [system] { header "__math/special_functions.h" }
 module std_private_math_traits                          [system] { header "__math/traits.h" }
 module std_private_math_trigonometric_functions         [system] { header "__math/trigonometric_functions.h" }
 

libcxx/modules/std/cmath.inc

@@ -268,12 +268,14 @@ export namespace std {
   using std::assoc_laguerre;
   using std::assoc_laguerref;
   using std::assoc_laguerrel;
+#endif
 
   // [sf.cmath.assoc.legendre], associated Legendre functions
   using std::assoc_legendre;
   using std::assoc_legendref;
   using std::assoc_legendrel;
 
+#if 0
   // [sf.cmath.beta], beta function
   using std::beta;
   using std::betaf;

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

@@ -0,0 +1,221 @@
+//===----------------------------------------------------------------------===//
+//
+// 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      assoc_legendre( unsigned l, unsigned m, double      x);
+// float       assoc_legendre( unsigned l, unsigned m, float       x);
+// long double assoc_legendre( unsigned l, unsigned m, long double x);
+// float       assoc_legendref(unsigned l, unsigned m, float       x);
+// long double assoc_legendrel(unsigned l, unsigned m, long double x);
+// template <class Integer>
+// double      assoc_legendre( unsigned l, unsigned m, Integer     x);
+
+#include <array>
+#include <cassert>
+#include <functional>
+#include <cmath>
+#include <limits>
+#include <iostream>
+#include <numeric>
+
+#include "type_algorithms.h"
+
+inline constexpr unsigned g_max_lm = 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 T, std::size_t N>
+std::array<T, N> range() {
+  std::array<T, N> arr;
+  std::iota(arr.begin(), arr.end(), T{0});
+  return arr;
+}
+
+template <class Real>
+void test() {
+#if 1
+  { // checks if NaNs are reported correctly (i.e. output == input for input == NaN)
+    using nl = std::numeric_limits<Real>;
+    for (unsigned l = 0; l < g_max_lm; ++l)
+      for (unsigned m = 0; m < g_max_lm; ++m)
+        for (Real NaN : {nl::quiet_NaN(), nl::signaling_NaN()})
+          assert(std::isnan(std::assoc_legendre(l, m, NaN)));
+  }
+
+  { // simple sample points for l,m=0..127 should not produce NaNs.
+    for (unsigned l = 0; l < g_max_lm; ++l)
+      for (unsigned m = 0; m < g_max_lm; ++m)
+        for (Real x : sample_points<Real>())
+          assert(!std::isnan(std::assoc_legendre(l, m, x)));
+  }
+
+  { // For x with abs(x) > 1 a domain_error exception should be thrown.
+#  ifndef _LIBCPP_NO_EXCEPTIONS
+    for (unsigned l = 0; l < g_max_lm; ++l)
+      for (unsigned m = 0; m < g_max_lm; ++m) {
+        const auto is_domain_error = [l, m](Real x) -> bool {
+          try {
+            std::assoc_legendre(l, m, x);
+          } catch (const std::domain_error&) {
+            return true;
+          }
+          return false;
+        };
+
+        Real inf = std::numeric_limits<Real>::infinity();
+        for (Real absX : {std::nextafter(Real(1), inf), Real(2.0), Real(7.77), Real(42.42), inf})
+          for (Real x : {-absX, absX})
+            assert(is_domain_error(x));
+      }
+#  endif // _LIBCPP_NO_EXCEPTIONS
+  }
+#endif
+
+#if 1
+  { // check against analytic polynoms for order n=0..6
+    using Func     = std::function<Real(Real)>;
+    const Func P00 = [](Real) { return 1; };
+
+    const Func P10 = [](Real y) -> Real { return y; };
+    const Func P11 = [](Real y) -> Real { return std::sqrt((1 - y) * (1 + y)); };
+
+    const Func P20 = [](Real y) -> Real { return (3 * y * y - 1) / 2; };
+    const Func P21 = [](Real y) -> Real { return 3 * y * std::sqrt((1 - y) * (1 + y)); };
+    const Func P22 = [](Real y) -> Real { return 3 * (1 - y) * (1 + y); };
+
+    const Func P30 = [](Real y) -> Real { return (5 * y * y - 3) * y / 2; };
+    const Func P31 = [](Real y) -> Real { return Real(1.5) * (5 * y * y - 1) * std::sqrt((1 - y) * (1 + y)); };
+    const Func P32 = [](Real y) -> Real { return 15 * y * (1 - y) * (1 + y); };
+    const Func P33 = [](Real y) -> Real { return 15 * std::pow((1 - y) * (1 + y), Real(1.5)); };
+
+    const Func P40 = [](Real y) -> Real { return (35 * std::pow(y, 4) - 30 * y * y + 3) / 8; };
+    const Func P41 = [](Real y) -> Real { return Real(2.5) * y * (7 * y * y - 3) * std::sqrt((1 - y) * (1 + y)); };
+    const Func P42 = [](Real y) -> Real { return Real(7.5) * (7 * y * y - 1) * (1 - y) * (1 + y); };
+    const Func P43 = [](Real y) -> Real { return 105 * y * std::pow((1 - y) * (1 + y), Real(1.5)); };
+    const Func P44 = [](Real y) -> Real { return 105 * std::pow((1 - y) * (1 + y), Real(2)); };
+
+    unsigned l = 0;
+    unsigned m = 0;
+    for (auto P : {P00, P10, P11, P20, P21, P22, P30, P31, P32, P33, P40, P41, P42, P43, P44}) {
+      for (Real x : sample_points<Real>()) {
+        const CompareFloatingValues<Real> compare;
+        assert(compare(std::assoc_legendre(l, m, x), P(x)));
+      }
+
+      if (l == m) {
+        ++l;
+        m = 0;
+      } else
+        ++m;
+    }
+  }
+#endif
+
+#if 1
+  { // checks std::assoc_legendref for bitwise equality with std::assoc_legendre(unsigned, float)
+    if constexpr (std::is_same_v<Real, float>)
+      for (unsigned l = 0; l < g_max_lm; ++l)
+        for (unsigned m = 0; m < g_max_lm; ++m)
+          for (float x : sample_points<float>())
+            assert(std::assoc_legendre(l, m, x) == std::assoc_legendref(l, m, x));
+  }
+
+  { // checks std::assoc_legendrel for bitwise equality with std::assoc_legendre(unsigned, long double)
+    if constexpr (std::is_same_v<Real, long double>)
+      for (unsigned l = 0; l < g_max_lm; ++l)
+        for (unsigned m = 0; m < g_max_lm; ++m)
+          for (long double x : sample_points<long double>())
+            assert(std::assoc_legendre(l, m, x) == std::assoc_legendrel(l, m, x));
+  }
+#endif
+
+#if 0
+  { // evaluation at x=1: P_n(1) = 1
+    const CompareFloatingValues<Real> compare;
+    for (unsigned n = 0; n < g_max_lm; ++n)
+      assert(compare(std::assoc_legendre(l, m, Real{1}), 1));
+  }
+
+  { // evaluation at x=-1: P_n(-1) = (-1)^n
+    const CompareFloatingValues<Real> compare;
+    for (unsigned n = 0; n < g_max_lm; ++n)
+      assert(compare(std::assoc_legendre(l, m, 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_lm; ++n) {
+      if (n & 1) // odd
+        doubleFactorials *= n;
+      else if (n != 0) // even and not zero
+        doubleFactorials /= n;
+
+      assert(compare(std::assoc_legendre(l, m, Real{0}), (n & 1) ? Real{0} : std::pow(-1, n / 2) * doubleFactorials));
+    }
+  }
+#endif
+}
+
+struct TestFloat {
+  template <class Real>
+  void operator()() {
+    test<Real>();
+  }
+};
+
+struct TestInt {
+  template <class Integer>
+  void operator()() {
+    // checks that std::assoc_legendre(unsigned, unsigned, Integer) actually wraps
+    // std::assoc_legendre(unsigned, unsigned, double)
+    for (unsigned l = 0; l < g_max_lm; ++l)
+      for (unsigned m = 0; m < g_max_lm; ++m)
+        for (Integer x : {-1, 0, 1})
+          assert(std::assoc_legendre(l, m, x) == std::assoc_legendre(l, m, 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());
+}