math/special_functions/legendre.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SPECIAL_LEGENDRE_HPP
0007 #define BOOST_MATH_SPECIAL_LEGENDRE_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <utility>
0014 #include <vector>
0015 #include <type_traits>
0016 #include <boost/math/special_functions/math_fwd.hpp>
0017 #include <boost/math/special_functions/factorials.hpp>
0018 #include <boost/math/tools/roots.hpp>
0019 #include <boost/math/tools/config.hpp>
0020 #include <boost/math/tools/cxx03_warn.hpp>
0021
0022 namespace boost{
0023 namespace math{
0024
0025 // Recurrence relation for legendre P and Q polynomials:
0026 template <class T1, class T2, class T3>
0027 inline typename tools::promote_args<T1, T2, T3>::type
0028    legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1)
0029 {
0030    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
0031    return ((2 * l + 1) * result_type(x) * result_type(Pl) - l * result_type(Plm1)) / (l + 1);
0032 }
0033
0034 namespace detail{
0035
0036 // Implement Legendre P and Q polynomials via recurrence:
0037 template <class T, class Policy>
0038 T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false)
0039 {
0040    static const char* function = "boost::math::legrendre_p<%1%>(unsigned, %1%)";
0041    // Error handling:
0042    if((x < -1) || (x > 1))
0043       return policies::raise_domain_error<T>(
0044          function,
0045          "The Legendre Polynomial is defined for"
0046          " -1 <= x <= 1, but got x = %1%.", x, pol);
0047
0048    T p0, p1;
0049    if(second)
0050    {
0051       // A solution of the second kind (Q):
0052       p0 = (boost::math::log1p(x, pol) - boost::math::log1p(-x, pol)) / 2;
0053       p1 = x * p0 - 1;
0054    }
0055    else
0056    {
0057       // A solution of the first kind (P):
0058       p0 = 1;
0059       p1 = x;
0060    }
0061    if(l == 0)
0062       return p0;
0063
0064    unsigned n = 1;
0065
0066    while(n < l)
0067    {
0068       std::swap(p0, p1);
0069       p1 = static_cast<T>(boost::math::legendre_next(n, x, p0, p1));
0070       ++n;
0071    }
0072    return p1;
0073 }
0074
0075 template <class T, class Policy>
0076 T legendre_p_prime_imp(unsigned l, T x, const Policy& pol, T* Pn
0077 #ifdef BOOST_NO_CXX11_NULLPTR
0078    = 0
0079 #else
0080    = nullptr
0081 #endif
0082 )
0083 {
0084    static const char* function = "boost::math::legrendre_p_prime<%1%>(unsigned, %1%)";
0085    // Error handling:
0086    if ((x < -1) || (x > 1))
0087       return policies::raise_domain_error<T>(
0088          function,
0089          "The Legendre Polynomial is defined for"
0090          " -1 <= x <= 1, but got x = %1%.", x, pol);
0091
0092    if (l == 0)
0093     {
0094         if (Pn)
0095         {
0096            *Pn = 1;
0097         }
0098         return 0;
0099     }
0100     T p0 = 1;
0101     T p1 = x;
0102     T p_prime;
0103     bool odd = ((l & 1) == 1);
0104     // If the order is odd, we sum all the even polynomials:
0105     if (odd)
0106     {
0107         p_prime = p0;
0108     }
0109     else // Otherwise we sum the odd polynomials * (2n+1)
0110     {
0111         p_prime = 3*p1;
0112     }
0113
0114     unsigned n = 1;
0115     while(n < l - 1)
0116     {
0117        std::swap(p0, p1);
0118        p1 = static_cast<T>(boost::math::legendre_next(n, x, p0, p1));
0119        ++n;
0120        if (odd)
0121        {
0122           p_prime += (2*n+1)*p1;
0123           odd = false;
0124        }
0125        else
0126        {
0127            odd = true;
0128        }
0129     }
0130     // This allows us to evaluate the derivative and the function for the same cost.
0131     if (Pn)
0132     {
0133         std::swap(p0, p1);
0134         *Pn = static_cast<T>(boost::math::legendre_next(n, x, p0, p1));
0135     }
0136     return p_prime;
0137 }
0138
0139 template <class T, class Policy>
0140 struct legendre_p_zero_func
0141 {
0142    int n;
0143    const Policy& pol;
0144
0145    legendre_p_zero_func(int n_, const Policy& p) : n(n_), pol(p) {}
0146
0147    std::pair<T, T> operator()(T x) const
0148    {
0149       T Pn;
0150       T Pn_prime = detail::legendre_p_prime_imp(n, x, pol, &Pn);
0151       return std::pair<T, T>(Pn, Pn_prime);
0152    }
0153 };
0154
0155 template <class T, class Policy>
0156 std::vector<T> legendre_p_zeros_imp(int n, const Policy& pol)
0157 {
0158     using std::cos;
0159     using std::sin;
0160     using std::ceil;
0161     using std::sqrt;
0162     using boost::math::constants::pi;
0163     using boost::math::constants::half;
0164     using boost::math::tools::newton_raphson_iterate;
0165
0166     BOOST_MATH_ASSERT(n >= 0);
0167     std::vector<T> zeros;
0168     if (n == 0)
0169     {
0170         // There are no zeros of P_0(x) = 1.
0171         return zeros;
0172     }
0173     int k;
0174     if (n & 1)
0175     {
0176         zeros.resize((n-1)/2 + 1, std::numeric_limits<T>::quiet_NaN());
0177         zeros[0] = 0;
0178         k = 1;
0179     }
0180     else
0181     {
0182         zeros.resize(n/2, std::numeric_limits<T>::quiet_NaN());
0183         k = 0;
0184     }
0185     T half_n = ceil(n*half<T>());
0186
0187     while (k < (int)zeros.size())
0188     {
0189         // Bracket the root: Szego:
0190         // Gabriel Szego, Inequalities for the Zeros of Legendre Polynomials and Related Functions, Transactions of the American Mathematical Society, Vol. 39, No. 1 (1936)
0191         T theta_nk =  ((half_n - half<T>()*half<T>() - static_cast<T>(k))*pi<T>())/(static_cast<T>(n)+half<T>());
0192         T lower_bound = cos( (half_n - static_cast<T>(k))*pi<T>()/static_cast<T>(n + 1));
0193         T cos_nk = cos(theta_nk);
0194         T upper_bound = cos_nk;
0195         // First guess follows from:
0196         //  F. G. Tricomi, Sugli zeri dei polinomi sferici ed ultrasferici, Ann. Mat. Pura Appl., 31 (1950), pp. 93-97;
0197         T inv_n_sq = 1/static_cast<T>(n*n);
0198         T sin_nk = sin(theta_nk);
0199         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;
0200
0201         std::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>();
0202
0203         legendre_p_zero_func<T, Policy> f(n, pol);
0204
0205         const T x_nk = newton_raphson_iterate(f, x_nk_guess,
0206                                               lower_bound, upper_bound,
0207                                               policies::digits<T, Policy>(),
0208                                               number_of_iterations);
0209         if (number_of_iterations >= policies::get_max_root_iterations<Policy>())
0210         {
0211            policies::raise_evaluation_error<T>("legendre_p_zeros<%1%>", "Unable to locate solution in a reasonable time:"  // LCOV_EXCL_LINE
0212               " either there is no answer or the answer is infinite.  Current best guess is %1%", x_nk, Policy()); // LCOV_EXCL_LINE
0213         }
0214
0215         BOOST_MATH_ASSERT(lower_bound < x_nk);
0216         BOOST_MATH_ASSERT(upper_bound > x_nk);
0217         zeros[k] = x_nk;
0218         ++k;
0219     }
0220     return zeros;
0221 }
0222
0223 } // namespace detail
0224
0225 template <class T, class Policy>
0226 inline typename std::enable_if<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
0227    legendre_p(int l, T x, const Policy& pol)
0228 {
0229    typedef typename tools::promote_args<T>::type result_type;
0230    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0231    static const char* function = "boost::math::legendre_p<%1%>(unsigned, %1%)";
0232    if(l < 0)
0233       return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(-l-1, static_cast<value_type>(x), pol, false), function);
0234    return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, false), function);
0235 }
0236
0237
0238 template <class T, class Policy>
0239 inline typename std::enable_if<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
0240    legendre_p_prime(int l, T x, const Policy& pol)
0241 {
0242    typedef typename tools::promote_args<T>::type result_type;
0243    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0244    static const char* function = "boost::math::legendre_p_prime<%1%>(unsigned, %1%)";
0245    if(l < 0)
0246       return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(-l-1, static_cast<value_type>(x), pol), function);
0247    return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(l, static_cast<value_type>(x), pol), function);
0248 }
0249
0250 template <class T>
0251 inline typename tools::promote_args<T>::type
0252    legendre_p(int l, T x)
0253 {
0254    return boost::math::legendre_p(l, x, policies::policy<>());
0255 }
0256
0257 template <class T>
0258 inline typename tools::promote_args<T>::type
0259    legendre_p_prime(int l, T x)
0260 {
0261    return boost::math::legendre_p_prime(l, x, policies::policy<>());
0262 }
0263
0264 template <class T, class Policy>
0265 inline std::vector<T> legendre_p_zeros(int l, const Policy& pol)
0266 {
0267     if(l < 0)
0268         return detail::legendre_p_zeros_imp<T>(-l-1, pol);
0269
0270     return detail::legendre_p_zeros_imp<T>(l, pol);
0271 }
0272
0273
0274 template <class T>
0275 inline std::vector<T> legendre_p_zeros(int l)
0276 {
0277    return boost::math::legendre_p_zeros<T>(l, policies::policy<>());
0278 }
0279
0280 template <class T, class Policy>
0281 inline typename std::enable_if<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
0282    legendre_q(unsigned l, T x, const Policy& pol)
0283 {
0284    typedef typename tools::promote_args<T>::type result_type;
0285    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0286    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%)");
0287 }
0288
0289 template <class T>
0290 inline typename tools::promote_args<T>::type
0291    legendre_q(unsigned l, T x)
0292 {
0293    return boost::math::legendre_q(l, x, policies::policy<>());
0294 }
0295
0296 // Recurrence for associated polynomials:
0297 template <class T1, class T2, class T3>
0298 inline typename tools::promote_args<T1, T2, T3>::type
0299    legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1)
0300 {
0301    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
0302    return ((2 * l + 1) * result_type(x) * result_type(Pl) - (l + m) * result_type(Plm1)) / (l + 1 - m);
0303 }
0304
0305 namespace detail{
0306 // Legendre P associated polynomial:
0307 template <class T, class Policy>
0308 T legendre_p_imp(int l, int m, T x, T sin_theta_power, const Policy& pol)
0309 {
0310    BOOST_MATH_STD_USING
0311    // Error handling:
0312    if((x < -1) || (x > 1))
0313       return policies::raise_domain_error<T>(
0314       "boost::math::legendre_p<%1%>(int, int, %1%)",
0315          "The associated Legendre Polynomial is defined for"
0316          " -1 <= x <= 1, but got x = %1%.", x, pol);
0317    // Handle negative arguments first:
0318    if(l < 0)
0319       return legendre_p_imp(-l-1, m, x, sin_theta_power, pol);
0320    if ((l == 0) && (m == -1))
0321    {
0322       return sqrt((1 - x) / (1 + x));
0323    }
0324    if ((l == 1) && (m == 0))
0325    {
0326       return x;
0327    }
0328    if (-m == l)
0329    {
0330       return pow((1 - x * x) / 4, T(l) / 2) / boost::math::tgamma<T>(l + 1, pol);
0331    }
0332    if(m < 0)
0333    {
0334       int sign = (m&1) ? -1 : 1;
0335       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);
0336    }
0337    // Special cases:
0338    if(m > l)
0339       return 0;
0340    if(m == 0)
0341       return boost::math::legendre_p(l, x, pol);
0342
0343    T p0 = boost::math::double_factorial<T>(2 * m - 1, pol) * sin_theta_power;
0344
0345    if(m&1)
0346       p0 *= -1;
0347    if(m == l)
0348       return p0;
0349
0350    T p1 = x * (2 * m + 1) * p0;
0351
0352    int n = m + 1;
0353
0354    while(n < l)
0355    {
0356       std::swap(p0, p1);
0357       p1 = boost::math::legendre_next(n, m, x, p0, p1);
0358       ++n;
0359    }
0360    return p1;
0361 }
0362
0363 template <class T, class Policy>
0364 inline T legendre_p_imp(int l, int m, T x, const Policy& pol)
0365 {
0366    BOOST_MATH_STD_USING
0367    // TODO: we really could use that mythical "pow1p" function here:
0368    return legendre_p_imp(l, m, x, static_cast<T>(pow(1 - x*x, T(abs(m))/2)), pol);
0369 }
0370
0371 }
0372
0373 template <class T, class Policy>
0374 inline typename tools::promote_args<T>::type
0375    legendre_p(int l, int m, T x, const Policy& pol)
0376 {
0377    typedef typename tools::promote_args<T>::type result_type;
0378    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0379    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%)");
0380 }
0381
0382 template <class T>
0383 inline typename tools::promote_args<T>::type
0384    legendre_p(int l, int m, T x)
0385 {
0386    return boost::math::legendre_p(l, m, x, policies::policy<>());
0387 }
0388
0389 } // namespace math
0390 } // namespace boost
0391
0392 #endif // BOOST_MATH_SPECIAL_LEGENDRE_HPP