1 // (C) Copyright Nick Thompson 2019.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 #ifndef BOOST_MATH_SPECIAL_GEGENBAUER_HPP
7 #define BOOST_MATH_SPECIAL_GEGENBAUER_HPP
8
9 #include <stdexcept>
10 #include <type_traits>
11
12 namespace boost { namespace math {
13
14 template<typename Real>
gegenbauer(unsigned n,Real lambda,Real x)15 Real gegenbauer(unsigned n, Real lambda, Real x)
16 {
17 static_assert(!std::is_integral<Real>::value, "Gegenbauer polynomials required floating point arguments.");
18 if (lambda <= -1/Real(2)) {
19 throw std::domain_error("lambda > -1/2 is required.");
20 }
21 // The only reason to do this is because of some instability that could be present for x < 0 that is not present for x > 0.
22 // I haven't observed this, but then again, I haven't managed to test an exhaustive number of parameters.
23 // In any case, the routine is distinctly faster without this test:
24 //if (x < 0) {
25 // if (n&1) {
26 // return -gegenbauer(n, lambda, -x);
27 // }
28 // return gegenbauer(n, lambda, -x);
29 //}
30
31 if (n == 0) {
32 return Real(1);
33 }
34 Real y0 = 1;
35 Real y1 = 2*lambda*x;
36
37 Real yk = y1;
38 Real k = 2;
39 Real k_max = n*(1+std::numeric_limits<Real>::epsilon());
40 Real gamma = 2*(lambda - 1);
41 while(k < k_max)
42 {
43 yk = ( (2 + gamma/k)*x*y1 - (1+gamma/k)*y0);
44 y0 = y1;
45 y1 = yk;
46 k += 1;
47 }
48 return yk;
49 }
50
51
52 template<typename Real>
gegenbauer_derivative(unsigned n,Real lambda,Real x,unsigned k)53 Real gegenbauer_derivative(unsigned n, Real lambda, Real x, unsigned k)
54 {
55 if (k > n) {
56 return Real(0);
57 }
58 Real gegen = gegenbauer<Real>(n-k, lambda + k, x);
59 Real scale = 1;
60 for (unsigned j = 0; j < k; ++j) {
61 scale *= 2*lambda;
62 lambda += 1;
63 }
64 return scale*gegen;
65 }
66
67 template<typename Real>
gegenbauer_prime(unsigned n,Real lambda,Real x)68 Real gegenbauer_prime(unsigned n, Real lambda, Real x) {
69 return gegenbauer_derivative<Real>(n, lambda, x, 1);
70 }
71
72
73 }}
74 #endif
75