polynomial.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2000 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_polynomial_h
#define dealii_polynomial_h
 
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/numbers.h>
#include <deal.II/base/point.h>
#include <deal.II/base/types.h>
 
#include <array>
#include <limits>
#include <memory>
#include <shared_mutex>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN


namespace Polynomials
{
  template <typename number>
  class Polynomial : public EnableObserverPointer
  {
  public:
    Polynomial(const std::vector<number> &coefficients);

    Polynomial(const unsigned int n);

    Polynomial(const std::vector<Point<1>> &lagrange_support_points,
               const unsigned int           j);

    Polynomial();

    number
    value(const number x) const;

    void
    value(const number x, std::vector<number> &values) const;

    template <typename Number2>
    void
    value(const Number2      x,
          const unsigned int n_derivatives,
          Number2           *values) const;

    template <std::size_t n_entries, typename Number2>
#ifndef DEBUG
    DEAL_II_ALWAYS_INLINE
#endif
      void
      values_of_array(const std::array<Number2, n_entries> &points,
                      const unsigned int                    n_derivatives,
                      std::array<Number2, n_entries>       *values) const;

    unsigned int
    degree() const;

    void
    scale(const number factor);

    template <typename number2>
    void
    shift(const number2 offset);

    Polynomial<number>
    derivative() const;

    Polynomial<number>
    primitive() const;

    Polynomial<number> &
    operator*=(const double s);

    Polynomial<number> &
    operator*=(const Polynomial<number> &p);

    Polynomial<number> &
    operator+=(const Polynomial<number> &p);

    Polynomial<number> &
    operator-=(const Polynomial<number> &p);

    bool
    operator==(const Polynomial<number> &p) const;

    void
    print(std::ostream &out) const;

    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int version);

    virtual std::size_t
    memory_consumption() const;
 
  protected:
    static void
    scale(std::vector<number> &coefficients, const number factor);

    template <typename number2>
    static void
    shift(std::vector<number> &coefficients, const number2 shift);

    static void
    multiply(std::vector<number> &coefficients, const number factor);

    void
    transform_into_standard_form();

    std::vector<number> coefficients;

    bool in_lagrange_product_form;

    std::vector<number> lagrange_support_points;

