// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2000 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------

#include <deal.II/base/exceptions.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/point.h>
#include <deal.II/base/polynomial.h>
#include <deal.II/base/quadrature_lib.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/utilities.h>

#include <algorithm>
#include <cmath>
#include <limits>
#include <shared_mutex>

DEAL_II_NAMESPACE_OPEN



namespace Polynomials
{
  // -------------------- class Polynomial ---------------- //


  template <typename number>
  Polynomial<number>::Polynomial(const std::vector<number> &a)
    : coefficients(a)
    , in_lagrange_product_form(false)
    , lagrange_weight(1.)
  {}



  template <typename number>
  Polynomial<number>::Polynomial(const unsigned int n)
    : coefficients(n + 1, 0.)
    , in_lagrange_product_form(false)
    , lagrange_weight(1.)
  {}



  template <typename number>
  Polynomial<number>::Polynomial(const std::vector<Point<1>> &supp,
                                 const unsigned int           center)
    : in_lagrange_product_form(true)
  {
    Assert(supp.size() > 0, ExcEmptyObject());
    AssertIndexRange(center, supp.size());

    lagrange_support_points.reserve(supp.size() - 1);
    number tmp_lagrange_weight = 1.;
    for (unsigned int i = 0; i < supp.size(); ++i)
      if (i != center)
        {
          lagrange_support_points.push_back(supp[i][0]);
          tmp_lagrange_weight *= supp[center][0] - supp[i][0];
        }

    // check for underflow and overflow
    Assert(std::fabs(tmp_lagrange_weight) > std::numeric_limits<number>::min(),
           ExcMessage("Underflow in computation of Lagrange denominator."));
    Assert(std::fabs(tmp_lagrange_weight) < std::numeric_limits<number>::max(),
           ExcMessage("Overflow in computation of Lagrange denominator."));

    lagrange_weight = static_cast<number>(1.) / tmp_lagrange_weight;
  }



  template <typename number>
  void
  Polynomial<number>::value(const number x, std::vector<number> &values) const
  {
    Assert(values.size() > 0, ExcZero());

    value(x, values.size() - 1, values.data());
  }



  template <typename number>
  void
  Polynomial<number>::transform_into_standard_form()
  {
    // should only be called when the product form is active
    Assert(in_lagrange_product_form == true, ExcInternalError());
    Assert(coefficients.empty(), ExcInternalError());

    // compute coefficients by expanding the product (x-x_i) term by term
    coefficients.resize(lagrange_support_points.size() + 1);
    if (lagrange_support_points.empty())
      coefficients[0] = 1.;
    else
      {
        coefficients[0] = -lagrange_support_points[0];
        coefficients[1] = 1.;
        for (unsigned int i = 1; i < lagrange_support_points.size(); ++i)
          {
            coefficients[i + 1] = 1.;
            for (unsigned int j = i; j > 0; --j)
              coefficients[j] = (-lagrange_support_points[i] * coefficients[j] +
                                 coefficients[j - 1]);
            coefficients[0] *= -lagrange_support_points[i];
          }
      }
    for (unsigned int i = 0; i < lagrange_support_points.size() + 1; ++i)
      coefficients[i] *= lagrange_weight;

    // delete the product form data
    std::vector<number> new_points;
    lagrange_support_points.swap(new_points);
    in_lagrange_product_form = false;
    lagrange_weight          = 1.;
  }



  template <typename number>
  void
  Polynomial<number>::scale(std::vector<number> &coefficients,
                            const number         factor)
  {
    number f = 1.;
    for (typename std::vector<number>::iterator c = coefficients.begin();
         c != coefficients.end();
         ++c)
      {
        *c *= f;
        f *= factor;
      }
  }



  template <typename number>
  void
  Polynomial<number>::scale(const number factor)
  {
    // to scale (x-x_0)*(x-x_1)*...*(x-x_n), scale
    // support points by 1./factor and the weight
    // likewise
    if (in_lagrange_product_form == true)
      {
        number inv_fact         = number(1.) / factor;
        number accumulated_fact = 1.;
        for (unsigned int i = 0; i < lagrange_support_points.size(); ++i)
          {
            lagrange_support_points[i] *= inv_fact;
            accumulated_fact *= factor;
          }
        lagrange_weight *= accumulated_fact;
      }
    // otherwise, use the function above
    else
      scale(coefficients, factor);
  }



  template <typename number>
  void
  Polynomial<number>::multiply(std::vector<number> &coefficients,
                               const number         factor)
  {
    for (typename std::vector<number>::iterator c = coefficients.begin();
         c != coefficients.end();
         ++c)
      *c *= factor;
  }



  template <typename number>
  Polynomial<number> &
  Polynomial<number>::operator*=(const double s)
  {
    if (in_lagrange_product_form == true)
      lagrange_weight *= s;
    else
      {
        for (typename std::vector<number>::iterator c = coefficients.begin();
             c != coefficients.end();
             ++c)
          *c *= s;
      }
    return *this;
  }



