LinearForm.h

/*
 * ########################################################################
 * The contents of this file is free and unencumbered software released into the
 * public domain. For more information, please refer to <http://unlicense.org/>
 * ########################################################################
 */
 
#ifndef BSPLINE_INTEGRATION_LINEARFORM_H
#define BSPLINE_INTEGRATION_LINEARFORM_H
 
#include <bspline/Spline.h>
#include <bspline/operators/GenericOperators.h>
 
namespace bspline::integration {
 
template <typename O,
          std::enable_if_t<operators::is_operator_v<O>, bool> = true>
class LinearForm final {
 private:
  O _o;
 
  template <typename T, size_t size>
  static T evaluateInterval(const std::array<T, size> &a, const T &dxhalf) {
    // Use Horner's scheme to evaluate.
    constexpr int endIndex =
        static_cast<int>(size) - static_cast<int>(size % 2 == 0 ? 2 : 1);
    static_assert(endIndex >= 0);
 
    T result = a.at(endIndex) / static_cast<T>(endIndex + 1);
    const T dxhalf_squared = dxhalf * dxhalf;
 
    for (int i = endIndex - 2; i >= 0; i -= 2) {
      result = dxhalf_squared * result + a[i] / static_cast<T>(i + 1);
    }
    return static_cast<T>(2) * dxhalf * result;
  }
 
 public:
  explicit LinearForm(O o) : _o(std::move(o)){};
 
  LinearForm() : _o(O{}){};
 
  template <typename T, size_t order>
  T evaluate(const Spline<T, order> &a) const {
    const size_t nintervals = a.getSupport().numberOfIntervals();
 
    T result = static_cast<T>(0);
 
    for (size_t i = 0; i < nintervals; i++) {
      const size_t absIndex = a.getSupport().absoluteFromRelative(i);
      const T dxhalf =
          (a.getSupport()[i + 1] - a.getSupport()[i]) / static_cast<T>(2);
 
      result +=
          evaluateInterval(_o.transform(a.getCoefficients()[i],
                                        a.getSupport().getGrid(), absIndex),
                           dxhalf);
    }
    return result;
  }
 
  template <typename T, size_t order>
  T operator()(const Spline<T, order> &a) const {
    return evaluate(a);
  }
};
 
LinearForm()->LinearForm<operators::IdentityOperator>;
 
}  // namespace bspline::integration
#endif  // BSPLINE_INTEGRATION_LINEARFORM_H