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

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