使用截断的 Whitaker-Shannon 级数（基数级数） 我们可以通过重复调用 sinc 例程来天真地计算总和，例如以下代码做：

$$f(t) = \sum_{j = 0}^{n-1} y_{j} \frac{\sin\left(\pi( \frac{t-t_0}{h} -j)\right)}{\pi\left(\frac{t-t_0}{h}-j\right)}$$

   Real naive_sum(Real t) const {
        using boost::math::constants::pi;
        Real f = 0;
        for (size_t i = 0; i < m_y.size(); ++i) {
            Real arg = pi<Real>()*( (t-m_t0)/m_h - i);
            f += m_y[i]*boost::math::sinc_pi(arg);
        }
        return y;
    }

然而，重复调用 sinc 代价高昂，因此我们可以使用恒等式将其写为 $\sin(\theta - j\pi) = (-1)^{j}\sin(\theta)$

$$f(t) = \frac{\sin(\pi(t-t_0)/h)}{\pi} \sum_{j=0}^{n-1} (-1)^{j}\frac{y_{j}}{(t-t_0)/h -j}$$

可以按如下方式实现：

    Real operator()(Real t) const {
        using boost::math::constants::pi;
        using std::sin;
        Real y = 0;
        Real x = (t - m_t0)/m_h;

        for (size_t i = 0; i < m_y.size(); ++i)
        {
            Real denom = (x - i);
            if (denom == 0) {
                return m_y[i];
            }
            if (i & 1) {
                y -= m_y[i]/denom;
            }
            else {
                y += m_y[i]/denom;
            }
        }
        return y*sin(pi<Real>()*x)/pi<Real>();
    }

但是，我观察到使用快速方法比使用慢速方法的准确性大大降低。能否在不大幅降低准确度的情况下保持快速方法的速度？

工作代码，对于那些关心的人：

#ifndef BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_HPP
#define BOOST_MATH_INTERPOLATORS_WHITAKKER_SHANNON_HPP
#include <boost/math/special_functions/sinc.hpp>
#include <boost/math/constants/constants.hpp>

namespace boost::math::interpolators {

template<class RandomAccessContainer>
class whittaker_shannon {
public:

    using Real = typename RandomAccessContainer::value_type;
    whittaker_shannon(RandomAccessContainer&& y, Real t0, Real h) : m_y{std::move(y)}, m_t0{t0}, m_h{h}
    {
    }

    Real operator()(Real t) const {
        using boost::math::constants::pi;
        using std::sin;
        Real y = 0;
        Real x = (t - m_t0)/m_h;

        for (size_t i = 0; i < m_y.size(); ++i)
        {
            Real denom = (x - i);
            if (denom == 0) {
                return m_y[i];
            }
            if (i & 1) {
                y -= m_y[i]/denom;
            }
            else {
                y += m_y[i]/denom;
            }
        }
        return y*sin(pi<Real>()*x)/pi<Real>();
    }


    Real naive_sum(Real t) const {
        using boost::math::constants::pi;
        Real y = 0;
        Real s = pi<Real>()*(t-m_t0)/m_h;
        for (size_t i = 0; i < m_y.size(); ++i) {
            Real arg = pi<Real>()*( (t-m_t0)/m_h - i);
            y += m_y[i]*sinc_pi(arg);
        }
        return y;
    }


    Real operator[](size_t i) const {
        return m_y[i];
    }


private:
    RandomAccessContainer m_y;
    Real m_t0;
    Real m_h;
};
}
#endif

这是一个重现该现象的测试：

template<class Real>
void test_bump()
{
    using std::exp;
    using std::abs;
    auto bump = [](Real x) { if (abs(x) >= 1) { return Real(0); } return exp(-Real(1)/(Real(1)-x*x)); };

    Real t0 = -1;
    size_t n = 2049;
    Real h = Real(2)/Real(n-1);

    std::vector<Real> v(n);
    for(size_t i = 0; i < n; ++i) {
        Real t = t0 + i*h;
        v[i] = bump(t);
    }


    auto ws = whittaker_shannon(std::move(v), t0, h);

    std::mt19937 gen(323723);
    std::uniform_real_distribution<long double> dis(-0.95, 0.95);

    size_t i = 0;
    while (i++ < 1000)
    {
        Real t = static_cast<Real>(dis(gen));
        Real expected = bump(t);
        if(!CHECK_MOLLIFIED_CLOSE(expected, ws(t), 50*std::numeric_limits<Real>::epsilon())) {
            std::cerr << "  Problem occured at abscissa " << t << "\n";
        }
    }
}