    number lagrange_weight;
  };
 

  template <typename number>
  class Monomial : public Polynomial<number>
  {
  public:
    Monomial(const unsigned int n, const double coefficient = 1.);

    static std::vector<Polynomial<number>>
    generate_complete_basis(const unsigned int degree);
 
  private:
    static std::vector<Point<1>>
    create_vector_of_roots(unsigned int n);
  };
 

  class LagrangeEquidistant : public Polynomial<double>
  {
  public:
    LagrangeEquidistant(const unsigned int n, const unsigned int support_point);

    static std::vector<Polynomial<double>>
    generate_complete_basis(const unsigned int degree);
 
  private:
    static void
    compute_coefficients(const unsigned int   n,
                         const unsigned int   support_point,
                         std::vector<double> &a);
  };
 
 

  std::vector<Polynomial<double>>
  generate_complete_Lagrange_basis(const std::vector<Point<1>> &points);
 
 

  class Legendre : public Polynomial<double>
  {
  public:
    Legendre(const unsigned int p);

    static std::vector<Polynomial<double>>
    generate_complete_basis(const unsigned int degree);
  };

  class Lobatto : public Polynomial<double>
  {
  public:
    Lobatto(const unsigned int p = 0);

    static std::vector<Polynomial<double>>
    generate_complete_basis(const unsigned int p);
 
  private:
    std::vector<double>
    compute_coefficients(const unsigned int p);
  };
 
 

  class Hierarchical : public Polynomial<double>
  {
  public:
    Hierarchical(const unsigned int p);

    static std::vector<Polynomial<double>>
    generate_complete_basis(const unsigned int degree);
 
  private:
    static void
    compute_coefficients(const unsigned int p);

    static const std::vector<double> &
    get_coefficients(const unsigned int p);

    static std::vector<std::unique_ptr<const std::vector<double>>>
      recursive_coefficients;

    static std::shared_mutex coefficients_lock;
  };
 
 

  class HermiteInterpolation : public Polynomial<double>
  {
  public:
    HermiteInterpolation(const unsigned int p);

    static std::vector<Polynomial<double>>
    generate_complete_basis(const unsigned int p);
  };
 
 

  class HermiteLikeInterpolation : public Polynomial<double>
  {
  public:
    HermiteLikeInterpolation(const unsigned int degree,
                             const unsigned int index);

    static std::vector<Polynomial<double>>
    generate_complete_basis(const unsigned int degree);
  };
 
 
 
  /*
   * Evaluate a Jacobi polynomial @f$ P_n^{\alpha, \beta}(x) @f$ specified by the
   * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
   * Jacobi polynomial.
   *
   * @note The Jacobi polynomials are not orthonormal and are defined on the
   * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
   * used in literature. @p x is the point of evaluation. If instead the point
   * @p x is given on the interval @f$[-1, +1]@f$ set
   * @p rescale_to_dealii_unit_interval to false.
   */
  template <typename Number>
  Number
  jacobi_polynomial_value(const unsigned int degree,
                          const int          alpha,
                          const int          beta,
                          const Number       x,
                          const bool rescale_to_dealii_unit_interval = true);
 
  /*
   * Evaluate the derivative of the Jacobi polynomial @f$ d P_n^{\alpha, \beta}(x)
   * / dx @f$ specified by the
   * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
   * Jacobi polynomial.
   *
   * @note The Jacobi polynomials are not orthonormal and are defined on the
   * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
   * used in literature. @p x is the point of evaluation. If instead the point @p
   * x is given on the interval @f$[-1, +1]@f$ set @p rescale_to_dealii_unit_interval
   * to false.
   */
  template <typename Number>
  Number
  jacobi_polynomial_derivative(const unsigned int degree,
                               const int          alpha,
                               const int          beta,
                               const Number       x,
                               const bool rescale_to_dealii_unit_interval);
 
  /*
   * Evaluate k-th derivative of the Jacobi polynomial
   * @f$ d^k P_n^{\alpha, \beta}(x) / dx^k @f$, where k equals @p derivative_order,
   * specified by the parameters @p alpha, @p beta, @p n, where @p n is the
   * degree of the Jacobi polynomial and @k the degree of the derivative.
   *
   * @note The Jacobi polynomials are not orthonormal and are defined on the
   * unit interval @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often
   * used in literature. @p x is the point of evaluation. If instead the point @p
   * x is given on the interval @f$[-1, +1]@f$ set @p rescale_to_dealii_unit_interval
   * to false.
   */
  template <typename Number>
  Number
  jacobi_polynomial_kth_derivative(const unsigned int derivative_order,
                                   const unsigned int degree,
                                   const int          alpha,
                                   const int          beta,
                                   const Number       x,
                                   const bool rescale_to_dealii_unit_interval);

  template <typename Number>
  std::vector<Number>
  jacobi_polynomial_roots(const unsigned int degree,
                          const int          alpha,
                          const int          beta);
 
  /*
   * Evaluate the homogenized Jacobi polynomial
   * @f$ Q_n^{\alpha, \beta}(x,s) = s^n P_n^{\alpha, \beta}(x/s) @f$
   * specified by the parameters @p alpha, @p beta, @p n, where @p n is the
   * degree of the Jacobi polynomial. The resulting polynomial has no
   * singularity, the computation of @f$x/s@f$ is explicitly avoided.
   *
   * @note The homogenized Jacobi polynomials are defined on the unit interval
   * @f$[0, 1]@f$ as usual for deal.II, rather than @f$[-1, +1]@f$ often used in
   * literature. @p x is the point of evaluation, the shift to the interval
   * @f$[-1, +1]@f$ is conducted as
   * @f$ Q_n^{\alpha, \beta}(x,s) = s^n P_n^{\alpha, \beta}(2 x/s - 1) @f$.
   */
  template <typename Number>
  Number
  jacobi_polynomial_homogenized_value(const unsigned int degree,
                                      const int          alpha,
                                      const int          beta,
                                      const Number       x,
                                      const Number       s);
 
  /*
   * Evaluate (mixed) derivatives of the homogenized Jacobi polynomial
   * @f$ d^k Q_n^{\alpha, \beta}(x,s) / dx^order_x ds^order_s @f$ specified by the
   * parameters @p alpha, @p beta, @p n, where @p n is the degree of the
   * Jacobi polynomial and @order_x and @order_s the degree of the derivative
   * with respect to @p x and @p s.
   */
  template <typename Number>
  Number
  jacobi_polynomial_homogenized_derivative(const unsigned int order_x,
                                           const unsigned int order_s,
                                           const unsigned int degree,
                                           const int          alpha,
                                           const int          beta,
                                           const Number       x,
                                           const Number       s);
} // namespace Polynomials
 

 
/* -------------------------- inline functions --------------------- */
 
