Linear Interpolation for sorted arrays

Question

I have made a linear interpolation functions as a side project of mine. It assumes everything is sorted before hand - x and f(x) are the same length.

I would like to ask for:

• general recommendations
• if I can improve the design (generalize it to all iterables like std::array etc...)
• and for further performance optimizations.

#include <algorithm>
#include <cmath>
#include <type_traits>
#include <vector>
#if (__cplusplus < 202002L)
template <typename T,
 typename = std::enable_if_t<std::is_floating_point<T>::value>>
T lerp(T a, T b, T t) noexcept {
 return a + t * (b - a);
}
#else
using std::lerp;
#endif
template <typename T,
 typename = std::enable_if_t<std::is_floating_point<T>::value>>
class ListLerp {
 public:
 explicit ListLerp(const std::vector<T>& xp, const std::vector<T>& yp)
 : xp(xp), yp(yp) {
 b = this->xp.begin() + 1;
 }
 [[nodiscard]] T interp(const T x) noexcept {
 // No extrapolation
 if (x <= xp.front()) return yp.front();
 if (x >= xp.back()) return yp.back();
 // b = xp.begin() + 1;
 // find b position
 while (x > *b) b++;
 // Calculate a, a_y, b_y from b
 const auto a = b - 1;
 const auto a_y = yp.begin() + (a - xp.begin());
 const auto b_y = yp.begin() + (b - xp.begin());
 const auto t = (x - *a) / (*b - *a);
 return lerp(*a_y, *b_y, t);
 };
 private:
 typename std::vector<T>::const_iterator b;
 const std::vector<T>& xp;
 const std::vector<T>& yp;
};
template <typename T,
 typename = std::enable_if_t<std::is_floating_point<T>::value>>
std::vector<T> interp(const std::vector<T>& xp, const std::vector<T>& yp,
 const std::vector<T>& x) {
 std::vector<T> out(x.size());
 auto lLerp = ListLerp(xp, yp);
 std::transform(x.begin(), x.end(), out.begin(),
 [&lLerp](const auto& xi) { return lLerp.interp(xi); });
 return out;
}
///////////////
int main() {
 std::vector<double> xp = {1, 2, 3};
 std::vector<double> fp = {3, 2, 0};
 for (auto &i : interp(xp, fp, {1, 1.5, 2.72})){
 std::cout << " " << i; 
 }
 std::cout << "\n";
}

Answer 1

template <typename T,
 typename = std::enable_if_t<std::is_floating_point<T>::value>>
T lerp(T a, T b, T t) noexcept;

Looks fine. I'd prefer to see static_assert than enable_if since there's no other overloads (and if there are I'd rather fine out with a compilation error).

T lerp(T a, T b, T t) noexcept {
 return a + t * (b - a);
}

How about std::fma?

Note that GCC and Clang both do this differently... They do b * t + a * (1-t) and have special cases for t at 0 and near 1. Maybe it's better to just copy their code exactly so the compiler version doesn't change the output?

---

template <typename T,
 typename = std::enable_if_t<std::is_floating_point<T>::value>>
std::vector<T> interp(const std::vector<T>& xp, const std::vector<T>& yp,
 const std::vector<T>& x);

There's nothing in this function that directly needs T to be a floating point type. I think it's better to remove the enable_if and let the compiler show a traceback with the "real" reason T needs to be an FP.

This function takes vectors and returns a new vector. That is a fine and very common pattern, but it's not optimally efficient. If this is a hot function, you could return an iterator/some sort of range that computes the next result when it is incremented. Boost example of this kind of thing: https://www.boost.org/doc/libs/1_66_0/libs/iterator/doc/html/iterator/specialized/transform.html.

If you want to generalize this to other iterables, then don't take std::vector... instead take a templated T const& and same for ListLerp.

---

[[nodiscard]] T interp(const T x) noexcept;

Suppose I want to discard the result? Maybe I only want to store every nth item? interp is not a const member function (one might want the mutation to happen...) and T is not an "undiscardable" type. I don't see why this needs to be here.

---

I don't think you've separated your interests very well between ListLerp::interp and the free function interp.

I think there are two good options:

• make ListLerp::interp a const member function and get rid of ListLerp::b. Then the free function interp just does a "pure" transformation of a single list.

• iterate all three lists in a single function (my preference).

E.g.

auto interp(vec const& a, vec const& b, vec const& t) {
 for idx = each index in a, b {
 ... Possibly emit a value corresponding to t ...
 }
}

This is very basic pseudo code, but the key idea is you iterate everything in one go rather than doing some iteration in a free function and some iteration in a class that holds potentially dangerous references.

---

const std::vector<T>& xp;

This (and the lines like it) is moderately dangerous since you need to ensure the vector that is passed in outlives the class.

E.g. you cannot do this:

auto fn() {
 vector<double> v;
 ListLerp ll(v); // take reference to local variable v
 return ll; // v goes out of scope; ll has dangling reference to dead v
}

Often times the best option is to take the const& and be careful about how you use it, but in this case I think it's easy enough to compute the full result in a function, not expose a class that can hold a reference to a ctor argument, and make it impossible to end up with a dangling reference.

---

It assumes everything is sorted before hand - x and f(x) are the same length.

How about asserting this?

---

// find b position
while (x > *b) b++;

If your input is sorted, this could be a binary search. May or may not be better.

---

I think your calculation at the end would be simpler if you use indices rather than iterators. Many of the lines are basically "find the iterator in A that has the same index as this iterator in B" ... easier to say "A[i]"

• Stavros Avramidis: Thanks for your answer. I will play around with the suggestions. I have 2 questions. One about your note in the separation between ListLerp::interp - if you could give me an example. The reason ListLerp exists as a class in general is not to do a linear search in every item. Also if you could explain the co st reference issue a bit.
• sudo rm -rf slash: I made a few edits in the answer. Lmk what you think.
• Stavros Avramidis: I see what you mean now, thx