1[section:legendre Legendre (and Associated) Polynomials] 2 3[h4 Synopsis] 4 5`` 6#include <boost/math/special_functions/legendre.hpp> 7`` 8 9 namespace boost{ namespace math{ 10 11 template <class T> 12 ``__sf_result`` legendre_p(int n, T x); 13 14 template <class T, class ``__Policy``> 15 ``__sf_result`` legendre_p(int n, T x, const ``__Policy``&); 16 17 template <class T> 18 ``__sf_result`` legendre_p_prime(int n, T x); 19 20 template <class T, class ``__Policy``> 21 ``__sf_result`` legendre_p_prime(int n, T x, const ``__Policy``&); 22 23 template <class T, class ``__Policy``> 24 std::vector<T> legendre_p_zeros(int l, const ``__Policy``&); 25 26 template <class T> 27 std::vector<T> legendre_p_zeros(int l); 28 29 template <class T> 30 ``__sf_result`` legendre_p(int n, int m, T x); 31 32 template <class T, class ``__Policy``> 33 ``__sf_result`` legendre_p(int n, int m, T x, const ``__Policy``&); 34 35 template <class T> 36 ``__sf_result`` legendre_q(unsigned n, T x); 37 38 template <class T, class ``__Policy``> 39 ``__sf_result`` legendre_q(unsigned n, T x, const ``__Policy``&); 40 41 template <class T1, class T2, class T3> 42 ``__sf_result`` legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1); 43 44 template <class T1, class T2, class T3> 45 ``__sf_result`` legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1); 46 47 48 }} // namespaces 49 50The return type of these functions is computed using the __arg_promotion_rules: 51note than when there is a single template argument the result is the same type 52as that argument or `double` if the template argument is an integer type. 53 54[optional_policy] 55 56[h4 Description] 57 58 template <class T> 59 ``__sf_result`` legendre_p(int l, T x); 60 61 template <class T, class ``__Policy``> 62 ``__sf_result`` legendre_p(int l, T x, const ``__Policy``&); 63 64Returns the Legendre Polynomial of the first kind: 65 66[equation legendre_0] 67 68Requires -1 <= x <= 1, otherwise returns the result of __domain_error. 69 70Negative orders are handled via the reflection formula: 71 72[:P[sub -l-1](x) = P[sub l](x)] 73 74The following graph illustrates the behaviour of the first few 75Legendre Polynomials: 76 77[graph legendre_p] 78 79 template <class T> 80 ``__sf_result`` legendre_p_prime(int n, T x); 81 82 template <class T, class ``__Policy``> 83 ``__sf_result`` legendre_p_prime(int n, T x, const ``__Policy``&); 84 85Returns the derivatives of the Legendre polynomials. 86 87 template <class T, class ``__Policy``> 88 std::vector<T> legendre_p_zeros(int l, const ``__Policy``&); 89 90 template <class T> 91 std::vector<T> legendre_p_zeros(int l); 92 93The zeros of the Legendre polynomials are calculated by Newton's method using an initial guess given by Tricomi with root bracketing provided by Szego. 94 95Since the Legendre polynomials are alternatively even and odd, only the non-negative zeros are returned. 96For the odd Legendre polynomials, the first zero is always zero. 97The rest of the zeros are returned in increasing order. 98 99Note that the argument to the routine is an integer, and the output is a floating-point type. 100Hence the template argument is mandatory. 101The time to extract a single root is linear in `l` (this is scaling to evaluate the Legendre polynomials), so recovering all roots is [bigo](`l`[super 2]). 102Algorithms with linear scaling [@ https://doi.org/10.1137/06067016X exist] for recovering all roots, but requires tooling not currently built into boost.math. 103This implementation proceeds under the assumption that calculating zeros of these functions will not be a bottleneck for any workflow. 104 105 template <class T> 106 ``__sf_result`` legendre_p(int l, int m, T x); 107 108 template <class T, class ``__Policy``> 109 ``__sf_result`` legendre_p(int l, int m, T x, const ``__Policy``&); 110 111Returns the associated Legendre polynomial of the first kind: 112 113[equation legendre_1] 114 115Requires -1 <= x <= 1, otherwise returns the result of __domain_error. 116 117Negative values of /l/ and /m/ are handled via the identity relations: 118 119[equation legendre_3] 120 121[caution The definition of the associated Legendre polynomial used here 122includes a leading Condon-Shortley phase term of (-1)[super m]. This 123matches the definition given by Abramowitz and Stegun (8.6.6) and that 124used by [@http://mathworld.wolfram.com/LegendrePolynomial.html Mathworld] 125and [@http://documents.wolfram.com/mathematica/functions/LegendreP 126Mathematica's LegendreP function]. However, uses in the literature 127do not always include this phase term, and strangely the specification 128for the associated Legendre function in the C++ TR1 (assoc_legendre) 129also omits it, in spite of stating that it uses Abramowitz and Stegun 130as the final arbiter on these matters. 131 132See: 133 134[@http://mathworld.wolfram.com/LegendrePolynomial.html 135Weisstein, Eric W. "Legendre Polynomial." 136From MathWorld--A Wolfram Web Resource]. 137 138Abramowitz, M. and Stegun, I. A. (Eds.). "Legendre Functions" and 139"Orthogonal Polynomials." Ch. 22 in Chs. 8 and 22 in Handbook of 140Mathematical Functions with Formulas, Graphs, and Mathematical Tables, 1419th printing. New York: Dover, pp. 331-339 and 771-802, 1972. 142 ] 143 144 template <class T> 145 ``__sf_result`` legendre_q(unsigned n, T x); 146 147 template <class T, class ``__Policy``> 148 ``__sf_result`` legendre_q(unsigned n, T x, const ``__Policy``&); 149 150Returns the value of the Legendre polynomial that is the second solution 151to the Legendre differential equation, for example: 152 153[equation legendre_2] 154 155Requires -1 <= x <= 1, otherwise __domain_error is called. 156 157The following graph illustrates the first few Legendre functions of the 158second kind: 159 160[graph legendre_q] 161 162 template <class T1, class T2, class T3> 163 ``__sf_result`` legendre_next(unsigned l, T1 x, T2 Pl, T3 Plm1); 164 165Implements the three term recurrence relation for the Legendre 166polynomials, this function can be used to create a sequence of 167values evaluated at the same /x/, and for rising /l/. This recurrence 168relation holds for Legendre Polynomials of both the first and second kinds. 169 170[equation legendre_4] 171 172For example we could produce a vector of the first 10 polynomial 173values using: 174 175 double x = 0.5; // Abscissa value 176 vector<double> v; 177 v.push_back(legendre_p(0, x)); 178 v.push_back(legendre_p(1, x)); 179 for(unsigned l = 1; l < 10; ++l) 180 v.push_back(legendre_next(l, x, v[l], v[l-1])); 181 // Double check values: 182 for(unsigned l = 1; l < 10; ++l) 183 assert(v[l] == legendre_p(l, x)); 184 185Formally the arguments are: 186 187[variablelist 188[[l][The degree of the last polynomial calculated.]] 189[[x][The abscissa value]] 190[[Pl][The value of the polynomial evaluated at degree /l/.]] 191[[Plm1][The value of the polynomial evaluated at degree /l-1/.]] 192] 193 194 template <class T1, class T2, class T3> 195 ``__sf_result`` legendre_next(unsigned l, unsigned m, T1 x, T2 Pl, T3 Plm1); 196 197Implements the three term recurrence relation for the Associated Legendre 198polynomials, this function can be used to create a sequence of 199values evaluated at the same /x/, and for rising /l/. 200 201[equation legendre_5] 202 203For example we could produce a vector of the first m+10 polynomial 204values using: 205 206 double x = 0.5; // Abscissa value 207 int m = 10; // order 208 vector<double> v; 209 v.push_back(legendre_p(m, m, x)); 210 v.push_back(legendre_p(1 + m, m, x)); 211 for(unsigned l = 1; l < 10; ++l) 212 v.push_back(legendre_next(l + 10, m, x, v[l], v[l-1])); 213 // Double check values: 214 for(unsigned l = 1; l < 10; ++l) 215 assert(v[l] == legendre_p(10 + l, m, x)); 216 217Formally the arguments are: 218 219[variablelist 220[[l][The degree of the last polynomial calculated.]] 221[[m][The order of the Associated Polynomial.]] 222[[x][The abscissa value]] 223[[Pl][The value of the polynomial evaluated at degree /l/.]] 224[[Plm1][The value of the polynomial evaluated at degree /l-1/.]] 225] 226 227[h4 Accuracy] 228 229The following table shows peak errors (in units of epsilon) 230for various domains of input arguments. 231Note that only results for the widest floating point type on the system are 232given as narrower types have __zero_error. 233 234[table_legendre_p] 235 236[table_legendre_q] 237 238[table_legendre_p_associated_] 239 240Note that the worst errors occur when the order increases, values greater than 241~120 are very unlikely to produce sensible results, especially in the associated 242polynomial case when the degree is also large. Further the relative errors 243are likely to grow arbitrarily large when the function is very close to a root. 244 245[h4 Testing] 246 247A mixture of spot tests of values calculated using functions.wolfram.com, 248and randomly generated test data are 249used: the test data was computed using 250[@http://shoup.net/ntl/doc/RR.txt NTL::RR] at 1000-bit precision. 251 252[h4 Implementation] 253 254These functions are implemented using the stable three term 255recurrence relations. These relations guarantee low absolute error 256but cannot guarantee low relative error near one of the roots of the 257polynomials. 258 259[endsect] [/section:beta_function The Beta Function] 260[/ 261 Copyright 2006 John Maddock and Paul A. Bristow. 262 Distributed under the Boost Software License, Version 1.0. 263 (See accompanying file LICENSE_1_0.txt or copy at 264 http://www.boost.org/LICENSE_1_0.txt). 265] 266