namespace Polynomials
{
  template <typename number>
  inline Polynomial<number>::Polynomial()
    : in_lagrange_product_form(false)
    , lagrange_weight(1.)
  {}
 
 
 
  template <typename number>
  inline unsigned int
  Polynomial<number>::degree() const
  {
    if (in_lagrange_product_form == true)
      {
        return lagrange_support_points.size();
      }
    else
      {
        Assert(coefficients.size() > 0, ExcEmptyObject());
        return coefficients.size() - 1;
      }
  }
 
 
 
  template <typename number>
  inline number
  Polynomial<number>::value(const number x) const
  {
    if (in_lagrange_product_form == false)
      {
        Assert(coefficients.size() > 0, ExcEmptyObject());
 
        // Horner scheme
        const unsigned int m     = coefficients.size();
        number             value = coefficients.back();
        for (int k = m - 2; k >= 0; --k)
          value = value * x + coefficients[k];
        return value;
      }
    else
      {
        // direct evaluation of Lagrange polynomial
        const unsigned int m     = lagrange_support_points.size();
        number             value = 1.;
        for (unsigned int j = 0; j < m; ++j)
          value *= x - lagrange_support_points[j];
        value *= lagrange_weight;
        return value;
      }
  }
 
 
 
  template <typename number>
  template <typename Number2>
  inline void
  Polynomial<number>::value(const Number2      x,
                            const unsigned int n_derivatives,
                            Number2           *values) const
  {
    values_of_array(std::array<Number2, 1ul>{{x}},
                    n_derivatives,
                    reinterpret_cast<std::array<Number2, 1ul> *>(values));
  }
 
 
 
