math/special_functions/jacobi.hpp

0001 //  (C) Copyright Nick Thompson 2019.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SPECIAL_JACOBI_HPP
0007 #define BOOST_MATH_SPECIAL_JACOBI_HPP
0008
0009 #include <limits>
0010 #include <stdexcept>
0011
0012 namespace boost { namespace math {
0013
0014 template<typename Real>
0015 Real jacobi(unsigned n, Real alpha, Real beta, Real x)
0016 {
0017     static_assert(!std::is_integral<Real>::value, "Jacobi polynomials do not work with integer arguments.");
0018
0019     if (n == 0) {
0020         return Real(1);
0021     }
0022     Real y0 = 1;
0023     Real y1 = (alpha+1) + (alpha+beta+2)*(x-1)/Real(2);
0024
0025     Real yk = y1;
0026     Real k = 2;
0027     Real k_max = n*(1+std::numeric_limits<Real>::epsilon());
0028     while(k < k_max)
0029     {
0030         // Hoping for lots of common subexpression elimination by the compiler:
0031         Real denom = 2*k*(k+alpha+beta)*(2*k+alpha+beta-2);
0032         Real gamma1 = (2*k+alpha+beta-1)*( (2*k+alpha+beta)*(2*k+alpha+beta-2)*x + alpha*alpha -beta*beta);
0033         Real gamma0 = -2*(k+alpha-1)*(k+beta-1)*(2*k+alpha+beta);
0034         yk = (gamma1*y1 + gamma0*y0)/denom;
0035         y0 = y1;
0036         y1 = yk;
0037         k += 1;
0038     }
0039     return yk;
0040 }
0041
0042 template<typename Real>
0043 Real jacobi_derivative(unsigned n, Real alpha, Real beta, Real x, unsigned k)
0044 {
0045     if (k > n) {
0046         return Real(0);
0047     }
0048     Real scale = 1;
0049     for(unsigned j = 1; j <= k; ++j) {
0050         scale *= (alpha + beta + n + j)/2;
0051     }
0052
0053     return scale*jacobi<Real>(n-k, alpha + k, beta+k, x);
0054 }
0055
0056 template<typename Real>
0057 Real jacobi_prime(unsigned n, Real alpha, Real beta, Real x)
0058 {
0059     return jacobi_derivative<Real>(n, alpha, beta, x, 1);
0060 }
0061
0062 template<typename Real>
0063 Real jacobi_double_prime(unsigned n, Real alpha, Real beta, Real x)
0064 {
0065     return jacobi_derivative<Real>(n, alpha, beta, x, 2);
0066 }
0067
0068 }}
0069 #endif