70 lines · plain
1// (C) Copyright Nick Thompson 2019.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_SPECIAL_JACOBI_HPP7#define BOOST_MATH_SPECIAL_JACOBI_HPP8 9#include <limits>10#include <stdexcept>11 12namespace boost { namespace math {13 14template<typename Real>15Real 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 42template<typename Real>43Real 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 56template<typename Real>57Real jacobi_prime(unsigned n, Real alpha, Real beta, Real x)58{59 return jacobi_derivative<Real>(n, alpha, beta, x, 1);60}61 62template<typename Real>63Real 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#endif70