]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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 | ||
1e59de90 | 9 | #include <limits> |
92f5a8d4 TL |
10 | #include <stdexcept> |
11 | #include <type_traits> | |
12 | ||
13 | namespace boost { namespace math { | |
14 | ||
15 | template<typename Real> | |
16 | Real gegenbauer(unsigned n, Real lambda, Real x) | |
17 | { | |
18 | static_assert(!std::is_integral<Real>::value, "Gegenbauer polynomials required floating point arguments."); | |
19 | if (lambda <= -1/Real(2)) { | |
20 | throw std::domain_error("lambda > -1/2 is required."); | |
21 | } | |
22 | // 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. | |
23 | // I haven't observed this, but then again, I haven't managed to test an exhaustive number of parameters. | |
24 | // In any case, the routine is distinctly faster without this test: | |
25 | //if (x < 0) { | |
26 | // if (n&1) { | |
27 | // return -gegenbauer(n, lambda, -x); | |
28 | // } | |
29 | // return gegenbauer(n, lambda, -x); | |
30 | //} | |
31 | ||
32 | if (n == 0) { | |
33 | return Real(1); | |
34 | } | |
35 | Real y0 = 1; | |
36 | Real y1 = 2*lambda*x; | |
37 | ||
38 | Real yk = y1; | |
39 | Real k = 2; | |
40 | Real k_max = n*(1+std::numeric_limits<Real>::epsilon()); | |
41 | Real gamma = 2*(lambda - 1); | |
42 | while(k < k_max) | |
43 | { | |
44 | yk = ( (2 + gamma/k)*x*y1 - (1+gamma/k)*y0); | |
45 | y0 = y1; | |
46 | y1 = yk; | |
47 | k += 1; | |
48 | } | |
49 | return yk; | |
50 | } | |
51 | ||
52 | ||
53 | template<typename Real> | |
54 | Real gegenbauer_derivative(unsigned n, Real lambda, Real x, unsigned k) | |
55 | { | |
56 | if (k > n) { | |
57 | return Real(0); | |
58 | } | |
59 | Real gegen = gegenbauer<Real>(n-k, lambda + k, x); | |
60 | Real scale = 1; | |
61 | for (unsigned j = 0; j < k; ++j) { | |
62 | scale *= 2*lambda; | |
63 | lambda += 1; | |
64 | } | |
65 | return scale*gegen; | |
66 | } | |
67 | ||
68 | template<typename Real> | |
69 | Real gegenbauer_prime(unsigned n, Real lambda, Real x) { | |
70 | return gegenbauer_derivative<Real>(n, lambda, x, 1); | |
71 | } | |
72 | ||
73 | ||
74 | }} | |
75 | #endif |