  template <typename number>
  template <std::size_t n_entries, typename Number2>
  inline
#ifndef DEBUG
    DEAL_II_ALWAYS_INLINE
#endif
    void
    Polynomial<number>::values_of_array(
      const std::array<Number2, n_entries> &x,
      const unsigned int                    n_derivatives,
      std::array<Number2, n_entries>       *values) const
  {
    // evaluate Lagrange polynomial and derivatives
    if (in_lagrange_product_form == true)
      {
        // to compute the value and all derivatives of a polynomial of the
        // form (x-x_1)*(x-x_2)*...*(x-x_n), expand the derivatives like
        // automatic differentiation does.
        const unsigned int n_supp = lagrange_support_points.size();
        const number       weight = lagrange_weight;
        switch (n_derivatives)
          {
            default:
              for (unsigned int e = 0; e < n_entries; ++e)
                values[0][e] = weight;
              for (unsigned int k = 1; k <= n_derivatives; ++k)
                for (unsigned int e = 0; e < n_entries; ++e)
                  values[k][e] = 0.;
              for (unsigned int i = 0; i < n_supp; ++i)
                {
                  std::array<Number2, n_entries> v = x;
                  for (unsigned int e = 0; e < n_entries; ++e)
                    v[e] -= lagrange_support_points[i];
 
                  // multiply by (x-x_i) and compute action on all derivatives,
                  // too (inspired from automatic differentiation: implement the
                  // product rule for the old value and the new variable 'v',
                  // i.e., expand value v and derivative one). since we reuse a
                  // value from the next lower derivative from the steps before,
                  // need to start from the highest derivative
                  for (unsigned int k = n_derivatives; k > 0; --k)
                    for (unsigned int e = 0; e < n_entries; ++e)
                      values[k][e] = (values[k][e] * v[e] + values[k - 1][e]);
                  for (unsigned int e = 0; e < n_entries; ++e)
                    values[0][e] *= v[e];
                }
              // finally, multiply derivatives by k! to transform the product
              // p_n = p^(n)(x)/k! into the actual form of the derivative
              {
                number k_factorial = 2;
                for (unsigned int k = 2; k <= n_derivatives; ++k)
                  {
                    for (unsigned int e = 0; e < n_entries; ++e)
                      values[k][e] *= k_factorial;
                    k_factorial *= static_cast<number>(k + 1);
                  }
              }
              break;
 
            // manually implement case 0 (values only), case 1 (value + first
            // derivative), and case 2 (up to second derivative) since they
            // might be called often. then, we can unroll the inner loop and
            // keep the temporary results as local variables to help the
            // compiler with the pointer aliasing analysis.
            case 0:
              {
                std::array<Number2, n_entries> value;
                for (unsigned int e = 0; e < n_entries; ++e)
                  value[e] = weight;
                for (unsigned int i = 0; i < n_supp; ++i)
                  for (unsigned int e = 0; e < n_entries; ++e)
                    value[e] *= (x[e] - lagrange_support_points[i]);
 
                for (unsigned int e = 0; e < n_entries; ++e)
                  values[0][e] = value[e];
                break;
              }
 
            case 1:
              {
                std::array<Number2, n_entries> value, derivative = {};
                for (unsigned int e = 0; e < n_entries; ++e)
                  value[e] = weight;
                for (unsigned int i = 0; i < n_supp; ++i)
                  for (unsigned int e = 0; e < n_entries; ++e)
                    {
                      const Number2 v = x[e] - lagrange_support_points[i];
                      derivative[e]   = derivative[e] * v + value[e];
                      value[e] *= v;
                    }
 
                for (unsigned int e = 0; e < n_entries; ++e)
                  {
                    values[0][e] = value[e];
                    values[1][e] = derivative[e];
                  }
                break;
              }
 
            case 2:
              {
                std::array<Number2, n_entries> value, derivative = {},
                                                      second = {};
                for (unsigned int e = 0; e < n_entries; ++e)
                  value[e] = weight;
                for (unsigned int i = 0; i < n_supp; ++i)
                  for (unsigned int e = 0; e < n_entries; ++e)
                    {
                      const Number2 v = x[e] - lagrange_support_points[i];
                      second[e]       = second[e] * v + derivative[e];
                      derivative[e]   = derivative[e] * v + value[e];
                      value[e] *= v;
                    }
 
                for (unsigned int e = 0; e < n_entries; ++e)
                  {
                    values[0][e] = value[e];
                    values[1][e] = derivative[e];
                    values[2][e] = static_cast<number>(2) * second[e];
                  }
                break;
              }
          }
        return;
      }
 
    Assert(coefficients.size() > 0, ExcEmptyObject());
 
    // if derivatives are needed, then do it properly by the full
    // Horner scheme
    const unsigned int                          m = coefficients.size();
    std::vector<std::array<Number2, n_entries>> a(coefficients.size());
    for (unsigned int i = 0; i < coefficients.size(); ++i)
      for (unsigned int e = 0; e < n_entries; ++e)
        a[i][e] = coefficients[i];
 
    unsigned int j_factorial = 1;
 
    // loop over all requested derivatives. note that derivatives @p{j>m} are
    // necessarily zero, as they differentiate the polynomial more often than
    // the highest power is
    const unsigned int min_valuessize_m = std::min(n_derivatives + 1, m);
    for (unsigned int j = 0; j < min_valuessize_m; ++j)
      {
        for (int k = m - 2; k >= static_cast<int>(j); --k)
          for (unsigned int e = 0; e < n_entries; ++e)
            a[k][e] += x[e] * a[k + 1][e];
        for (unsigned int e = 0; e < n_entries; ++e)
          values[j][e] = static_cast<number>(j_factorial) * a[j][e];
 
        j_factorial *= j + 1;
      }
 
    // fill higher derivatives by zero
    for (unsigned int j = min_valuessize_m; j <= n_derivatives; ++j)
      for (unsigned int e = 0; e < n_entries; ++e)
        values[j][e] = 0.;
  }
 
 
 
  template <typename number>
  template <class Archive>
  inline void
  Polynomial<number>::serialize(Archive &ar, const unsigned int)
  {
    // forward to serialization function in the base class.
    ar &static_cast<EnableObserverPointer &>(*this);
    ar &coefficients;
    ar &in_lagrange_product_form;
    ar &lagrange_support_points;
    ar &lagrange_weight;
  }
 
 
 
