381 lines · plain
1// (C) Copyright John Maddock 2006.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_SPECIAL_LEGENDRE_HPP7#define BOOST_MATH_SPECIAL_LEGENDRE_HPP8 9#ifdef _MSC_VER10#pragma once11#endif12 13#include <utility>14#include <vector>15#include <type_traits>16#include <boost/math/special_functions/math_fwd.hpp>17#include <boost/math/special_functions/factorials.hpp>18#include <boost/math/tools/roots.hpp>19#include <boost/math/tools/config.hpp>20#include <boost/math/tools/cxx03_warn.hpp>21 22namespace boost{23namespace math{24 25// Recurrence relation for legendre P and Q polynomials:26template <class T1, class T2, class T3>27inline typename tools::promote_args<T1, T2, T3>::type28 legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1)29{30 typedef typename tools::promote_args<T1, T2, T3>::type result_type;31 return ((2 * l + 1) * result_type(x) * result_type(Pl) - l * result_type(Plm1)) / (l + 1);32}33 34namespace detail{35 36// Implement Legendre P and Q polynomials via recurrence:37template <class T, class Policy>38T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false)39{40 static const char* function = "boost::math::legrendre_p<%1%>(unsigned, %1%)";41 // Error handling:42 if((x < -1) || (x > 1))43 return policies::raise_domain_error<T>(function, "The Legendre Polynomial is defined for -1 <= x <= 1, but got x = %1%.", x, pol);44 45 T p0, p1;46 if(second)47 {48 // A solution of the second kind (Q):49 p0 = (boost::math::log1p(x, pol) - boost::math::log1p(-x, pol)) / 2;50 p1 = x * p0 - 1;51 }52 else53 {54 // A solution of the first kind (P):55 p0 = 1;56 p1 = x;57 }58 if(l == 0)59 return p0;60 61 unsigned n = 1;62 63 while(n < l)64 {65 std::swap(p0, p1);66 p1 = static_cast<T>(boost::math::legendre_next(n, x, p0, p1));67 ++n;68 }69 return p1;70}71 72template <class T, class Policy>73T legendre_p_prime_imp(unsigned l, T x, const Policy& pol, T* Pn 74#ifdef BOOST_NO_CXX11_NULLPTR75 = 076#else77 = nullptr78#endif79)80{81 static const char* function = "boost::math::legrendre_p_prime<%1%>(unsigned, %1%)";82 // Error handling:83 if ((x < -1) || (x > 1))84 return policies::raise_domain_error<T>(function, "The Legendre Polynomial is defined for -1 <= x <= 1, but got x = %1%.", x, pol);85 86 if (l == 0)87 {88 BOOST_MATH_ASSERT(Pn == nullptr); // There are no zeros of P_0 so we shoud never call this with l = 0 and Pn non-null.89 return 0;90 }91 T p0 = 1;92 T p1 = x;93 T p_prime;94 bool odd = ((l & 1) == 1);95 // If the order is odd, we sum all the even polynomials:96 if (odd)97 {98 p_prime = p0;99 }100 else // Otherwise we sum the odd polynomials * (2n+1)101 {102 p_prime = 3*p1;103 }104 105 unsigned n = 1;106 while(n < l - 1)107 {108 std::swap(p0, p1);109 p1 = static_cast<T>(boost::math::legendre_next(n, x, p0, p1));110 ++n;111 if (odd)112 {113 p_prime += (2*n+1)*p1;114 odd = false;115 }116 else117 {118 odd = true;119 }120 }121 // This allows us to evaluate the derivative and the function for the same cost.122 if (Pn)123 {124 std::swap(p0, p1);125 *Pn = static_cast<T>(boost::math::legendre_next(n, x, p0, p1));126 }127 return p_prime;128}129 130template <class T, class Policy>131struct legendre_p_zero_func132{133 int n;134 const Policy& pol;135 136 legendre_p_zero_func(int n_, const Policy& p) : n(n_), pol(p) {}137 138 std::pair<T, T> operator()(T x) const139 { 140 T Pn;141 T Pn_prime = detail::legendre_p_prime_imp(n, x, pol, &Pn);142 return std::pair<T, T>(Pn, Pn_prime); 143 }144};145 146template <class T, class Policy>147std::vector<T> legendre_p_zeros_imp(int n, const Policy& pol)148{149 using std::cos;150 using std::sin;151 using std::ceil;152 using std::sqrt;153 using boost::math::constants::pi;154 using boost::math::constants::half;155 using boost::math::tools::newton_raphson_iterate;156 157 BOOST_MATH_ASSERT(n >= 0);158 std::vector<T> zeros;159 if (n == 0)160 {161 // There are no zeros of P_0(x) = 1.162 return zeros;163 }164 int k;165 if (n & 1)166 {167 zeros.resize((n-1)/2 + 1, std::numeric_limits<T>::quiet_NaN());168 zeros[0] = 0;169 k = 1;170 }171 else172 {173 zeros.resize(n/2, std::numeric_limits<T>::quiet_NaN());174 k = 0;175 }176 T half_n = ceil(n*half<T>());177 178 while (k < (int)zeros.size())179 {180 // Bracket the root: Szego:181 // Gabriel Szego, Inequalities for the Zeros of Legendre Polynomials and Related Functions, Transactions of the American Mathematical Society, Vol. 39, No. 1 (1936)182 T theta_nk = ((half_n - half<T>()*half<T>() - static_cast<T>(k))*pi<T>())/(static_cast<T>(n)+half<T>());183 T lower_bound = cos( (half_n - static_cast<T>(k))*pi<T>()/static_cast<T>(n + 1));184 T cos_nk = cos(theta_nk);185 T upper_bound = cos_nk;186 // First guess follows from:187 // F. G. Tricomi, Sugli zeri dei polinomi sferici ed ultrasferici, Ann. Mat. Pura Appl., 31 (1950), pp. 93-97;188 T inv_n_sq = 1/static_cast<T>(n*n);189 T sin_nk = sin(theta_nk);190 T x_nk_guess = (1 - inv_n_sq/static_cast<T>(8) + inv_n_sq /static_cast<T>(8*n) - (inv_n_sq*inv_n_sq/384)*(39 - 28 / (sin_nk*sin_nk) ) )*cos_nk;191 192 std::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();193 194 legendre_p_zero_func<T, Policy> f(n, pol);195 196 const T x_nk = newton_raphson_iterate(f, x_nk_guess,197 lower_bound, upper_bound,198 policies::digits<T, Policy>(),199 number_of_iterations);200 if (number_of_iterations >= policies::get_max_root_iterations<Policy>())201 {202 policies::raise_evaluation_error<T>("legendre_p_zeros<%1%>", "Unable to locate solution in a reasonable time:" // LCOV_EXCL_LINE203 " either there is no answer or the answer is infinite. Current best guess is %1%", x_nk, Policy()); // LCOV_EXCL_LINE204 }205 206 BOOST_MATH_ASSERT(lower_bound < x_nk);207 BOOST_MATH_ASSERT(upper_bound > x_nk);208 zeros[k] = x_nk;209 ++k;210 }211 return zeros;212} // LCOV_EXCL_LINE213 214} // namespace detail215 216template <class T, class Policy>217inline typename std::enable_if<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type218 legendre_p(int l, T x, const Policy& pol)219{220 typedef typename tools::promote_args<T>::type result_type;221 typedef typename policies::evaluation<result_type, Policy>::type value_type;222 static const char* function = "boost::math::legendre_p<%1%>(unsigned, %1%)";223 if(l < 0)224 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(-l-1, static_cast<value_type>(x), pol, false), function);225 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, false), function);226}227 228 229template <class T, class Policy>230inline typename std::enable_if<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type231 legendre_p_prime(int l, T x, const Policy& pol)232{233 typedef typename tools::promote_args<T>::type result_type;234 typedef typename policies::evaluation<result_type, Policy>::type value_type;235 static const char* function = "boost::math::legendre_p_prime<%1%>(unsigned, %1%)";236 if(l < 0)237 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(-l-1, static_cast<value_type>(x), pol), function);238 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(l, static_cast<value_type>(x), pol), function);239}240 241template <class T>242inline typename tools::promote_args<T>::type243 legendre_p(int l, T x)244{245 return boost::math::legendre_p(l, x, policies::policy<>());246}247 248template <class T>249inline typename tools::promote_args<T>::type250 legendre_p_prime(int l, T x)251{252 return boost::math::legendre_p_prime(l, x, policies::policy<>());253}254 255template <class T, class Policy>256inline std::vector<T> legendre_p_zeros(int l, const Policy& pol)257{258 if(l < 0)259 return detail::legendre_p_zeros_imp<T>(-l-1, pol);260 261 return detail::legendre_p_zeros_imp<T>(l, pol);262}263 264 265template <class T>266inline std::vector<T> legendre_p_zeros(int l)267{268 return boost::math::legendre_p_zeros<T>(l, policies::policy<>());269}270 271template <class T, class Policy>272inline typename std::enable_if<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type273 legendre_q(unsigned l, T x, const Policy& pol)274{275 typedef typename tools::promote_args<T>::type result_type;276 typedef typename policies::evaluation<result_type, Policy>::type value_type;277 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, true), "boost::math::legendre_q<%1%>(unsigned, %1%)");278}279 280template <class T>281inline typename tools::promote_args<T>::type282 legendre_q(unsigned l, T x)283{284 return boost::math::legendre_q(l, x, policies::policy<>());285}286 287// Recurrence for associated polynomials:288template <class T1, class T2, class T3>289inline typename tools::promote_args<T1, T2, T3>::type290 legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1)291{292 typedef typename tools::promote_args<T1, T2, T3>::type result_type;293 return ((2 * l + 1) * result_type(x) * result_type(Pl) - (l + m) * result_type(Plm1)) / (l + 1 - m);294}295 296namespace detail{297// Legendre P associated polynomial:298template <class T, class Policy>299T legendre_p_imp(int l, int m, T x, T sin_theta_power, const Policy& pol)300{301 BOOST_MATH_STD_USING302 // Error handling:303 if((x < -1) || (x > 1))304 return policies::raise_domain_error<T>("boost::math::legendre_p<%1%>(int, int, %1%)", "The associated Legendre Polynomial is defined for -1 <= x <= 1, but got x = %1%.", x, pol);305 // Handle negative arguments first:306 if(l < 0)307 return legendre_p_imp(-l-1, m, x, sin_theta_power, pol);308 if ((l == 0) && (m == -1))309 {310 return sqrt((1 - x) / (1 + x));311 }312 if ((l == 1) && (m == 0))313 {314 return x;315 }316 if (-m == l)317 {318 return pow((1 - x * x) / 4, T(l) / 2) / boost::math::tgamma<T>(l + 1, pol);319 }320 if(m < 0)321 {322 int sign = (m&1) ? -1 : 1;323 return sign * boost::math::tgamma_ratio(static_cast<T>(l+m+1), static_cast<T>(l+1-m), pol) * legendre_p_imp(l, -m, x, sin_theta_power, pol);324 }325 // Special cases:326 if(m > l)327 return 0;328 if(m == 0)329 return boost::math::legendre_p(l, x, pol);330 331 T p0 = boost::math::double_factorial<T>(2 * m - 1, pol) * sin_theta_power;332 333 if(m&1)334 p0 *= -1;335 if(m == l)336 return p0;337 338 T p1 = x * (2 * m + 1) * p0;339 340 int n = m + 1;341 342 while(n < l)343 {344 std::swap(p0, p1);345 p1 = boost::math::legendre_next(n, m, x, p0, p1);346 ++n;347 }348 return p1;349}350 351template <class T, class Policy>352inline T legendre_p_imp(int l, int m, T x, const Policy& pol)353{354 BOOST_MATH_STD_USING355 // TODO: we really could use that mythical "pow1p" function here:356 return legendre_p_imp(l, m, x, static_cast<T>(pow(1 - x*x, T(abs(m))/2)), pol);357}358 359}360 361template <class T, class Policy>362inline typename tools::promote_args<T>::type363 legendre_p(int l, int m, T x, const Policy& pol)364{365 typedef typename tools::promote_args<T>::type result_type;366 typedef typename policies::evaluation<result_type, Policy>::type value_type;367 return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_imp(l, m, static_cast<value_type>(x), pol), "boost::math::legendre_p<%1%>(int, int, %1%)");368}369 370template <class T>371inline typename tools::promote_args<T>::type372 legendre_p(int l, int m, T x)373{374 return boost::math::legendre_p(l, m, x, policies::policy<>());375}376 377} // namespace math378} // namespace boost379 380#endif // BOOST_MATH_SPECIAL_LEGENDRE_HPP381