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

@@ -20,6 +20,12 @@
    template <class T, class ``__Policy``>
    ``__sf_result`` legendre_p_prime(int n, T x, const ``__Policy``&);
 
+   template <class T, class ``__Policy``>
+   T legendre_p_zeros(int l, int k, const ``__Policy``&);
+
+   template <class T>
+   T legendre_p_zeros(int l, int k);
+
    template <class T>
    ``__sf_result`` legendre_p(int n, int m, T x);
 
@@ -78,6 +84,26 @@ Legendre Polynomials:
 
 Returns the derivatives of the Legendre polynomials.
 
+    template <class T, class ``__Policy``>
+    T legendre_p_zeros(int l, int k, const ``__Policy``&);
+
+    template <class T>
+    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<double>(2*l+1, 0) == 0` holds.
+The rest of the zeros are returned in increasing order, so that `legendre_p_zeros<double>(l, k) < legendre_p_zeros<double>(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 <class T>
    ``__sf_result`` legendre_p(int l, int m, T x);
 

include/boost/math/special_functions/legendre.hpp

@@ -11,8 +11,10 @@
 #pragma once
 #endif
 
+#include <utility>
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/special_functions/factorials.hpp>
+#include <boost/math/tools/roots.hpp>
 #include <boost/math/tools/config.hpp>
 
 namespace boost{
@@ -69,10 +71,14 @@ T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false)
 }
 
 template <class T, class Policy>
-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 <class T, class Policy>
+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<T>(
+          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<T>());
+    if (k > half_n - 1)
+    {
+        return policies::raise_domain_error<T>(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<T>()*half<T>() - k)*pi<T>())/(n+half<T>());
+    T lower_bound = cos( (half_n - k)*pi<T>()/(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<Policy>();
+
+    auto f = [&] (T x) { T Pn;
+                         T Pn_prime = detail::legendre_p_prime_imp(n_in, x, pol, &Pn);
+                         return std::pair<T, T>(Pn, Pn_prime); };
+
+    const T x_nk = newton_raphson_iterate(f, x_nk_guess,
+                                          lower_bound, upper_bound,
+                                          policies::digits<T, Policy>(),
+                                          number_of_iterations);
+
+    if(number_of_iterations >= policies::get_max_root_iterations<Policy>())
+    {
+       return policies::raise_evaluation_error<T>(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 <class T, class Policy>
@@ -149,6 +239,21 @@ inline typename tools::promote_args<T>::type
    return boost::math::legendre_p_prime(l, x, policies::policy<>());
 }
 
+template <class T, class Policy>
+inline T legendre_p_zeros(int l, int k, const Policy& pol)
+{
+    if(l < 0)
+        return detail::legendre_p_zeros_imp<T>(-l-1, k, pol);
+
+    return detail::legendre_p_zeros_imp<T>(l, k, pol);
+}
+
+
+template <class T>
+inline T legendre_p_zeros(int l, int k)
+{
+   return boost::math::legendre_p_zeros<T>(l, k, policies::policy<>());
+}
 
 template <class T, class Policy>
 inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type

include/boost/math/special_functions/math_fwd.hpp

@@ -185,6 +185,13 @@ namespace boost
    typename tools::promote_args<T>::type
           legendre_p_prime(int l, T x);
 
+
+   template <class T, class Policy>
+   inline T legendre_p_zeros(int l, int k, const Policy& pol);
+
+   template <class T>
+   inline T legendre_p_zeros(int l, int k);
+
 #if !BOOST_WORKAROUND(BOOST_MSVC, <= 1310)
    template <class T, class Policy>
    typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type

test/test_legendre.cpp

@@ -222,4 +222,13 @@ BOOST_AUTO_TEST_CASE( test_main )
    test_legendre_p_prime<float>();
    test_legendre_p_prime<double>();
    test_legendre_p_prime<long double>();
+
+   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<float>();
+   test_legendre_p_zeros<double>();
+   test_legendre_p_zeros<long double>();
 }

test/test_legendre.hpp

@@ -15,6 +15,7 @@
 #include <boost/math/special_functions/math_fwd.hpp>
 #include <boost/math/special_functions/legendre.hpp>
 #include <boost/math/constants/constants.hpp>
+#include <boost/multiprecision/cpp_bin_float.hpp>
 #include <boost/array.hpp>
 #include "functor.hpp"
 
@@ -258,3 +259,81 @@ void test_legendre_p_prime()
         ++n;
     }
 }
+
+template<class Real>
+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<Real>::epsilon();
+
+    // Check the trivial cases:
+    BOOST_CHECK_CLOSE_FRACTION(legendre_p_zeros<Real>(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<Real>(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<Real>(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<Real>()); ++k)
+        {
+            Real next_zero = legendre_p_zeros<Real>(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<double>()); ++k)
+        {
+            double double_zero = legendre_p_zeros<double>(n, k);
+            cpp_bin_float_quad quad_zero = legendre_p_zeros<cpp_bin_float_quad>(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;
+}