  template <typename Number>
  Number
  jacobi_polynomial_value(const unsigned int degree,
                          const int          alpha,
                          const int          beta,
                          const Number       x,
                          const bool         rescale_to_dealii_unit_interval)
  {
    Assert(alpha >= 0 && beta >= 0,
           ExcNotImplemented("Negative alpha/beta coefficients not supported"));
    // the Jacobi polynomial is evaluated using a recursion formula.
    Number p0, p1;
 
    // The recursion formula is defined for the interval [-1, 1], so rescale
    // to that interval here
    const Number xeval =
      rescale_to_dealii_unit_interval ? Number(-1) + 2. * x : x;
 
    // initial values P_0(x), P_1(x):
    p0 = 1.0;
    if (degree == 0)
      return p0;
    p1 = ((alpha + beta + 2) * xeval + (alpha - beta)) / 2;
    if (degree == 1)
      return p1;
 
    for (unsigned int i = 1; i < degree; ++i)
      {
        const Number v  = 2 * i + (alpha + beta);
        const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
        const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
        const Number a3 = v * (v + 1) * (v + 2);
        const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
 
        const Number pn = ((a2 + a3 * xeval) * p1 - a4 * p0) / a1;
        p0              = p1;
        p1              = pn;
      }
    return p1;
  }
 
 
 
  template <typename Number>
  Number
  jacobi_polynomial_derivative(const unsigned int degree,
                               const int          alpha,
                               const int          beta,
                               const Number       x,
                               const bool rescale_to_dealii_unit_interval)
  {
    return jacobi_polynomial_kth_derivative(
      1, degree, alpha, beta, x, rescale_to_dealii_unit_interval);
  }
 
 
 
  template <typename Number>
  Number
  jacobi_polynomial_kth_derivative(const unsigned int derivative_order,
                                   const unsigned int degree,
                                   const int          alpha,
                                   const int          beta,
                                   const Number       x,
                                   const bool rescale_to_dealii_unit_interval)
  {
    Assert(alpha >= 0 && beta >= 0,
           ExcNotImplemented("Negative alpha/beta coefficients not supported"));
 
    if (derivative_order > degree)
      return 0.0;
    // P_0 = 1, so the derivative is 0
    // special case when derivative_order = 0, then the value should be returned
    // and not 0
    if (degree == 0 && derivative_order != 0)
      return 0.0;
 
    // The derivative of the Jacobi polynomial is evaluated using the recurrence
    // relations
    Number pre_factor = 1.0;
    for (unsigned int i = 1; i < derivative_order + 1; ++i)
      pre_factor *= (alpha + beta + degree + i);
 
    if (rescale_to_dealii_unit_interval)
      return pre_factor * jacobi_polynomial_value(degree - derivative_order,
                                                  alpha + derivative_order,
                                                  beta + derivative_order,
                                                  x,
                                                  true);
 
    return std::pow(0.5, derivative_order) * pre_factor *
           jacobi_polynomial_value(degree - derivative_order,
                                   alpha + derivative_order,
                                   beta + derivative_order,
                                   x,
                                   false);
  }
 
 
 
  template <typename Number>
  std::vector<Number>
  jacobi_polynomial_roots(const unsigned int degree,
                          const int          alpha,
                          const int          beta)
  {
    std::vector<Number> x(degree, 0.5);
 
    // compute zeros with a Newton algorithm.
 
    // Set tolerance. For long double we might not always get the additional
    // precision in a run time environment (e.g. with valgrind), so we must
    // limit the tolerance to double. Since we do a Newton iteration, doing
    // one more iteration after the residual has indicated convergence will be
    // enough for all number types due to the quadratic convergence of
    // Newton's method
 
    const Number tolerance =
      4 * std::max(static_cast<Number>(std::numeric_limits<double>::epsilon()),
                   std::numeric_limits<Number>::epsilon());
 
    // The following implementation follows closely the one given in the
    // appendix of the book by Karniadakis and Sherwin: Spectral/hp element
    // methods for computational fluid dynamics (Oxford University Press,
    // 2005)
 
    // If symmetric, we only need to compute the half of points
    const unsigned int n_points = (alpha == beta ? degree / 2 : degree);
    for (unsigned int k = 0; k < n_points; ++k)
      {
        // we take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as
        // initial values, corrected by the initial value
        Number r = 0.5 - 0.5 * std::cos(static_cast<Number>(2 * k + 1) /
                                        (2 * degree) * numbers::PI);
        if (k > 0)
          r = (r + x[k - 1]) / 2;
 
        unsigned int converged = numbers::invalid_unsigned_int;
        for (unsigned int it = 1; it < 1000; ++it)
          {
            Number s = 0.;
            for (unsigned int i = 0; i < k; ++i)
              s += 1. / (r - x[i]);
 
            // derivative of P_n^{alpha,beta}, rescaled to [0, 1]
            const Number J_x =
              (alpha + beta + degree + 1) *
              jacobi_polynomial_value(degree - 1, alpha + 1, beta + 1, r);
 
            // value of P_n^{alpha,beta}
            const Number f = jacobi_polynomial_value(degree, alpha, beta, r);
            const Number delta = f / (f * s - J_x);
            r += delta;
            if (converged == numbers::invalid_unsigned_int &&
                std::abs(delta) < tolerance)
              converged = it;
 
            // do one more iteration to ensure accuracy also for tighter
            // types than double (e.g. long double)
            if (it == converged + 1)
              break;
          }
 
        Assert(converged != numbers::invalid_unsigned_int,
               ExcMessage("Newton iteration for zero of Jacobi polynomial "
                          "did not converge."));
 
        x[k] = r;
      }
 
    // in case we assumed symmetry, fill up the missing values
    for (unsigned int k = n_points; k < degree; ++k)
      x[k] = 1.0 - x[degree - k - 1];
 
    return x;
  }
 
 
 
