Line data    Source code

       1           1 : // Distributed under the MIT License.
       2             : // See LICENSE.txt for details.
       3             : 
       4             : /// \file
       5             : /// Defines functions lagrange_polynomial
       6             : 
       7             : #pragma once
       8             : 
       9             : #include <iterator>
      10             : 
      11             : #include "Utilities/ErrorHandling/Assert.hpp"
      12             : 
      13             : /// Evaluate the jth Lagrange interpolating polynomial with the given
      14             : /// control points where j is the index of index_point.
      15             : template <typename Iterator>
      16           1 : double lagrange_polynomial(const Iterator& index_point,
      17             :                            double x,
      18             :                            const Iterator& control_points_begin,
      19             :                            const Iterator& control_points_end) {
      20             :   ASSERT(control_points_begin != control_points_end, "No control points");
      21             : 
      22             :   const double x_j = *index_point;
      23             :   double result = 1.;
      24             :   for (auto m_it = control_points_begin; m_it != control_points_end; ++m_it) {
      25             :     if (m_it == index_point) { continue; }
      26             :     const double x_m = *m_it;
      27             :     result *= (x_m - x) / (x_m - x_j);
      28             :   }
      29             :   return result;
      30             : }
      31             : 
      32             : /// Evaluate the jth (zero-indexed) Lagrange interpolating polynomial
      33             : /// with the given control points.
      34             : template <typename Iterator>
      35           1 : double lagrange_polynomial(size_t j, double x,
      36             :                            const Iterator& control_points_begin,
      37             :                            const Iterator& control_points_end) {
      38             :   const auto j_diff =
      39             :       static_cast<typename std::iterator_traits<Iterator>::difference_type>(j);
      40             :   ASSERT(j_diff < std::distance(control_points_begin, control_points_end),
      41             :          "Polynomial number out of range " << j << " > "
      42             :          << (std::distance(control_points_begin, control_points_end) - 1));
      43             :   return lagrange_polynomial(std::next(control_points_begin, j_diff),
      44             :                              x, control_points_begin, control_points_end);
      45             : }