fix for macos

test/CMakeLists.txt

@@ -54,6 +54,9 @@ target_link_libraries( gausschebyshev2 PUBLIC integratorxx_common_ut )
 add_executable( composite_quadratures composite_quadratures.cxx )
 target_link_libraries( composite_quadratures PUBLIC integratorxx_common_ut )
 
+add_executable( assoc_legendre assoc_legendre.cxx )
+target_link_libraries( assoc_legendre PUBLIC integratorxx_common_ut )
+
 add_test( NAME QUADRATURE_MANIP         COMMAND quadrature_manipulation )
 add_test( NAME QUADRATURES_1D           COMMAND 1d_quadratures          )
 add_test( NAME QUADRATURES_COMPOSITE    COMMAND composite_quadratures   )
@@ -63,3 +66,4 @@ add_test( NAME QUADRATURES_LOBATTO      COMMAND gausslobatto            )
 add_test( NAME QUADRATURES_MURAKNOWLES  COMMAND muraknowles             )
 add_test( NAME QUADRATURES_CHEBYSHEV1   COMMAND gausschebyshev1         )
 add_test( NAME QUADRATURES_CHEBYSHEV2   COMMAND gausschebyshev2         )
+add_test( NAME ASSOC_LEGENDRE           COMMAND assoc_legendre          )

test/assoc_legendre.cxx

@@ -0,0 +1,85 @@
+#include "catch2/catch_all.hpp"
+#include <integratorxx/util/legendre.hpp>
+#include "test_functions.hpp"
+
+/* Parity check for the hand-rolled associated Legendre polynomial in
+ * test_functions.hpp.
+ *
+ * AssociatedLegendre::evaluate_fallback exists because libc++ does not
+ * implement the C++17 mathematical special functions ([sf.cmath]). Where the
+ * standard library DOES provide std::assoc_legendre, we check the fallback
+ * against it, so that the code path the libc++ platforms depend on is
+ * continuously verified rather than exercised only where it is the sole
+ * option.
+ */
+
+#ifdef IXX_TEST_HAS_STD_ASSOC_LEGENDRE
+
+TEST_CASE( "Associated Legendre matches std::assoc_legendre", "[assoc-legendre]" ) {
+
+  constexpr int    max_l   = 12;
+  constexpr int    npts    = 201;
+  // Both sides carry rounding error from their own recurrences; measured
+  // against an exact rational reference the fallback alone reaches ~1.2e-13
+  // relative near l=8. 1e-11 leaves room for the standard library's error and
+  // for compiler-to-compiler variation while staying orders of magnitude
+  // tighter than any defect would be.
+  constexpr double rel_tol = 1e-11;
+  constexpr double abs_tol = 1e-13;
+
+  for( int l = 0; l <= max_l; ++l )
+  for( int m = 0; m <= l;     ++m ) {
+
+    for( int i = 0; i < npts; ++i ) {
+      const double x   = -1.0 + 2.0 * i / (npts - 1.0);
+      const double ref = std::assoc_legendre(l, m, x);
+      const double val = AssociatedLegendre::evaluate_fallback(l, m, x);
+
+      const std::string msg = "P_l^m (L,M,X) = (" + std::to_string(l) + "," +
+        std::to_string(m) + "," + std::to_string(x) + ")";
+
+      // Relative agreement, with an absolute floor so values near a zero of the
+      // polynomial are not held to an unattainable relative tolerance
+      INFO( msg );
+      REQUIRE_THAT( val, Catch::Matchers::WithinRel(ref, rel_tol) ||
+                         Catch::Matchers::WithinAbs(ref, abs_tol) );
+    }
+
+  }
+
+}
+
+#endif
+
+TEST_CASE( "Associated Legendre special values", "[assoc-legendre]" ) {
+
+  // P_l^0 == P_l, checked against the recurrence in the library proper
+  for( int l = 0; l <= 12; ++l )
+  for( int i = 0; i < 51;  ++i ) {
+    const double x = -1.0 + 2.0 * i / 50.0;
+    double p_n;
+    std::tie(p_n, std::ignore, std::ignore) = IntegratorXX::eval_Pn(x, l);
+    INFO( "P_l^0 vs P_l, (L,X) = (" << l << "," << x << ")" );
+    REQUIRE_THAT( AssociatedLegendre::evaluate_fallback(l, 0, x),
+                  Catch::Matchers::WithinRel(p_n, 1e-11) ||
+                  Catch::Matchers::WithinAbs(p_n, 1e-13) );
+  }
+
+  // P_l^m(+/-1) == 0 for m > 0
+  for( int l = 1; l <= 12; ++l )
+  for( int m = 1; m <= l;  ++m ) {
+    REQUIRE( AssociatedLegendre::evaluate_fallback(l, m,  1.0) == 0.0 );
+    REQUIRE( AssociatedLegendre::evaluate_fallback(l, m, -1.0) == 0.0 );
+  }
+
+  // Closed forms for the low orders
+  const double x = 0.375;
+  const double s = std::sqrt(1.0 - x*x);
+  REQUIRE_THAT( AssociatedLegendre::evaluate_fallback(1, 1, x),
+                Catch::Matchers::WithinRel(s, 1e-14) );
+  REQUIRE_THAT( AssociatedLegendre::evaluate_fallback(2, 1, x),
+                Catch::Matchers::WithinRel(3.0 * x * s, 1e-14) );
+  REQUIRE_THAT( AssociatedLegendre::evaluate_fallback(2, 2, x),
+                Catch::Matchers::WithinRel(3.0 * (1.0 - x*x), 1e-14) );
+
+}