  template <typename Number>
  Number
  jacobi_polynomial_homogenized_value(const unsigned int degree,
                                      const int          alpha,
                                      const int          beta,
                                      const Number       x,
                                      const Number       s)
  {
    Assert(alpha >= 0 && beta >= 0,
           ExcNotImplemented("Negative alpha/beta coefficients not supported"));
 
    // the homogenized Jacobi polynomial is evaluated using a recursion formula
    // to get the recursion formula, the recursion rule for the Jacobi
    // polynomial is multiplied on both sides by s^degree
    Number p0, p1;
 
    // initial values p0 = s^0 * P_0(2x/s-1) = 1 * 1 = 1, p1 = s * P_1(2x/s-1):
    p0 = 1.0;
    if (degree == 0)
      return p0;
    p1 = ((alpha + beta + 2.0) * (2.0 * x - s) + s * (alpha - beta)) / 2.0;
    if (degree == 1)
      return p1;
 
    for (unsigned int i = 1; i < degree; ++i)
      {
        const Number v  = 2 * i + (alpha + beta);
        const Number a1 = 2 * (i + 1) * (i + (alpha + beta + 1)) * v;
        const Number a2 = (v + 1) * (alpha * alpha - beta * beta);
        const Number a3 = v * (v + 1) * (v + 2);
        const Number a4 = 2 * (i + alpha) * (i + beta) * (v + 2);
 
        const Number pn =
          ((a2 * s + a3 * (2.0 * x - s)) * p1 - a4 * s * s * p0) / a1;
        p0 = p1;
        p1 = pn;
      }
    return p1;
  }
 
 
 
  template <typename Number>
  Number
  jacobi_polynomial_homogenized_derivative(const unsigned int order_x,
                                           const unsigned int order_s,
                                           const unsigned int degree,
                                           const int          alpha,
                                           const int          beta,
                                           const Number       x,
                                           const Number       s)
  {
    Assert(alpha >= 0 && beta >= 0,
           ExcNotImplemented("Negative alpha/beta coefficients not supported"));
 
    if (order_x + order_s > degree)
      return 0.0;
    if (degree == 0 && (order_x + order_s != 0))
      return 0.0;
 
    // The derivative of the homogenized Jacobi polynomial is evaluated using
    // the recurrence relations
    Number pre_factor = 1.0;
    for (unsigned int i = 1; i < order_x + 1; ++i)
      pre_factor *= (alpha + beta + degree + i);
 
    for (unsigned int i = 0; i < order_s; ++i)
      pre_factor *= (degree + beta - i);
 
    const Number derivative =
      std::pow(-1.0, order_s) * pre_factor *
      jacobi_polynomial_homogenized_value<Number>(degree - order_x - order_s,
                                                  alpha + order_x + order_s,
                                                  beta + order_x,
                                                  x,
                                                  s);
 
    return derivative;
  }
} // namespace Polynomials
DEAL_II_NAMESPACE_CLOSE
 
#endif