ceph/src/boost/boost/math/special_functions/legendre.hpp

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