From 6f3f92a89505def73aa696fd380c28b12cebee6a Mon Sep 17 00:00:00 2001 From: "Ryan M. Richard" Date: Mon, 31 Aug 2026 16:49:35 -0500 Subject: [PATCH] fix for macos --- test/CMakeLists.txt | 4 ++ test/assoc_legendre.cxx | 85 +++++++++++++++++++++++++++++++++++++++++ test/test_functions.hpp | 57 ++++++++++++++++++++++++++- 3 files changed, 144 insertions(+), 2 deletions(-) create mode 100644 test/assoc_legendre.cxx diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index e1a2460..ae6ff99 100644 --- a/test/CMakeLists.txt +++ b/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 ) diff --git a/test/assoc_legendre.cxx b/test/assoc_legendre.cxx new file mode 100644 index 0000000..dc61bed --- /dev/null +++ b/test/assoc_legendre.cxx @@ -0,0 +1,85 @@ +#include "catch2/catch_all.hpp" +#include +#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) ); + +} diff --git a/test/test_functions.hpp b/test/test_functions.hpp index 4a88729..6b4a824 100644 --- a/test/test_functions.hpp +++ b/test/test_functions.hpp @@ -1,8 +1,20 @@ #pragma once #include +#include #include +#include #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) { @@ -10,8 +22,49 @@ namespace detail { } } 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( std::cos(m*phi), std::sin(m*phi) ); } static auto evaluate(int l, int m, double x, double y, double z) {