  template <typename number>
  Polynomial<number> &
  Polynomial<number>::operator*=(const Polynomial<number> &p)
  {
    // if we are in Lagrange form, just append the
    // new points
    if (in_lagrange_product_form == true && p.in_lagrange_product_form == true)
      {
        lagrange_weight *= p.lagrange_weight;
        lagrange_support_points.insert(lagrange_support_points.end(),
                                       p.lagrange_support_points.begin(),
                                       p.lagrange_support_points.end());
      }

    // cannot retain product form, recompute...
    else if (in_lagrange_product_form == true)
      transform_into_standard_form();

    // need to transform p into standard form as
    // well if necessary. copy the polynomial to
    // do this
    std::unique_ptr<Polynomial<number>> q_data;
    const Polynomial<number>           *q = nullptr;
    if (p.in_lagrange_product_form == true)
      {
        q_data = std::make_unique<Polynomial<number>>(p);
        q_data->transform_into_standard_form();
        q = q_data.get();
      }
    else
      q = &p;

    // Degree of the product
    unsigned int new_degree = this->degree() + q->degree();

    std::vector<number> new_coefficients(new_degree + 1, 0.);

    for (unsigned int i = 0; i < q->coefficients.size(); ++i)
      for (unsigned int j = 0; j < this->coefficients.size(); ++j)
        new_coefficients[i + j] += this->coefficients[j] * q->coefficients[i];
    this->coefficients = new_coefficients;

    return *this;
  }



  template <typename number>
  Polynomial<number> &
  Polynomial<number>::operator+=(const Polynomial<number> &p)
  {
    // Lagrange product form cannot reasonably be
    // retained after polynomial addition. we
    // could in theory check if either this
    // polynomial or the other is a zero
    // polynomial and retain it, but we actually
    // currently (r23974) assume that the addition
    // of a zero polynomial changes the state and
    // tests equivalence.
    if (in_lagrange_product_form == true)
      transform_into_standard_form();

    // need to transform p into standard form as
    // well if necessary. copy the polynomial to
    // do this
    std::unique_ptr<Polynomial<number>> q_data;
    const Polynomial<number>           *q = nullptr;
    if (p.in_lagrange_product_form == true)
      {
        q_data = std::make_unique<Polynomial<number>>(p);
        q_data->transform_into_standard_form();
        q = q_data.get();
      }
    else
      q = &p;

    // if necessary expand the number
    // of coefficients we store
    if (q->coefficients.size() > coefficients.size())
      coefficients.resize(q->coefficients.size(), 0.);

    for (unsigned int i = 0; i < q->coefficients.size(); ++i)
      coefficients[i] += q->coefficients[i];

    return *this;
  }



  template <typename number>
  Polynomial<number> &
  Polynomial<number>::operator-=(const Polynomial<number> &p)
  {
    // Lagrange product form cannot reasonably be
    // retained after polynomial addition
    if (in_lagrange_product_form == true)
      transform_into_standard_form();

    // need to transform p into standard form as
    // well if necessary. copy the polynomial to
    // do this
    std::unique_ptr<Polynomial<number>> q_data;
    const Polynomial<number>           *q = nullptr;
    if (p.in_lagrange_product_form == true)
      {
        q_data = std::make_unique<Polynomial<number>>(p);
        q_data->transform_into_standard_form();
        q = q_data.get();
      }
    else
      q = &p;

    // if necessary expand the number
    // of coefficients we store
    if (q->coefficients.size() > coefficients.size())
      coefficients.resize(q->coefficients.size(), 0.);

    for (unsigned int i = 0; i < q->coefficients.size(); ++i)
      coefficients[i] -= q->coefficients[i];

    return *this;
  }



  template <typename number>
  bool
  Polynomial<number>::operator==(const Polynomial<number> &p) const
  {
    // need to distinguish a few cases based on
    // whether we are in product form or not. two
    // polynomials can still be the same when they
    // are on different forms, but the expansion
    // is the same
    if (in_lagrange_product_form == true && p.in_lagrange_product_form == true)
      return ((lagrange_weight == p.lagrange_weight) &&
              (lagrange_support_points == p.lagrange_support_points));
    else if (in_lagrange_product_form == true)
      {
        Polynomial<number> q = *this;
        q.transform_into_standard_form();
        return (q.coefficients == p.coefficients);
      }
    else if (p.in_lagrange_product_form == true)
      {
        Polynomial<number> q = p;
        q.transform_into_standard_form();
        return (q.coefficients == coefficients);
      }
    else
      return (p.coefficients == coefficients);
  }



  template <typename number>
  template <typename number2>
  void
  Polynomial<number>::shift(std::vector<number> &coefficients,
                            const number2        offset)
  {
    // too many coefficients cause overflow in
    // the binomial coefficient used below
    Assert(coefficients.size() < 31, ExcNotImplemented());

    // Copy coefficients to a vector of
    // accuracy given by the argument
    std::vector<number2> new_coefficients(coefficients.begin(),
                                          coefficients.end());

    // Traverse all coefficients from
    // c_1. c_0 will be modified by
    // higher degrees, only.
    for (unsigned int d = 1; d < new_coefficients.size(); ++d)
      {
        const unsigned int n = d;
        // Binomial coefficients are
        // needed for the
        // computation. The rightmost
        // value is unity.
        unsigned int binomial_coefficient = 1;

        // Powers of the offset will be
        // needed and computed
        // successively.
        number2 offset_power = offset;

        // Compute (x+offset)^d
        // and modify all values c_k
        // with k<d.
        // The coefficient in front of
        // x^d is not modified in this step.
        for (unsigned int k = 0; k < d; ++k)
          {
            // Recursion from Bronstein
            // Make sure no remainders
            // occur in integer
            // division.
            binomial_coefficient = (binomial_coefficient * (n - k)) / (k + 1);

            new_coefficients[d - k - 1] +=
              new_coefficients[d] * binomial_coefficient * offset_power;
            offset_power *= offset;
          }
        // The binomial coefficient
        // should have gone through a
        // whole row of Pascal's
        // triangle.
        Assert(binomial_coefficient == 1, ExcInternalError());
      }

    // copy new elements to old vector
    coefficients.assign(new_coefficients.begin(), new_coefficients.end());
  }



