diff options
Diffstat (limited to 'boost/math/special_functions/legendre.hpp')
-rw-r--r-- | boost/math/special_functions/legendre.hpp | 191 |
1 files changed, 185 insertions, 6 deletions
diff --git a/boost/math/special_functions/legendre.hpp b/boost/math/special_functions/legendre.hpp index 1a2ef5d615..6028b377d3 100644 --- a/boost/math/special_functions/legendre.hpp +++ b/boost/math/special_functions/legendre.hpp @@ -11,8 +11,11 @@ #pragma once #endif +#include <utility> +#include <vector> #include <boost/math/special_functions/math_fwd.hpp> #include <boost/math/special_functions/factorials.hpp> +#include <boost/math/tools/roots.hpp> #include <boost/math/tools/config.hpp> namespace boost{ @@ -68,6 +71,149 @@ T legendre_imp(unsigned l, T x, const Policy& pol, bool second = false) return p1; } +template <class T, class Policy> +T legendre_p_prime_imp(unsigned l, T x, const Policy& pol, T* Pn +#ifdef BOOST_NO_CXX11_NULLPTR + = 0 +#else + = nullptr +#endif +) +{ + static const char* function = "boost::math::legrendre_p_prime<%1%>(unsigned, %1%)"; + // Error handling: + if ((x < -1) || (x > 1)) + return policies::raise_domain_error<T>( + function, + "The Legendre Polynomial is defined for" + " -1 <= x <= 1, but got x = %1%.", x, pol); + + if (l == 0) + { + if (Pn) + { + *Pn = 1; + } + return 0; + } + T p0 = 1; + T p1 = x; + T p_prime; + bool odd = l & 1; + // If the order is odd, we sum all the even polynomials: + if (odd) + { + p_prime = p0; + } + else // Otherwise we sum the odd polynomials * (2n+1) + { + p_prime = 3*p1; + } + + unsigned n = 1; + while(n < l - 1) + { + std::swap(p0, p1); + p1 = boost::math::legendre_next(n, x, p0, p1); + ++n; + if (odd) + { + p_prime += (2*n+1)*p1; + odd = false; + } + else + { + odd = true; + } + } + // This allows us to evaluate the derivative and the function for the same cost. + if (Pn) + { + std::swap(p0, p1); + *Pn = boost::math::legendre_next(n, x, p0, p1); + } + return p_prime; +} + +template <class T, class Policy> +struct legendre_p_zero_func +{ + int n; + const Policy& pol; + + legendre_p_zero_func(int n_, const Policy& p) : n(n_), pol(p) {} + + std::pair<T, T> operator()(T x) const + { + T Pn; + T Pn_prime = detail::legendre_p_prime_imp(n, x, pol, &Pn); + return std::pair<T, T>(Pn, Pn_prime); + }; +}; + +template <class T, class Policy> +std::vector<T> legendre_p_zeros_imp(int n, const Policy& pol) +{ + using std::cos; + using std::sin; + using std::ceil; + using std::sqrt; + using boost::math::constants::pi; + using boost::math::constants::half; + using boost::math::tools::newton_raphson_iterate; + + BOOST_ASSERT(n >= 0); + std::vector<T> zeros; + if (n == 0) + { + // There are no zeros of P_0(x) = 1. + return zeros; + } + int k; + if (n & 1) + { + zeros.resize((n-1)/2 + 1, std::numeric_limits<T>::quiet_NaN()); + zeros[0] = 0; + k = 1; + } + else + { + zeros.resize(n/2, std::numeric_limits<T>::quiet_NaN()); + k = 0; + } + T half_n = ceil(n*half<T>()); + + while (k < (int)zeros.size()) + { + // Bracket the root: Szego: + // Gabriel Szego, Inequalities for the Zeros of Legendre Polynomials and Related Functions, Transactions of the American Mathematical Society, Vol. 39, No. 1 (1936) + T theta_nk = ((half_n - half<T>()*half<T>() - static_cast<T>(k))*pi<T>())/(static_cast<T>(n)+half<T>()); + T lower_bound = cos( (half_n - static_cast<T>(k))*pi<T>()/static_cast<T>(n + 1)); + T cos_nk = cos(theta_nk); + T upper_bound = cos_nk; + // First guess follows from: + // F. G. Tricomi, Sugli zeri dei polinomi sferici ed ultrasferici, Ann. Mat. Pura Appl., 31 (1950), pp. 93–97; + T inv_n_sq = 1/static_cast<T>(n*n); + T sin_nk = sin(theta_nk); + 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; + + boost::uintmax_t number_of_iterations = policies::get_max_root_iterations<Policy>(); + + legendre_p_zero_func<T, Policy> f(n, pol); + + const T x_nk = newton_raphson_iterate(f, x_nk_guess, + lower_bound, upper_bound, + policies::digits<T, Policy>(), + number_of_iterations); + + BOOST_ASSERT(lower_bound < x_nk); + BOOST_ASSERT(upper_bound > x_nk); + zeros[k] = x_nk; + ++k; + } + return zeros; +} + } // namespace detail template <class T, class Policy> @@ -82,13 +228,49 @@ inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_imp(l, static_cast<value_type>(x), pol, false), function); } + +template <class T, class Policy> +inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type + legendre_p_prime(int l, T x, const Policy& pol) +{ + typedef typename tools::promote_args<T>::type result_type; + typedef typename policies::evaluation<result_type, Policy>::type value_type; + static const char* function = "boost::math::legendre_p_prime<%1%>(unsigned, %1%)"; + if(l < 0) + return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(-l-1, static_cast<value_type>(x), pol), function); + return policies::checked_narrowing_cast<result_type, Policy>(detail::legendre_p_prime_imp(l, static_cast<value_type>(x), pol), function); +} + template <class T> -inline typename tools::promote_args<T>::type +inline typename tools::promote_args<T>::type legendre_p(int l, T x) { return boost::math::legendre_p(l, x, policies::policy<>()); } +template <class T> +inline typename tools::promote_args<T>::type + legendre_p_prime(int l, T x) +{ + return boost::math::legendre_p_prime(l, x, policies::policy<>()); +} + +template <class T, class Policy> +inline std::vector<T> legendre_p_zeros(int l, const Policy& pol) +{ + if(l < 0) + return detail::legendre_p_zeros_imp<T>(-l-1, pol); + + return detail::legendre_p_zeros_imp<T>(l, pol); +} + + +template <class T> +inline std::vector<T> legendre_p_zeros(int l) +{ + return boost::math::legendre_p_zeros<T>(l, policies::policy<>()); +} + template <class T, class Policy> inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename tools::promote_args<T>::type>::type legendre_q(unsigned l, T x, const Policy& pol) @@ -99,7 +281,7 @@ inline typename boost::enable_if_c<policies::is_policy<Policy>::value, typename } template <class T> -inline typename tools::promote_args<T>::type +inline typename tools::promote_args<T>::type legendre_q(unsigned l, T x) { return boost::math::legendre_q(l, x, policies::policy<>()); @@ -107,7 +289,7 @@ inline typename tools::promote_args<T>::type // Recurrence for associated polynomials: template <class T1, class T2, class T3> -inline typename tools::promote_args<T1, T2, T3>::type +inline typename tools::promote_args<T1, T2, T3>::type legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1) { typedef typename tools::promote_args<T1, T2, T3>::type result_type; @@ -189,6 +371,3 @@ inline typename tools::promote_args<T>::type } // namespace boost #endif // BOOST_MATH_SPECIAL_LEGENDRE_HPP - - - |