From 4110a69416b7919412d06d87c954696400a67f29 Mon Sep 17 00:00:00 2001 From: Nick Thompson Date: Fri, 14 Apr 2017 15:51:39 -0500 Subject: [PATCH] Zeros of Legendre polynomials. This uses a root bracketing given by Szego with an asymptotic by Tricomi to get a domain and an initial guess for the root, then refines it via Newton's method. --- doc/sf/legendre.qbk | 26 +++++ .../boost/math/special_functions/legendre.hpp | 107 +++++++++++++++++- .../boost/math/special_functions/math_fwd.hpp | 7 ++ test/test_legendre.cpp | 9 ++ test/test_legendre.hpp | 79 +++++++++++++ 5 files changed, 227 insertions(+), 1 deletion(-) diff --git a/doc/sf/legendre.qbk b/doc/sf/legendre.qbk index 7f09bd7d0..83be46f8f 100644 --- a/doc/sf/legendre.qbk +++ b/doc/sf/legendre.qbk @@ -20,6 +20,12 @@ template ``__sf_result`` legendre_p_prime(int n, T x, const ``__Policy``&); + template + T legendre_p_zeros(int l, int k, const ``__Policy``&); + + template + T legendre_p_zeros(int l, int k); + template ``__sf_result`` legendre_p(int n, int m, T x); @@ -78,6 +84,26 @@ Legendre Polynomials: Returns the derivatives of the Legendre polynomials. + template + T legendre_p_zeros(int l, int k, const ``__Policy``&); + + template + T legendre_p_zeros(int l, int k); + +The zeros of the Legendre polynomials are calculated by Newton's method using an initial guess given by Tricomi with root bracketing provided by Szego. +Although the use of Newton's method is believed to be less robust than (say) the Golub & Welsch algorithm, +all zeros are computed to 1 ULP up to `l = 100`, and 2 ULP up to `l = 350` in double precision. + +Since the Legendre polynomials are alternatively even and odd, only the non-negative zeros are returned. +For the odd Legendre polynomials, the identity `legendre_p_zeros(2*l+1, 0) == 0` holds. +The rest of the zeros are returned in increasing order, so that `legendre_p_zeros(l, k) < legendre_p_zeros(l, k+1)`. +For each `l`, `k` may take the values 0 to `ceil(l*0.5)`, as the other zeros are symmetric about zero. + +Note that the arguments to the routine are integers, and the output is a real-type. Hence the template argument is mandatory. +The time to extract a single root is linear in `l` (this is scaling to evaluate the Legendre polynomials), so recovering all roots is [bigo](`l`[super 2]). +Algorithms with better linear scaling [@ http://dx.doi.org/10.1137/06067016X exist] for recovering all roots, but they are complicated. +This implementation proceeds under the assumption that calculating zeros of these functions will not be a bottleneck for any workflow. + template ``__sf_result`` legendre_p(int l, int m, T x); diff --git a/include/boost/math/special_functions/legendre.hpp b/include/boost/math/special_functions/legendre.hpp index 8b3257313..a685de42d 100644 --- a/include/boost/math/special_functions/legendre.hpp +++ b/include/boost/math/special_functions/legendre.hpp @@ -11,8 +11,10 @@ #pragma once #endif +#include #include #include +#include #include namespace boost{ @@ -69,10 +71,14 @@ T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false) } template -T legendre_p_prime_imp(unsigned l, T x, const Policy& pol) +T legendre_p_prime_imp(unsigned l, T x, const Policy& pol, T* Pn = nullptr) { if (l == 0) { + if (Pn) + { + *Pn = 1; + } return 0; } T p0 = 1; @@ -105,9 +111,93 @@ T legendre_p_prime_imp(unsigned l, T x, const Policy& pol) odd = true; } } + // This allows us to evaluate the derivative and the function for the same cost. + if (Pn) + { + std::swap(p0, p1); + *Pn = boost::math::legendre_next(n, x, p0, p1); + } return p_prime; } +template +T legendre_p_zeros_imp(int n_in, int k_in, const Policy& pol) +{ + using std::cos; + using std::sin; + using std::floor; + using std::ceil; + using std::sqrt; + using boost::math::constants::pi; + using boost::math::constants::half; + using boost::math::tools::newton_raphson_iterate; + + static const char* function = "boost::math::legrendre_p_zeros<%1%>(int, int)"; + if (n_in == 0) + { + return policies::raise_domain_error( + function, + "The %1%-th Legendre polynomial has no zeros.\n", n_in, pol); + } + // If n is odd, and k is zero, then the root is zero. + // Taking this branch as a special case is reasonable, because + // the root bracketing gives us an interval [0, 0] with floating point rounding + // error on both positive and negative. This can lead to weird behavior, such as + // upper_bound < lower_bound. + if ( (n_in & 1) && k_in == 0) + { + return (T) 0; + } + + if (n_in == 2 && k_in == 0) + { + return (T) 1/sqrt(3); + } + + T n = n_in; + T k = k_in; + T half_n = ceil(n*half()); + if (k > half_n - 1) + { + return policies::raise_domain_error(function, + "k must be in range {0, 1, ..., ceil(n/2) - 1}," + " requested root %1% of a Legendre polynomial.\n", + k, pol); + } + + // Bracket the root: Szego: + // Gabriel Szego, Inequalities for the Zeros of Legendre Polynomials and Related Functions, Transactions of the American Mathematical Society, Vol. 39, No. 1 (1936) + T theta_nk = ((half_n - half()*half() - k)*pi())/(n+half()); + T lower_bound = cos( (half_n - k)*pi()/(n + 1)); + T cos_nk = cos(theta_nk); + T upper_bound = cos_nk; + // First guess follows from: + // F. G. Tricomi, Sugli zeri dei polinomi sferici ed ultrasferici, Ann. Mat. Pura Appl., 31 (1950), pp. 93–97; + T inv_n_sq = (T) 1/(n*n); + T sin_nk = sin(theta_nk); + T x_nk_guess = (1 - inv_n_sq/8 + inv_n_sq /(8*n) - (inv_n_sq*inv_n_sq/384)*(39 - 28 / (sin_nk*sin_nk) ) )*cos_nk; + + boost::uintmax_t number_of_iterations = policies::get_max_root_iterations(); + + auto f = [&] (T x) { T Pn; + T Pn_prime = detail::legendre_p_prime_imp(n_in, x, pol, &Pn); + return std::pair(Pn, Pn_prime); }; + + const T x_nk = newton_raphson_iterate(f, x_nk_guess, + lower_bound, upper_bound, + policies::digits(), + number_of_iterations); + + if(number_of_iterations >= policies::get_max_root_iterations()) + { + return policies::raise_evaluation_error(function, "Unable to locate root in a reasonable time:" + " Current best guess is %1%", x_nk, Policy()); + } + BOOST_ASSERT(lower_bound < x_nk); + BOOST_ASSERT(upper_bound > x_nk); + return x_nk; +} + } // namespace detail template @@ -149,6 +239,21 @@ inline typename tools::promote_args::type return boost::math::legendre_p_prime(l, x, policies::policy<>()); } +template +inline T legendre_p_zeros(int l, int k, const Policy& pol) +{ + if(l < 0) + return detail::legendre_p_zeros_imp(-l-1, k, pol); + + return detail::legendre_p_zeros_imp(l, k, pol); +} + + +template +inline T legendre_p_zeros(int l, int k) +{ + return boost::math::legendre_p_zeros(l, k, policies::policy<>()); +} template inline typename boost::enable_if_c::value, typename tools::promote_args::type>::type diff --git a/include/boost/math/special_functions/math_fwd.hpp b/include/boost/math/special_functions/math_fwd.hpp index cd01b88ef..2fccc1656 100644 --- a/include/boost/math/special_functions/math_fwd.hpp +++ b/include/boost/math/special_functions/math_fwd.hpp @@ -185,6 +185,13 @@ namespace boost typename tools::promote_args::type legendre_p_prime(int l, T x); + + template + inline T legendre_p_zeros(int l, int k, const Policy& pol); + + template + inline T legendre_p_zeros(int l, int k); + #if !BOOST_WORKAROUND(BOOST_MSVC, <= 1310) template typename boost::enable_if_c::value, typename tools::promote_args::type>::type diff --git a/test/test_legendre.cpp b/test/test_legendre.cpp index e5fd3e3a4..9053d887d 100644 --- a/test/test_legendre.cpp +++ b/test/test_legendre.cpp @@ -222,4 +222,13 @@ BOOST_AUTO_TEST_CASE( test_main ) test_legendre_p_prime(); test_legendre_p_prime(); test_legendre_p_prime(); + + int distance = test_legendre_p_zeros_double_ulp(1, 107); + BOOST_CHECK(distance <= 1); + // This test is very expensive; the total runtime grows cubically with n. + //distance = test_legendre_p_zeros_double_ulp(108, 350); + //BOOST_CHECK(distance <= 2); + test_legendre_p_zeros(); + test_legendre_p_zeros(); + test_legendre_p_zeros(); } diff --git a/test/test_legendre.hpp b/test/test_legendre.hpp index 81593ff16..447664cfa 100644 --- a/test/test_legendre.hpp +++ b/test/test_legendre.hpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include "functor.hpp" @@ -258,3 +259,81 @@ void test_legendre_p_prime() ++n; } } + +template +void test_legendre_p_zeros() +{ + using std::sqrt; + using std::abs; + using boost::math::legendre_p_zeros; + using boost::math::legendre_p; + Real tol = std::numeric_limits::epsilon(); + + // Check the trivial cases: + BOOST_CHECK_CLOSE_FRACTION(legendre_p_zeros(2, 0), (Real) 1/ sqrt(3), tol); + + // The first zero of the odd Legendre Polynomials is obviously zero. + for (int n = 1; n < 5000; n += 2) + { + BOOST_CHECK_SMALL(legendre_p_zeros(n, 0), tol); + } + + // Don't take the tolerances too seriously. + // The other test shows that the zeros are estimated more accurately than the function! + for (int n = 1; n < 130; ++n) + { + Real previous_zero = legendre_p_zeros(n, 0); + if (!(n & 1)) + { + // Zero is not a zero of the odd Legendre polynomials + BOOST_CHECK(previous_zero > 0); + BOOST_CHECK_SMALL(legendre_p(n, previous_zero), 550*tol); + } + for (int k = 1; k < ceil(n*boost::math::constants::half()); ++k) + { + Real next_zero = legendre_p_zeros(n, k); + BOOST_CHECK(next_zero > previous_zero); + + std::string err = "Tolerance failed for (n, k) = (" + std::to_string(n) + "," + std::to_string(k) + ")\n"; + if (n < 40) + { + BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 100*tol, + err); + } + else + { + BOOST_CHECK_MESSAGE( abs(legendre_p(n, next_zero)) < 1000*tol, + err); + } + previous_zero = next_zero; + } + // The zeros of orthogonal polynomials are contained strictly in (a, b). + BOOST_CHECK(previous_zero < 1); + } + return; +} + +int test_legendre_p_zeros_double_ulp(int min_x, int max_n) +{ + using std::abs; + using boost::math::legendre_p_zeros; + using boost::math::float_distance; + using boost::multiprecision::cpp_bin_float_quad; + + double max_float_distance = 0; + for (int n = min_x; n < max_n; ++n) + { + for (int k = 1; k < ceil(n*boost::math::constants::half()); ++k) + { + double double_zero = legendre_p_zeros(n, k); + cpp_bin_float_quad quad_zero = legendre_p_zeros(n, k); + double d = abs(float_distance(double_zero, (double) quad_zero)); + if (d > max_float_distance) + { + max_float_distance = d; + } + } + } + + return (int) max_float_distance; +}