  template <typename number>
  template <typename number2>
  void
  Polynomial<number>::shift(const number2 offset)
  {
    // shift is simple for a polynomial in product
    // form, (x-x_0)*(x-x_1)*...*(x-x_n). just add
    // offset to all shifts
    if (in_lagrange_product_form == true)
      {
        for (unsigned int i = 0; i < lagrange_support_points.size(); ++i)
          lagrange_support_points[i] -= offset;
      }
    else
      // do the shift in any case
      shift(coefficients, offset);
  }



  template <typename number>
  Polynomial<number>
  Polynomial<number>::derivative() const
  {
    // no simple form possible for Lagrange
    // polynomial on product form
    if (degree() == 0)
      return Monomial<number>(0, 0.);

    std::unique_ptr<Polynomial<number>> q_data;
    const Polynomial<number>           *q = nullptr;
    if (in_lagrange_product_form == true)
      {
        q_data = std::make_unique<Polynomial<number>>(*this);
        q_data->transform_into_standard_form();
        q = q_data.get();
      }
    else
      q = this;

    std::vector<number> newcoefficients(q->coefficients.size() - 1);
    for (unsigned int i = 1; i < q->coefficients.size(); ++i)
      newcoefficients[i - 1] = number(i) * q->coefficients[i];

    return Polynomial<number>(newcoefficients);
  }



  template <typename number>
  Polynomial<number>
  Polynomial<number>::primitive() const
  {
    // no simple form possible for Lagrange
    // polynomial on product form
    std::unique_ptr<Polynomial<number>> q_data;
    const Polynomial<number>           *q = nullptr;
    if (in_lagrange_product_form == true)
      {
        q_data = std::make_unique<Polynomial<number>>(*this);
        q_data->transform_into_standard_form();
        q = q_data.get();
      }
    else
      q = this;

    std::vector<number> newcoefficients(q->coefficients.size() + 1);
    newcoefficients[0] = 0.;
    for (unsigned int i = 0; i < q->coefficients.size(); ++i)
      newcoefficients[i + 1] = q->coefficients[i] / number(i + 1.);

    return Polynomial<number>(newcoefficients);
  }



  template <typename number>
  void
  Polynomial<number>::print(std::ostream &out) const
  {
    if (in_lagrange_product_form == true)
      {
        out << lagrange_weight;
        for (unsigned int i = 0; i < lagrange_support_points.size(); ++i)
          out << " (x-" << lagrange_support_points[i] << ")";
        out << std::endl;
      }
    else
      for (int i = degree(); i >= 0; --i)
        {
          out << coefficients[i] << " x^" << i << std::endl;
        }
  }


  template <typename number>
  std::size_t
  Polynomial<number>::memory_consumption() const
  {
    return (MemoryConsumption::memory_consumption(coefficients) +
            MemoryConsumption::memory_consumption(in_lagrange_product_form) +
            MemoryConsumption::memory_consumption(lagrange_support_points) +
            MemoryConsumption::memory_consumption(lagrange_weight));
  }



  // ------------------ class Monomial -------------------------- //

  template <typename number>
  std::vector<number>
  Monomial<number>::make_vector(unsigned int n, double coefficient)
  {
    std::vector<number> result(n + 1, 0.);
    result[n] = coefficient;
    return result;
  }



  template <typename number>
  Monomial<number>::Monomial(unsigned int n, double coefficient)
    : Polynomial<number>(make_vector(n, coefficient))
  {}



  template <typename number>
  std::vector<Polynomial<number>>
  Monomial<number>::generate_complete_basis(const unsigned int degree)
  {
    std::vector<Polynomial<number>> v;
    v.reserve(degree + 1);
    for (unsigned int i = 0; i <= degree; ++i)
      v.push_back(Monomial<number>(i));
    return v;
  }



  // ------------------ class LagrangeEquidistant --------------- //

  namespace internal
  {
    namespace LagrangeEquidistantImplementation
    {
      std::vector<Point<1>>
      generate_equidistant_unit_points(const unsigned int n)
      {
        std::vector<Point<1>> points(n + 1);
        const double          one_over_n = 1. / n;
        for (unsigned int k = 0; k <= n; ++k)
          points[k][0] = static_cast<double>(k) * one_over_n;
        return points;
      }
    } // namespace LagrangeEquidistantImplementation
  }   // namespace internal



  LagrangeEquidistant::LagrangeEquidistant(const unsigned int n,
                                           const unsigned int support_point)
    : Polynomial<double>(internal::LagrangeEquidistantImplementation::
                           generate_equidistant_unit_points(n),
                         support_point)
  {
    Assert(coefficients.empty(), ExcInternalError());

    // For polynomial order up to 3, we have precomputed weights. Use these
    // weights instead of the product form
    if (n <= 3)
      {
        this->in_lagrange_product_form = false;
        this->lagrange_weight          = 1.;
        std::vector<double> new_support_points;
        this->lagrange_support_points.swap(new_support_points);
        this->coefficients.resize(n + 1);
        compute_coefficients(n, support_point, this->coefficients);
      }
  }



