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 <boost/math/special_functions/math_fwd.hpp>
  15 #include <boost/math/special_functions/factorials.hpp>
  16 #include <boost/math/tools/config.hpp>
  17
  18 namespace boost{
  19 namespace math{
  20
  21 // Recurrance relation for legendre P and Q polynomials:
  22 template <class T1, class T2, class T3>
  23 inline typename tools::promote_args<T1, T2, T3>::type
  24    legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1)
  25 {
  26    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
  27    return ((2 * l + 1) * result_type(x) * result_type(Pl) - l * result_type(Plm1)) / (l + 1);
  28 }
  29
  30 namespace detail{
  31
  32 // Implement Legendre P and Q polynomials via recurrance:
  33 template <class T, class Policy>
  34 T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false)
  35 {
  36    static const char* function = "boost::math::legrendre_p<%1%>(unsigned, %1%)";
  37    // Error handling:
  38    if((x < -1) || (x > 1))
  39       return policies::raise_domain_error<T>(
  40          function,
  41          "The Legendre Polynomial is defined for"
  42          " -1 <= x <= 1, but got x = %1%.", x, pol);
  43
  44    T p0, p1;
  45    if(second)
  46    {
  47       // A solution of the second kind (Q):
  48       p0 = (boost::math::log1p(x, pol) - boost::math::log1p(-x, pol)) / 2;
  49       p1 = x * p0 - 1;
  50    }
  51    else
  52    {
  53       // A solution of the first kind (P):
  54       p0 = 1;
  55       p1 = x;
  56    }
  57    if(l == 0)
  58       return p0;
  59
  60    unsigned n = 1;
  61
  62    while(n < l)
  63    {
  64       std::swap(p0, p1);
  65       p1 = boost::math::legendre_next(n, x, p0, p1);
  66       ++n;
  67    }
  68    return p1;
  69 }
  70
  71 } // namespace detail
  72
  73 template <class T, class Policy>
  74 inline typename tools::promote_args<T>::type
  75    legendre_p(int l, T x, const Policy& pol)
  76 {
  77    typedef typename tools::promote_args<T>::type result_type;
  78    typedef typename policies::evaluation<result_type, Policy>::type value_type;
  79    static const char* function = "boost::math::legendre_p<%1%>(unsigned, %1%)";
  80    if(l < 0)
  81       return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(-l-1, static_cast<value_type>(x), pol, false), function);
  82    return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, false), function);
  83 }
  84
  85 template <class T>
  86 inline typename tools::promote_args<T>::type
  87    legendre_p(int l, T x)
  88 {
  89    return boost::math::legendre_p(l, x, policies::policy<>());
  90 }
  91
  92 template <class T, class Policy>
  93 inline typename tools::promote_args<T>::type
  94    legendre_q(unsigned l, T x, const Policy& pol)
  95 {
  96    typedef typename tools::promote_args<T>::type result_type;
  97    typedef typename policies::evaluation<result_type, Policy>::type value_type;
  98    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%)");
  99 }
 100
 101 template <class T>
 102 inline typename tools::promote_args<T>::type
 103    legendre_q(unsigned l, T x)
 104 {
 105    return boost::math::legendre_q(l, x, policies::policy<>());
 106 }
 107
 108 // Recurrence for associated polynomials:
 109 template <class T1, class T2, class T3>
 110 inline typename tools::promote_args<T1, T2, T3>::type
 111    legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1)
 112 {
 113    typedef typename tools::promote_args<T1, T2, T3>::type result_type;
 114    return ((2 * l + 1) * result_type(x) * result_type(Pl) - (l + m) * result_type(Plm1)) / (l + 1 - m);
 115 }
 116
 117 namespace detail{
 118 // Legendre P associated polynomial:
 119 template <class T, class Policy>
 120 T legendre_p_imp(int l, int m, T x, T sin_theta_power, const Policy& pol)
 121 {
 122    // Error handling:
 123    if((x < -1) || (x > 1))
 124       return policies::raise_domain_error<T>(
 125       "boost::math::legendre_p<%1%>(int, int, %1%)",
 126          "The associated Legendre Polynomial is defined for"
 127          " -1 <= x <= 1, but got x = %1%.", x, pol);
 128    // Handle negative arguments first:
 129    if(l < 0)
 130       return legendre_p_imp(-l-1, m, x, sin_theta_power, pol);
 131    if(m < 0)
 132    {
 133       int sign = (m&1) ? -1 : 1;
 134       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);
 135    }
 136    // Special cases:
 137    if(m > l)
 138       return 0;
 139    if(m == 0)
 140       return boost::math::legendre_p(l, x, pol);
 141
 142    T p0 = boost::math::double_factorial<T>(2 * m - 1, pol) * sin_theta_power;
 143
 144    if(m&1)
 145       p0 *= -1;
 146    if(m == l)
 147       return p0;
 148
 149    T p1 = x * (2 * m + 1) * p0;
 150
 151    int n = m + 1;
 152
 153    while(n < l)
 154    {
 155       std::swap(p0, p1);
 156       p1 = boost::math::legendre_next(n, m, x, p0, p1);
 157       ++n;
 158    }
 159    return p1;
 160 }
 161
 162 template <class T, class Policy>
 163 inline T legendre_p_imp(int l, int m, T x, const Policy& pol)
 164 {
 165    BOOST_MATH_STD_USING
 166    // TODO: we really could use that mythical "pow1p" function here:
 167    return legendre_p_imp(l, m, x, static_cast<T>(pow(1 - x*x, T(abs(m))/2)), pol);
 168 }
 169
 170 }
 171
 172 template <class T, class Policy>
 173 inline typename tools::promote_args<T>::type
 174    legendre_p(int l, int m, T x, const Policy& pol)
 175 {
 176    typedef typename tools::promote_args<T>::type result_type;
 177    typedef typename policies::evaluation<result_type, Policy>::type value_type;
 178    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%)");
 179 }
 180
 181 template <class T>
 182 inline typename tools::promote_args<T>::type
 183    legendre_p(int l, int m, T x)
 184 {
 185    return boost::math::legendre_p(l, m, x, policies::policy<>());
 186 }
 187
 188 } // namespace math
 189 } // namespace boost
 190
 191 #endif // BOOST_MATH_SPECIAL_LEGENDRE_HPP
 192
 193
 194