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              : }