  void
  LagrangeEquidistant::compute_coefficients(const unsigned int   n,
                                            const unsigned int   support_point,
                                            std::vector<double> &a)
  {
    AssertIndexRange(support_point, n + 1);

    unsigned int n_functions = n + 1;
    AssertIndexRange(support_point, n_functions);
    const double *x = nullptr;

    switch (n)
      {
        case 1:
          {
            static const double x1[4] = {1.0, -1.0, 0.0, 1.0};
            x                         = &x1[0];
            break;
          }
        case 2:
          {
            static const double x2[9] = {
              1.0, -3.0, 2.0, 0.0, 4.0, -4.0, 0.0, -1.0, 2.0};
            x = &x2[0];
            break;
          }
        case 3:
          {
            static const double x3[16] = {1.0,
                                          -11.0 / 2.0,
                                          9.0,
                                          -9.0 / 2.0,
                                          0.0,
                                          9.0,
                                          -45.0 / 2.0,
                                          27.0 / 2.0,
                                          0.0,
                                          -9.0 / 2.0,
                                          18.0,
                                          -27.0 / 2.0,
                                          0.0,
                                          1.0,
                                          -9.0 / 2.0,
                                          9.0 / 2.0};
            x                          = &x3[0];
            break;
          }
        default:
          DEAL_II_ASSERT_UNREACHABLE();
      }

    Assert(x != nullptr, ExcInternalError());
    for (unsigned int i = 0; i < n_functions; ++i)
      a[i] = x[support_point * n_functions + i];
  }



  std::vector<Polynomial<double>>
  LagrangeEquidistant::generate_complete_basis(const unsigned int degree)
  {
    if (degree == 0)
      // create constant polynomial
      return std::vector<Polynomial<double>>(
        1, Polynomial<double>(std::vector<double>(1, 1.)));
    else
      {
        // create array of Lagrange
        // polynomials
        std::vector<Polynomial<double>> v;
        for (unsigned int i = 0; i <= degree; ++i)
          v.push_back(LagrangeEquidistant(degree, i));
        return v;
      }
  }



  //----------------------------------------------------------------------//


  std::vector<Polynomial<double>>
  generate_complete_Lagrange_basis(const std::vector<Point<1>> &points)
  {
    std::vector<Polynomial<double>> p;
    p.reserve(points.size());

    for (unsigned int i = 0; i < points.size(); ++i)
      p.emplace_back(points, i);
    return p;
  }



  // ------------------ class Legendre --------------- //



  Legendre::Legendre(const unsigned int k)
    : Polynomial<double>(0)
  {
    this->coefficients.clear();
    this->in_lagrange_product_form = true;
    this->lagrange_support_points.resize(k);

    // the roots of a Legendre polynomial are exactly the points in the
    // Gauss-Legendre quadrature formula
    if (k > 0)
      {
        const QGauss<1> gauss(k);
        for (unsigned int i = 0; i < k; ++i)
          this->lagrange_support_points[i] = gauss.get_points()[i][0];
      }

    // compute the abscissa in zero of the product of monomials. The exact
    // value should be sqrt(2*k+1), so set the weight to that value.
    double prod = 1.;
    for (unsigned int i = 0; i < k; ++i)
      prod *= this->lagrange_support_points[i];
    this->lagrange_weight = std::sqrt(double(2 * k + 1)) / prod;
  }



  std::vector<Polynomial<double>>
  Legendre::generate_complete_basis(const unsigned int degree)
  {
    std::vector<Polynomial<double>> v;
    v.reserve(degree + 1);
    for (unsigned int i = 0; i <= degree; ++i)
      v.push_back(Legendre(i));
    return v;
  }



  // ------------------ class Lobatto -------------------- //


  Lobatto::Lobatto(const unsigned int p)
    : Polynomial<double>(compute_coefficients(p))
  {}

