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_JACOBI_HPP
7 #define BOOST_MATH_SPECIAL_JACOBI_HPP
8
9 #include <limits>
10 #include <stdexcept>
11
12 namespace boost { namespace math {
13
14 template<typename Real>
jacobi(unsigned n,Real alpha,Real beta,Real x)15 Real jacobi(unsigned n, Real alpha, Real beta, Real x)
16 {
17 static_assert(!std::is_integral<Real>::value, "Jacobi polynomials do not work with integer arguments.");
18
19 if (n == 0) {
20 return Real(1);
21 }
22 Real y0 = 1;
23 Real y1 = (alpha+1) + (alpha+beta+2)*(x-1)/Real(2);
24
25 Real yk = y1;
26 Real k = 2;
27 Real k_max = n*(1+std::numeric_limits<Real>::epsilon());
28 while(k < k_max)
29 {
30 // Hoping for lots of common subexpression elimination by the compiler:
31 Real denom = 2*k*(k+alpha+beta)*(2*k+alpha+beta-2);
32 Real gamma1 = (2*k+alpha+beta-1)*( (2*k+alpha+beta)*(2*k+alpha+beta-2)*x + alpha*alpha -beta*beta);
33 Real gamma0 = -2*(k+alpha-1)*(k+beta-1)*(2*k+alpha+beta);
34 yk = (gamma1*y1 + gamma0*y0)/denom;
35 y0 = y1;
36 y1 = yk;
37 k += 1;
38 }
39 return yk;
40 }
41
42 template<typename Real>
jacobi_derivative(unsigned n,Real alpha,Real beta,Real x,unsigned k)43 Real jacobi_derivative(unsigned n, Real alpha, Real beta, Real x, unsigned k)
44 {
45 if (k > n) {
46 return Real(0);
47 }
48 Real scale = 1;
49 for(unsigned j = 1; j <= k; ++j) {
50 scale *= (alpha + beta + n + j)/2;
51 }
52
53 return scale*jacobi<Real>(n-k, alpha + k, beta+k, x);
54 }
55
56 template<typename Real>
jacobi_prime(unsigned n,Real alpha,Real beta,Real x)57 Real jacobi_prime(unsigned n, Real alpha, Real beta, Real x)
58 {
59 return jacobi_derivative<Real>(n, alpha, beta, x, 1);
60 }
61
62 template<typename Real>
jacobi_double_prime(unsigned n,Real alpha,Real beta,Real x)63 Real jacobi_double_prime(unsigned n, Real alpha, Real beta, Real x)
64 {
65 return jacobi_derivative<Real>(n, alpha, beta, x, 2);
66 }
67
68 }}
69 #endif
70