Simple Lagrange interpolation

Question

I wrote a Lagrange polynomial interpolation class. It works fine. At the beginning some definitions: $$ \begin{align} &nodes:\qquad x_0, \dots, x_n \\ &values:\qquad f(x_0), \dots, f(x_n) \\ &weights:\qquad w_0, \dots, w_n\qquad\text{where}~~~ w_j = \prod_{k=0, k \neq j}^{n}(x_j-x_k)^{-1} \\ &coeff:\qquad f(x_0)w_0, \dots, f(x_n)w_n. \end{align} $$

Let's denote $$c_i = f(x_i)w_i.$$ Value of approximation at some point is equal to $$ \begin{align} p(x) = l(x) \sum_{i=0}^{n} c_i(x-x_i)^{-1}, \qquad \text{where}~~l(x) = \prod_{i=0}^{n}(x-x_i). \end{align} $$

Code (firstly lagrange_interpolation.h):

#include <vector>
#include <valarray>
#include <algorithm>
#include <iostream>
#include <stdexcept>
#include <numeric>
template<typename T, typename V>
int binary_search_find_index(const std::vector<T> &v, V data) {
 std::size_t index = 0;
 auto it = std::lower_bound(v.begin(), v.end(), data);
 if (it == v.end() || *it != data) {
 return -1;
 } else {
 index = std::distance(v.begin(), it);
 } 
 return index;
}
class LInter {
typedef std::vector<double> vector;
public:
 LInter(vector nodes, vector values) {
 if (nodes.empty() || values.empty()) {
 throw std::invalid_argument("one of the vector is empty");
 }
 
 std::vector<int> indices(nodes.size());
 std::iota(indices.begin(), indices.end(), 0);
 std::sort(indices.begin(), indices.end(),[&](int a, int b) -> bool {
 if (nodes[a] == nodes[b]) {
 throw std::invalid_argument("nodes are not unique");
 }
 return nodes[a] < nodes[b]; }
 );
 
 if (nodes.size() != values.size()) {
 throw std::invalid_argument("sizes of vectors differ");
 }
 
 for (int i = 0; i < nodes.size(); i++) {
 std::swap(nodes[i], nodes[indices[i]]);
 std::swap(values[i], values[indices[i]]);
 }
 nodes_ = nodes;
 values_ = values;
 poly_deg_ = nodes_.size() - 1;
 GenerateWeights();
 GenerateCoefficients();
 }
 template<typename T>
 std::enable_if_t<std::is_arithmetic_v<T>, double> operator()(T);
 
 template<typename T>
 std::enable_if_t<!std::is_arithmetic_v<T>, std::valarray<double>> 
 operator()(T const &);
 
private:
 vector nodes_;
 vector values_;
 vector weights_;
 vector coeff_;
 unsigned poly_deg_;
 void GenerateWeights(void);
 void GenerateCoefficients(void);
};

and lagrange_interpolation.cc:

#include "lagrange_interpolation.h"
#include <type_traits>
void LInter::GenerateWeights(void) {
 int size = poly_deg_ + 1;
 weights_.insert(weights_.begin(), size, 1);
 for (int i = 0; i < size; i++){
 for (int j = 0; j < size; j++) {
 if (i != j) {
 weights_[i] *= (nodes_[i] - nodes_[j]);
 }
 }
 weights_[i] = 1/weights_[i];
 }
}
void LInter::GenerateCoefficients(void) {
 int size = poly_deg_ + 1;
 coeff_.insert(coeff_.begin(), size, 1);
 std::transform(values_.begin(), values_.end(), weights_.begin(), coeff_.begin(), [](double v, double w) { return v*w;});
}
template<typename T>
std::enable_if_t<std::is_arithmetic_v<T>, double> LInter::operator()(T x) {
 auto index = binary_search_find_index(nodes_, x);
 if (index != -1) { 
 return values_[index];
 }
 double polynomial = 1;
 for (auto node : nodes_) {
 polynomial *= (x - node);
 }
 double result = std::inner_product(coeff_.begin(), coeff_.end(), nodes_.begin(), 0.0, std::plus<>(),[&](double &a, double &b) { return a/(x - b); });
 result *= polynomial;
 return result;
}
template<typename T>
std::enable_if_t<!std::is_arithmetic_v<T>, std::valarray<double>> 
 LInter::operator()(T const &input) {
 std::valarray<double> result(input.size());
 std::transform(std::begin(input), std::end(input), std::begin(result), 
 [this](double x) { return operator()(x); } );
 return result;
}
```

• Interpolated: Could someone edit my post and change coeff to \text{coeff} (and other) in align enviroment? Unfortunately I can't do this....

Answer 1

binary_search_find_index

I'm not sure why you are using different template types for the vector and data.

It looks like you are trying to find the index of the closest element in the vector but then you discard it if it's not exactly equal to data

if (it == v.end() || *it != data) {
 return -1;
} else ...

Perhaps this is closer to what you initially intended to do:

if (it == v.end() || std::abs(*it - data) > someVerySmallNumber {
 return -1;
} else...

Since I cannot really guess what your intentions were, I will just assume that you are trying to get the index of an element in an array.

First I would change the return type to std::size_t which is the correct type for array size and indexes. You can also use std::optional<std::size_t> and std::find like the following:

template<typename T>
std::optional<std::size_t> binary_search_find_index(const std::vector<T>& v, T data) {
 auto index = std::find(v.begin(), v.end(), data);
 if (index == v.end())
 return std::nullopt;
 return index - v.begin();
}

Then you can use binary_search_find_index like so:

void test() {
 std::vector<int> v = { 1, 2, 3, 4, 5 };
 auto index = binary_search_find_index(v, 3);
 if (index == std::nullopt)
 std::cout << "item is not in array.\n";
 else
 std::cout << "index = " << *index << '\n';
}

Since you are only using binary_search_find_index in one location, I think it might be better to actually remove it from the program and rather change ::operator()(T x) to the following:

template<class T>
std::enable_if_t<std::is_arithmetic_v<T>, double> LInter<T>::operator()(T x) {
 auto node = std::find(nodes_.begin(), nodes_.end(), x);
 if (node != nodes_.end()) {
 return *node;
 }
 // ...
 return result;
}

LInter

I would have rather templated the class instead of using typedef std::vector<double> vector;.

Consider changing to the following:

template<class T>
class LInter {
public:
 LInter(std::vector<T> nodes, std::vector<T> values) {
 if (nodes.empty() || values.empty()) {
 throw std::invalid_argument("one of the vector is empty");
 }
 }
 // ...
}

Also note that poly_deg_ should be std::size_t poly_deg_ instead of simply unsigned poly_deg_

Test cases

Reviewing the code would have been easier and more productive if you included some test cases.