  std::vector<double>
  Lobatto::compute_coefficients(const unsigned int p)
  {
    switch (p)
      {
        case 0:
          {
            std::vector<double> coefficients(2);

            coefficients[0] = 1.0;
            coefficients[1] = -1.0;
            return coefficients;
          }

        case 1:
          {
            std::vector<double> coefficients(2);

            coefficients[0] = 0.0;
            coefficients[1] = 1.0;
            return coefficients;
          }

        case 2:
          {
            std::vector<double> coefficients(3);

            coefficients[0] = 0.0;
            coefficients[1] = -1.0 * std::sqrt(3.);
            coefficients[2] = std::sqrt(3.);
            return coefficients;
          }

        default:
          {
            std::vector<double> coefficients(p + 1);
            std::vector<double> legendre_coefficients_tmp1(p);
            std::vector<double> legendre_coefficients_tmp2(p - 1);

            coefficients[0]               = -1.0 * std::sqrt(3.);
            coefficients[1]               = 2.0 * std::sqrt(3.);
            legendre_coefficients_tmp1[0] = 1.0;

            for (unsigned int i = 2; i < p; ++i)
              {
                for (unsigned int j = 0; j < i - 1; ++j)
                  legendre_coefficients_tmp2[j] = legendre_coefficients_tmp1[j];

                for (unsigned int j = 0; j < i; ++j)
                  legendre_coefficients_tmp1[j] = coefficients[j];

                coefficients[0] =
                  std::sqrt(2 * i + 1.) *
                  ((1.0 - 2 * i) * legendre_coefficients_tmp1[0] /
                     std::sqrt(2 * i - 1.) +
                   (1.0 - i) * legendre_coefficients_tmp2[0] /
                     std::sqrt(2 * i - 3.)) /
                  i;

                for (unsigned int j = 1; j < i - 1; ++j)
                  coefficients[j] =
                    std::sqrt(2 * i + 1.) *
                    (std::sqrt(2 * i - 1.) *
                       (2.0 * legendre_coefficients_tmp1[j - 1] -
                        legendre_coefficients_tmp1[j]) +
                     (1.0 - i) * legendre_coefficients_tmp2[j] /
                       std::sqrt(2 * i - 3.)) /
                    i;

                coefficients[i - 1] = std::sqrt(4 * i * i - 1.) *
                                      (2.0 * legendre_coefficients_tmp1[i - 2] -
                                       legendre_coefficients_tmp1[i - 1]) /
                                      i;
                coefficients[i] = 2.0 * std::sqrt(4 * i * i - 1.) *
                                  legendre_coefficients_tmp1[i - 1] / i;
              }

            for (int i = p; i > 0; --i)
              coefficients[i] = coefficients[i - 1] / i;

            coefficients[0] = 0.0;
            return coefficients;
          }
      }
  }

  std::vector<Polynomial<double>>
  Lobatto::generate_complete_basis(const unsigned int p)
  {
    std::vector<Polynomial<double>> basis(p + 1);

    for (unsigned int i = 0; i <= p; ++i)
      basis[i] = Lobatto(i);

    return basis;
  }



  // ------------------ class Hierarchical --------------- //

  // Reserve space for polynomials up to degree 19. Should be sufficient
  // for the start.
  std::vector<std::unique_ptr<const std::vector<double>>>
    Hierarchical::recursive_coefficients(20);

  std::shared_mutex Hierarchical::coefficients_lock;


  Hierarchical::Hierarchical(const unsigned int k)
    : Polynomial<double>(get_coefficients(k))
  {}



  void
  Hierarchical::compute_coefficients(const unsigned int k_)
  {
    unsigned int k = k_;
    // The first 2 coefficients
    // are hard-coded
    if (k == 0)
      k = 1;

    // First see whether the coefficients we need have already been
    // computed. This is a read operation, and so we can do that
    // with a shared lock.
    //
    // (We could have gotten away without any lock at all if the
    // inner pointers were std::atomic<std::unique_ptr<...>>, but
    // first, there is no such specialization of std::atomic that
    // is mutex-free, and then there is also the issue that the
    // outer vector may be resized and that can definitely not
    // be guarded against without a mutex of some sort.)
    {
      std::shared_lock<std::shared_mutex> lock(coefficients_lock);

      if ((recursive_coefficients.size() >= k + 1) &&
          (recursive_coefficients[k].get() != nullptr))
        return;
    }

    // Having gotten here, we know that we need to compute a new set
    // of coefficients. This has to happen under a unique lock because
    // we're not only reading, but writing into the data structures:
    std::unique_lock<std::shared_mutex> lock(coefficients_lock);

    // First make sure that there is enough
    // space in the array for the
    // coefficients, so we have to resize
    // it to size k+1

    // but it's more complicated than
    // that: we call this function
    // recursively, so if we simply
    // resize it to k+1 here, then
    // compute the coefficients for
    // degree k-1 by calling this
    // function recursively, then it will
    // reset the size to k -- not enough
    // for what we want to do below. the
    // solution therefore is to only
    // resize the size if we are going to
    // *increase* it
    if (recursive_coefficients.size() < k + 1)
      recursive_coefficients.resize(k + 1);

    if (k <= 1)
      {
        // create coefficients
        // vectors for k=0 and k=1
        //
        // allocate the respective
        // amount of memory and
        // later assign it to the
        // coefficients array to
        // make it const
        std::vector<double> c0(2);
        c0[0] = 1.;
        c0[1] = -1.;

        std::vector<double> c1(2);
        c1[0] = 0.;
        c1[1] = 1.;

        // now make these arrays
        // const
        recursive_coefficients[0] =
          std::make_unique<const std::vector<double>>(std::move(c0));
        recursive_coefficients[1] =
          std::make_unique<const std::vector<double>>(std::move(c1));
      }
    else if (k == 2)
      {
        coefficients_lock.unlock();
        compute_coefficients(1);
        coefficients_lock.lock();

        std::vector<double> c2(3);

        const double a = 1.; // 1./8.;

        c2[0] = 0. * a;
        c2[1] = -4. * a;
        c2[2] = 4. * a;

        recursive_coefficients[2] =
          std::make_unique<const std::vector<double>>(std::move(c2));
      }
    else
      {
        // for larger numbers,
        // compute the coefficients
        // recursively. to do so,
        // we have to release the
        // lock temporarily to
        // allow the called
        // function to acquire it
        // itself
        coefficients_lock.unlock();
        compute_coefficients(k - 1);
        coefficients_lock.lock();

        std::vector<double> ck(k + 1);

        const double a = 1.; // 1./(2.*k);

        ck[0] = -a * (*recursive_coefficients[k - 1])[0];

        for (unsigned int i = 1; i <= k - 1; ++i)
          ck[i] = a * (2. * (*recursive_coefficients[k - 1])[i - 1] -
                       (*recursive_coefficients[k - 1])[i]);

        ck[k] = a * 2. * (*recursive_coefficients[k - 1])[k - 1];
        // for even degrees, we need
        // to add a multiple of
        // basis fcn phi_2
        if ((k % 2) == 0)
          {
            double b = 1.; // 8.;
            // for (unsigned int i=1; i<=k; ++i)
            //  b /= 2.*i;

            ck[1] += b * (*recursive_coefficients[2])[1];
            ck[2] += b * (*recursive_coefficients[2])[2];
          }
        // finally assign the newly
        // created vector to the
        // const pointer in the
        // coefficients array
        recursive_coefficients[k] =
          std::make_unique<const std::vector<double>>(std::move(ck));
      }
  }



