1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374 |
- #ifndef BOOST_MATH_SPECIAL_GEGENBAUER_HPP
- #define BOOST_MATH_SPECIAL_GEGENBAUER_HPP
- #include <stdexcept>
- #include <type_traits>
- namespace boost { namespace math {
- template<typename Real>
- Real gegenbauer(unsigned n, Real lambda, Real x)
- {
- static_assert(!std::is_integral<Real>::value, "Gegenbauer polynomials required floating point arguments.");
- if (lambda <= -1/Real(2)) {
- throw std::domain_error("lambda > -1/2 is required.");
- }
-
-
-
-
-
-
-
-
-
- if (n == 0) {
- return Real(1);
- }
- Real y0 = 1;
- Real y1 = 2*lambda*x;
- Real yk = y1;
- Real k = 2;
- Real k_max = n*(1+std::numeric_limits<Real>::epsilon());
- Real gamma = 2*(lambda - 1);
- while(k < k_max)
- {
- yk = ( (2 + gamma/k)*x*y1 - (1+gamma/k)*y0);
- y0 = y1;
- y1 = yk;
- k += 1;
- }
- return yk;
- }
- template<typename Real>
- Real gegenbauer_derivative(unsigned n, Real lambda, Real x, unsigned k)
- {
- if (k > n) {
- return Real(0);
- }
- Real gegen = gegenbauer<Real>(n-k, lambda + k, x);
- Real scale = 1;
- for (unsigned j = 0; j < k; ++j) {
- scale *= 2*lambda;
- lambda += 1;
- }
- return scale*gegen;
- }
- template<typename Real>
- Real gegenbauer_prime(unsigned n, Real lambda, Real x) {
- return gegenbauer_derivative<Real>(n, lambda, x, 1);
- }
- }}
- #endif
|