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)
7 #ifndef BOOST_MATH_SPECIAL_LEGENDRE_HPP
8 #define BOOST_MATH_SPECIAL_LEGENDRE_HPP
14 #include <boost/math/special_functions/math_fwd.hpp>
15 #include <boost/math/special_functions/factorials.hpp>
16 #include <boost/math/tools/config.hpp>
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)
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);
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)
36 static const char* function = "boost::math::legrendre_p<%1%>(unsigned, %1%)";
38 if((x < -1) || (x > 1))
39 return policies::raise_domain_error<T>(
41 "The Legendre Polynomial is defined for"
42 " -1 <= x <= 1, but got x = %1%.", x, pol);
47 // A solution of the second kind (Q):
48 p0 = (boost::math::log1p(x, pol) - boost::math::log1p(-x, pol)) / 2;
53 // A solution of the first kind (P):
65 p1 = boost::math::legendre_next(n, x, p0, p1);
73 template <class T, class Policy>
74 inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
75 legendre_p(int l, T x, const Policy& pol)
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%)";
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);
86 inline typename tools::promote_args<T>::type
87 legendre_p(int l, T x)
89 return boost::math::legendre_p(l, x, policies::policy<>());
92 template <class T, class Policy>
93 inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type
94 legendre_q(unsigned l, T x, const Policy& pol)
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%)");
102 inline typename tools::promote_args<T>::type
103 legendre_q(unsigned l, T x)
105 return boost::math::legendre_q(l, x, policies::policy<>());
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)
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);
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)
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:
130 return legendre_p_imp(-l-1, m, x, sin_theta_power, pol);
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);
140 return boost::math::legendre_p(l, x, pol);
142 T p0 = boost::math::double_factorial<T>(2 * m - 1, pol) * sin_theta_power;
149 T p1 = x * (2 * m + 1) * p0;
156 p1 = boost::math::legendre_next(n, m, x, p0, p1);
162 template <class T, class Policy>
163 inline T legendre_p_imp(int l, int m, T x, const Policy& pol)
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);
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)
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%)");
182 inline typename tools::promote_args<T>::type
183 legendre_p(int l, int m, T x)
185 return boost::math::legendre_p(l, m, x, policies::policy<>());
191 #endif // BOOST_MATH_SPECIAL_LEGENDRE_HPP