  const std::vector<double> &
  Hierarchical::get_coefficients(const unsigned int k)
  {
    // First make sure the coefficients get computed if so necessary
    compute_coefficients(k);

    // Then get a pointer to the array of coefficients. Do that in a MT
    // safe way, but since we're only reading information we can do
    // that with a shared lock
    std::shared_lock<std::shared_mutex> lock(coefficients_lock);
    return *recursive_coefficients[k];
  }



  std::vector<Polynomial<double>>
  Hierarchical::generate_complete_basis(const unsigned int degree)
  {
    if (degree == 0)
      // create constant
      // polynomial. note that we
      // can't use the other branch
      // of the if-statement, since
      // calling the constructor of
      // this class with argument
      // zero does _not_ create the
      // constant polynomial, but
      // rather 1-x
      return std::vector<Polynomial<double>>(
        1, Polynomial<double>(std::vector<double>(1, 1.)));
    else
      {
        std::vector<Polynomial<double>> v;
        v.reserve(degree + 1);
        for (unsigned int i = 0; i <= degree; ++i)
          v.push_back(Hierarchical(i));
        return v;
      }
  }



  // ------------------ HermiteInterpolation --------------- //

  HermiteInterpolation::HermiteInterpolation(const unsigned int p)
    : Polynomial<double>(0)
  {
    this->coefficients.clear();
    this->in_lagrange_product_form = true;

    this->lagrange_support_points.resize(3);
    if (p == 0)
      {
        this->lagrange_support_points[0] = -0.5;
        this->lagrange_support_points[1] = 1.;
        this->lagrange_support_points[2] = 1.;
        this->lagrange_weight            = 2.;
      }
    else if (p == 1)
      {
        this->lagrange_support_points[0] = 0.;
        this->lagrange_support_points[1] = 0.;
        this->lagrange_support_points[2] = 1.5;
        this->lagrange_weight            = -2.;
      }
    else if (p == 2)
      {
        this->lagrange_support_points[0] = 0.;
        this->lagrange_support_points[1] = 1.;
        this->lagrange_support_points[2] = 1.;
      }
    else if (p == 3)
      {
        this->lagrange_support_points[0] = 0.;
        this->lagrange_support_points[1] = 0.;
        this->lagrange_support_points[2] = 1.;
      }
    else
      {
        this->lagrange_support_points.resize(4);
        this->lagrange_support_points[0] = 0.;
        this->lagrange_support_points[1] = 0.;
        this->lagrange_support_points[2] = 1.;
        this->lagrange_support_points[3] = 1.;
        this->lagrange_weight            = 16.;

        if (p > 4)
          {
            Legendre legendre(p - 4);
            (*this) *= legendre;
          }
      }
  }


  std::vector<Polynomial<double>>
  HermiteInterpolation::generate_complete_basis(const unsigned int n)
  {
    Assert(n >= 3,
           ExcNotImplemented("Hermite interpolation makes no sense for "
                             "degrees less than three"));
    std::vector<Polynomial<double>> basis(n + 1);

    for (unsigned int i = 0; i <= n; ++i)
      basis[i] = HermiteInterpolation(i);

    return basis;
  }


