Stencil.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// SPDX-FileCopyrightText: The Eigen Authors
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_NUMERICALDIFF_STENCIL_H
#define EIGEN_NUMERICALDIFF_STENCIL_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {

template <int Derivative, typename Scalar, int Size>
class Stencil {
  EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
  static_assert(Size >= 0, "Stencil requires a fixed-size point set");
  static_assert(Derivative >= 0 && Derivative < Size, "Stencil requires at least `Derivative + 1` points");
 
 public:
  static constexpr int DerivativeOrder = Derivative;
  static constexpr int PointCount = Size;

  using ArrayType = Array<Scalar, Size, 1>;

  template <typename Derived>
  explicit constexpr Stencil(const DenseBase<Derived>& points) {
    static_assert(Derived::IsVectorAtCompileTime, "Stencil requires a 1D dense expression");
    static_assert(Derived::SizeAtCompileTime == Size, "Stencil requires exactly `Size` points");
    for (Index i = 0; i < Size; ++i) m_points[i] = points.coeff(i);
    computeWeights();
  }

  explicit constexpr Stencil(const Scalar (&points)[Size]) {
    for (Index i = 0; i < Size; ++i) m_points[i] = points[i];
    computeWeights();
  }

  ArrayType weights() const { return Map<const ArrayType>(m_weights); }

  ArrayType points() const { return Map<const ArrayType>(m_points); }

  constexpr const Scalar& weight(Index i) const { return m_weights[i]; }

  constexpr const Scalar& point(Index i) const { return m_points[i]; }
 
 private:
  // Fornberg's recursion (1988), in the in-place algorithmic form given by Fornberg (1998).
  // `delta[nu][m]` is the weight of `points[nu]` in the order-`m` formula built from `points[0..n]`
  // so far; each column `nu` only ever reads/writes its own entries as `n` grows, so a single working
  // table can be updated in place.
  constexpr void computeWeights() {
    Scalar delta[Size][Derivative + 1]{};
    delta[0][0] = Scalar(1);
 
    for (Index n = 1; n < Size; ++n) {
      const Index mn = numext::mini(n, Index(Derivative));
      const Scalar c4 = m_points[n];
      // c1/c2 = product_{j<n-1} ((x_{n-1}-x_j)/(x_n-x_j)) / (x_n-x_{n-1}).
      // Form the ratio directly to avoid overflow/underflow of the individual products.
      Scalar ratio(1);
      for (Index nu = 0; nu < n; ++nu) {
        const Scalar c3 = m_points[n] - m_points[nu];
        eigen_assert(!(c3 == Scalar(0)) && "Stencil points must be pairwise distinct");
        if (nu == n - 1) {
          ratio /= c3;
          delta[n][0] = -ratio * (m_points[n - 1] * delta[n - 1][0]);
          for (Index m = mn; m >= 1; --m)
            delta[n][m] = ratio * (Scalar(m) * delta[n - 1][m - 1] - m_points[n - 1] * delta[n - 1][m]);
        } else {
          ratio *= (m_points[n - 1] - m_points[nu]) / c3;
        }
        for (Index m = mn; m >= 1; --m) delta[nu][m] = (c4 * delta[nu][m] - Scalar(m) * delta[nu][m - 1]) / c3;
        delta[nu][0] = c4 * delta[nu][0] / c3;
      }
    }
 
    for (Index nu = 0; nu < Size; ++nu) m_weights[nu] = delta[nu][Derivative];
  }
 
  // NOTE: C arrays instead of `std::array` (through `Eigen::array`) is intentional since
  // non-`const` `operator[]` is only `constexpr` in C++17, not C++14.
  Scalar m_points[Size]{};
  Scalar m_weights[Size]{};
};
 
// Redundant out-of-class definitions are required pre-C++17 but deprecated since.
#if EIGEN_COMP_CXXVER < 17
template <int Derivative, typename Scalar, int Size>
constexpr int Stencil<Derivative, Scalar, Size>::DerivativeOrder;
template <int Derivative, typename Scalar, int Size>
constexpr int Stencil<Derivative, Scalar, Size>::PointCount;
#endif

template <int Derivative, typename Derived>
constexpr Stencil<Derivative, typename Derived::Scalar, Derived::SizeAtCompileTime> makeStencil(
    const DenseBase<Derived>& points) {
  return Stencil<Derivative, typename Derived::Scalar, Derived::SizeAtCompileTime>(points);
}

template <int Derivative, typename Scalar, int Size>
constexpr Stencil<Derivative, Scalar, Size> makeStencil(const Scalar (&points)[Size]) {
  return Stencil<Derivative, Scalar, Size>(points);
}
 
}  // namespace Eigen
 
#endif