...one of the most highly
regarded and expertly designed C++ library projects in the
world.

— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards

#include <boost/math/interpolators/whittaker_shannon.hpp>

namespace boost { namespace math { namespace interpolators { template <class RandomAccessContainer> class whittaker_shannon { public: using Real = RandomAccessContainer::value_type; whittaker_shannon(RandomAccessContainer&& v, Real left_endpoint, Real step_size); Real operator()(Real x) const; Real prime(Real x) const; }; }}} // namespaces

The Whittaker-Shannon interpolator takes equispaced data and interpolates between
them via a sum of sinc functions. This interpolation is stable and infinitely
smooth, but has linear complexity in the data, making it slow relative to compactly-supported
b-splines. In addition, we cannot pass an infinite amount of data into the
class, and must truncate the (perhaps) infinite sinc series to a finite number
of terms. Since the sinc function has slow *1/x* decay,
the truncation of the series can incur large error. Hence this interpolator
works best when operating on samples of compactly-supported functions. Here
is an example of interpolating a smooth "bump function":

auto bump = [](double x) { if (std::abs(x) >= 1) { return 0.0; } return std::exp(-1.0/(1.0-x*x)); }; double t0 = -1; size_t n = 2049; double h = 2.0/(n-1.0); std::vector<double> v(n); for(size_t i = 0; i < n; ++i) { double t = t0 + i*h; v[i] = bump(t); } auto ws = whittaker_shannon(std::move(v), t0, h); double y = ws(0.3);

The derivative of the interpolant can also be evaluated, but the accuracy is not as high:

double yp = ws.prime(0.3);

The call to the constructor requires 𝑶(1) operations, simply moving data into
the class. Each call to the interpolant is 𝑶(*n*), where
*n* is the number of points to interpolate.