  // ------------------ HermiteLikeInterpolation --------------- //
  namespace
  {
    // Finds the zero position x_star such that the mass matrix entry (0,1)
    // with the Hermite polynomials evaluates to zero. The function has
    // originally been derived by a secant method for the integral entry
    // l_0(x) * l_1(x) but we only need to do one iteration because the zero
    // x_star is linear in the integral value.
    double
    find_support_point_x_star(const std::vector<double> &jacobi_roots)
    {
      // Initial guess for the support point position values: The zero turns
      // out to be between zero and the first root of the Jacobi polynomial,
      // but the algorithm is agnostic about that, so simply choose two points
      // that are sufficiently far apart.
      double             guess_left  = 0;
      double             guess_right = 0.5;
      const unsigned int degree      = jacobi_roots.size() + 3;

      // Compute two integrals of the product of l_0(x) * l_1(x)
      // l_0(x) =
      // (x-y)*(x-jacobi_roots(0))*...*(x-jacobi_roos(degree-4))*(x-1)*(x-1)
      // l_1(x) =
      // (x-0)*(x-jacobi_roots(0))*...*(x-jacobi_roots(degree-4))*(x-1)*(x-1)
      // where y is either guess_left or guess_right for the two integrals.
      // Note that the polynomials are not yet normalized here, which is not
      // necessary because we are only looking for the x_star where the matrix
      // entry is zero, for which the constants do not matter.
      const QGauss<1> gauss(degree + 1);
      double          integral_left = 0, integral_right = 0;
      for (unsigned int q = 0; q < gauss.size(); ++q)
        {
          const double x               = gauss.point(q)[0];
          double       poly_val_common = x;
          for (unsigned int j = 0; j < degree - 3; ++j)
            poly_val_common *= Utilities::fixed_power<2>(x - jacobi_roots[j]);
          poly_val_common *= Utilities::fixed_power<4>(x - 1.);
          integral_left +=
            gauss.weight(q) * (poly_val_common * (x - guess_left));
          integral_right +=
            gauss.weight(q) * (poly_val_common * (x - guess_right));
        }

      // compute guess by secant method. Due to linearity in the root x_star,
      // this is the correct position after this single step
      return guess_right - (guess_right - guess_left) /
                             (integral_right - integral_left) * integral_right;
    }
  } // namespace