test/test_functions.hpp

@@ -1,17 +1,70 @@
 #pragma once
 #include <cmath>
+#include <complex>
 #include <random>
+#include <version>
 #include "quad_matcher.hpp"
 
+// The C++17 mathematical special functions ([sf.cmath]) are an optional part of
+// the standard library. libstdc++ and the MSVC STL provide them; libc++ does
+// not, so std::assoc_legendre is unavailable on macOS and on any clang build
+// against libc++. Detect it with the standard feature-test macro rather than by
+// sniffing the standard library, so this keeps working if libc++ implements
+// [sf.cmath] later.
+#if defined(__cpp_lib_math_special_functions) && \
+    __cpp_lib_math_special_functions >= 201603L
+#  define IXX_TEST_HAS_STD_ASSOC_LEGENDRE 1
+#endif
 
 namespace detail {
   constexpr size_t factorial(size_t n) {
     return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
   }
 }
 struct AssociatedLegendre {
+  /**
+   *  Associated Legendre polynomial P_l^m(x), evaluated with the standard
+   *  upward recurrence. This is compiled on every platform, whether or not
+   *  std::assoc_legendre exists, so that the parity test in assoc_legendre.cxx
+   *  can check it against the standard library wherever that is available.
+   *
+   *  Follows the [sf.cmath] convention, P_l^m(x) = (1-x^2)^(m/2) d^m/dx^m
+   *  P_l(x), i.e. WITHOUT the Condon-Shortley phase (-1)^m. This is the same
+   *  convention std::assoc_legendre uses, which is what makes the two
+   *  interchangeable.
+   */
+  static inline double evaluate_fallback(int l, int m, double x) {
+    // Preconditions of std::assoc_legendre; callers pass std::abs(m) and a cos
+    if( m < 0 or m > l or std::abs(x) > 1.0 ) return 0.0;
+
+    // P_m^m(x) = (2m-1)!! (1-x^2)^(m/2), accumulated as a product of the exact
+    // integer coefficients (2j-1) and sqrt(1-x^2) to avoid forming (2m-1)!! and
+    // the fractional power separately
+    const double s = std::sqrt( (1.0 - x) * (1.0 + x) );
+    double p_mm = 1.0;
+    for( int j = 1; j <= m; ++j ) p_mm *= (2 * j - 1) * s;
+    if( l == m ) return p_mm;
+
+    // P_{m+1}^m(x) = x (2m+1) P_m^m(x)
+    double p_lm_m1 = p_mm;
+    double p_lm    = x * (2 * m + 1) * p_mm;
+
+    // (l-m) P_l^m(x) = x (2l-1) P_{l-1}^m(x) - (l+m-1) P_{l-2}^m(x)
+    for( int ll = m + 2; ll <= l; ++ll ) {
+      const double p_lm_m2 = p_lm_m1;
+      p_lm_m1 = p_lm;
+      p_lm = ( x * (2 * ll - 1) * p_lm_m1 - (ll + m - 1) * p_lm_m2 ) / (ll - m);
+    }
+
+    return p_lm;
+  }
+
   static inline double evaluate(int l, int m, double x) {
+#ifdef IXX_TEST_HAS_STD_ASSOC_LEGENDRE
     return std::assoc_legendre(l, m, x);
+#else
+    return evaluate_fallback(l, m, x);
+#endif
   }
 };
 
@@ -27,8 +80,8 @@ struct SphericalHarmonic {
     }
 
     return prefactor * 
-           std::assoc_legendre(l, std::abs(m), std::cos(theta)) *
-           std::complex( std::cos(m*phi), std::sin(m*phi) );
+           AssociatedLegendre::evaluate(l, std::abs(m), std::cos(theta)) *
+           std::complex<double>( std::cos(m*phi), std::sin(m*phi) );
   }
 
   static auto evaluate(int l, int m, double x, double y, double z) {