Linear Interpolation for sorted arrays 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 &lt;algorithm&gt;
#include &lt;cmath&gt;
#include &lt;type_traits&gt;
#include &lt;vector&gt;

#if (__cplusplus &lt; 202002L)
template &lt;typename T,
          typename = std::enable_if_t&lt;std::is_floating_point&lt;T&gt;::value&gt;&gt;
T lerp(T a, T b, T t) noexcept {
  return a + t * (b - a);
}
#else
using std::lerp;
#endif

template &lt;typename T,
          typename = std::enable_if_t&lt;std::is_floating_point&lt;T&gt;::value&gt;&gt;
class ListLerp {
 public:
  explicit ListLerp(const std::vector&lt;T&gt;&amp; xp, const std::vector&lt;T&gt;&amp; yp)
      : xp(xp), yp(yp) {
    b = this-&gt;xp.begin() + 1;
  }

  [[nodiscard]] T interp(const T x) noexcept {
    // No extrapolation
    if (x &lt;= xp.front()) return yp.front();
    if (x &gt;= xp.back()) return yp.back();

    // b = xp.begin() + 1;

    // find b position
    while (x &gt; *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&lt;T&gt;::const_iterator b;
  const std::vector&lt;T&gt;&amp; xp;
  const std::vector&lt;T&gt;&amp; yp;
};

template &lt;typename T,
          typename = std::enable_if_t&lt;std::is_floating_point&lt;T&gt;::value&gt;&gt;
std::vector&lt;T&gt; interp(const std::vector&lt;T&gt;&amp; xp, const std::vector&lt;T&gt;&amp; yp,
                      const std::vector&lt;T&gt;&amp; x) {
  std::vector&lt;T&gt; out(x.size());

  auto lLerp = ListLerp(xp, yp);

  std::transform(x.begin(), x.end(), out.begin(),
                 [&amp;lLerp](const auto&amp; xi) { return lLerp.interp(xi); });

  return out;
}


///////////////
int main() {
  std::vector&lt;double&gt; xp = {1, 2, 3};
  std::vector&lt;double&gt; fp = {3, 2, 0};

  for (auto &amp;i : interp(xp, fp, {1, 1.5, 2.72})){
    std::cout &lt;&lt; " " &lt;&lt; i; 
  }
  std::cout &lt;&lt; "\n";

}
template &lt;typename T,
          typename = std::enable_if_t&lt;std::is_floating_point&lt;T&gt;::value&gt;&gt;
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 &lt;typename T,
          typename = std::enable_if_t&lt;std::is_floating_point&lt;T&gt;::value&gt;&gt;
std::vector&lt;T&gt; interp(const std::vector&lt;T&gt;&amp; xp, const std::vector&lt;T&gt;&amp; yp,
                      const std::vector&lt;T&gt;&amp; 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&amp; 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&amp; a, vec const&amp; b, vec const&amp; 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&lt;T&gt;&amp; 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&lt;double&gt; 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&amp; 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 &gt; *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]"