  HermiteLikeInterpolation::HermiteLikeInterpolation(const unsigned int degree,
                                                     const unsigned int index)
    : Polynomial<double>(0)
  {
    AssertIndexRange(index, degree + 1);

    this->coefficients.clear();
    this->in_lagrange_product_form = true;

    this->lagrange_support_points.resize(degree);

    if (degree == 0)
      this->lagrange_weight = 1.;
    else if (degree == 1)
      {
        if (index == 0)
          {
            this->lagrange_support_points[0] = 1.;
            this->lagrange_weight            = -1.;
          }
        else
          {
            this->lagrange_support_points[0] = 0.;
            this->lagrange_weight            = 1.;
          }
      }
    else if (degree == 2)
      {
        if (index == 0)
          {
            this->lagrange_support_points[0] = 1.;
            this->lagrange_support_points[1] = 1.;
            this->lagrange_weight            = 1.;
          }
        else if (index == 1)
          {
            this->lagrange_support_points[0] = 0;
            this->lagrange_support_points[1] = 1;
            this->lagrange_weight            = -2.;
          }
        else
          {
            this->lagrange_support_points[0] = 0.;
            this->lagrange_support_points[1] = 0.;
            this->lagrange_weight            = 1.;
          }
      }
    else if (degree == 3)
      {
        // 4 Polynomials with degree 3
        // entries (1,0) and (3,2) of the mass matrix will be equal to 0
        //
        //     | x  0  x  x |
        //     | 0  x  x  x |
        // M = | x  x  x  0 |
        //     | x  x  0  x |
        //
        if (index == 0)
          {
            this->lagrange_support_points[0] = 2. / 7.;
            this->lagrange_support_points[1] = 1.;
            this->lagrange_support_points[2] = 1.;
            this->lagrange_weight            = -3.5;
          }
        else if (index == 1)
          {
            this->lagrange_support_points[0] = 0.;
            this->lagrange_support_points[1] = 1.;
            this->lagrange_support_points[2] = 1.;

            // this magic value 5.5 is obtained when evaluating the general
            // formula below for the degree=3 case
            this->lagrange_weight = 5.5;
          }
        else if (index == 2)
          {
            this->lagrange_support_points[0] = 0.;
            this->lagrange_support_points[1] = 0.;
            this->lagrange_support_points[2] = 1.;
            this->lagrange_weight            = -5.5;
          }
        else if (index == 3)
          {
            this->lagrange_support_points[0] = 0.;
            this->lagrange_support_points[1] = 0.;
            this->lagrange_support_points[2] = 5. / 7.;
            this->lagrange_weight            = 3.5;
          }
      }
    else
      {
        // Higher order Polynomials degree>=4: the entries (1,0) and
        // (degree,degree-1) of the mass matrix will be equal to 0
        //
        //     | x  0  x  x         x  x  x |
        //     | 0  x  x  x  . . .  x  x  x |
        //     | x  x  x  0         0  x  x |
        //     | x  x  0  x         0  x  x |
        //     |     .       .         .    |
        // M = |     .         .       .    |
        //     |     .           .     .    |
        //     | x  x  0  0         x  x  x |
        //     | x  x  x  x  . . .  x  x  0 |
        //     | x  x  x  x         x  0  x |
        //
        // We find the inner points as the zeros of the Jacobi polynomials
        // with alpha = beta = 4 which is the polynomial with the kernel
        // (1-x)^4 (1+x)^4. Since polynomials (1-x)^2 (1+x)^2 are contained
        // in every interior polynomial (bubble function), their product
        // leads us to the orthogonality condition of the Jacobi(4,4)
        // polynomials.

        std::vector<double> jacobi_roots =
          jacobi_polynomial_roots<double>(degree - 3, 4, 4);
        AssertDimension(jacobi_roots.size(), degree - 3);

        // iteration from variable support point N with secant method
        // initial values

        this->lagrange_support_points.resize(degree);
        if (index == 0)
          {
            const double auxiliary_zero =
              find_support_point_x_star(jacobi_roots);
            this->lagrange_support_points[0] = auxiliary_zero;
            for (unsigned int m = 0; m < degree - 3; ++m)
              this->lagrange_support_points[m + 1] = jacobi_roots[m];
            this->lagrange_support_points[degree - 2] = 1.;
            this->lagrange_support_points[degree - 1] = 1.;

            // ensure that the polynomial evaluates to one at x=0
            this->lagrange_weight = 1. / this->value(0);
          }
        else if (index == 1)
          {
            this->lagrange_support_points[0] = 0.;
            for (unsigned int m = 0; m < degree - 3; ++m)
              this->lagrange_support_points[m + 1] = jacobi_roots[m];
            this->lagrange_support_points[degree - 2] = 1.;
            this->lagrange_support_points[degree - 1] = 1.;

            // Select the weight to make the derivative of the sum of P_0 and
            // P_1 in zero to be 0. The derivative in x=0 is simply given by
            // p~(0)/auxiliary_zero+p~'(0) + a*p~(0), where p~(x) is the
            // Lagrange polynomial in all points except the first one which is
            // the same for P_0 and P_1, and a is the weight we seek here. If
            // we solve this for a, we obtain the desired property. Since the
            // basis is nodal for all interior points, this property ensures
            // that the sum of all polynomials with weight 1 is one.
            std::vector<Point<1>> points(degree);
            double                ratio = 1.;
            for (unsigned int i = 0; i < degree; ++i)
              {
                points[i][0] = this->lagrange_support_points[i];
                if (i > 0)
                  ratio *= -this->lagrange_support_points[i];
              }
            Polynomial<double>  helper(points, 0);
            std::vector<double> value_and_grad(2);
            helper.value(0., value_and_grad);
            Assert(std::abs(value_and_grad[0]) > 1e-10,
                   ExcInternalError("There should not be a zero at x=0."));

            const double auxiliary_zero =
              find_support_point_x_star(jacobi_roots);
            this->lagrange_weight =
              (1. / auxiliary_zero - value_and_grad[1] / value_and_grad[0]) /
              ratio;
          }
        else if (index >= 2 && index < degree - 1)
          {
            this->lagrange_support_points[0] = 0.;
            this->lagrange_support_points[1] = 0.;
            for (unsigned int m = 0, c = 2; m < degree - 3; ++m)
              if (m + 2 != index)
                this->lagrange_support_points[c++] = jacobi_roots[m];
            this->lagrange_support_points[degree - 2] = 1.;
            this->lagrange_support_points[degree - 1] = 1.;

            // ensure that the polynomial evaluates to one at the respective
            // nodal point
            this->lagrange_weight = 1. / this->value(jacobi_roots[index - 2]);
          }
        else if (index == degree - 1)
          {
            this->lagrange_support_points[0] = 0.;
            this->lagrange_support_points[1] = 0.;
            for (unsigned int m = 0; m < degree - 3; ++m)
              this->lagrange_support_points[m + 2] = jacobi_roots[m];
            this->lagrange_support_points[degree - 1] = 1.;

            std::vector<Point<1>> points(degree);
            double                ratio = 1.;
            for (unsigned int i = 0; i < degree; ++i)
              {
                points[i][0] = this->lagrange_support_points[i];
                if (i < degree - 1)
                  ratio *= 1. - this->lagrange_support_points[i];
              }
            Polynomial<double>  helper(points, degree - 1);
            std::vector<double> value_and_grad(2);
            helper.value(1., value_and_grad);
            Assert(std::abs(value_and_grad[0]) > 1e-10,
                   ExcInternalError("There should not be a zero at x=1."));

            const double auxiliary_zero =
              find_support_point_x_star(jacobi_roots);
            this->lagrange_weight =
              (-1. / auxiliary_zero - value_and_grad[1] / value_and_grad[0]) /
              ratio;
          }
        else if (index == degree)
          {
            const double auxiliary_zero =
              find_support_point_x_star(jacobi_roots);
            this->lagrange_support_points[0] = 0.;
            this->lagrange_support_points[1] = 0.;
            for (unsigned int m = 0; m < degree - 3; ++m)
              this->lagrange_support_points[m + 2] = jacobi_roots[m];
            this->lagrange_support_points[degree - 1] = 1. - auxiliary_zero;

            // ensure that the polynomial evaluates to one at x=1
            this->lagrange_weight = 1. / this->value(1.);
          }
      }
  }



  std::vector<Polynomial<double>>
  HermiteLikeInterpolation::generate_complete_basis(const unsigned int degree)
  {
    std::vector<Polynomial<double>> basis(degree + 1);

    for (unsigned int i = 0; i <= degree; ++i)
      basis[i] = HermiteLikeInterpolation(degree, i);

    return basis;
  }

} // namespace Polynomials

// ------------------ explicit instantiations --------------- //

#ifndef DOXYGEN
namespace Polynomials
{
  template class Polynomial<float>;
  template class Polynomial<double>;
  template class Polynomial<long double>;

  template void
  Polynomial<float>::shift(const float offset);
  template void
  Polynomial<float>::shift(const double offset);
  template void
  Polynomial<double>::shift(const double offset);
  template void
  Polynomial<long double>::shift(const long double offset);
  template void
  Polynomial<float>::shift(const long double offset);
  template void
  Polynomial<double>::shift(const long double offset);

  template class Monomial<float>;
  template class Monomial<double>;
  template class Monomial<long double>;
} // namespace Polynomials
#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE