Whittaker-Shannon interpolation: Accuracy dies with speedup; can it be fixed?

Question

With a truncated Whitaker-Shannon series (cardinal series) $$ 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)} $$ we can naively evaluate the sum by repeated calls to sinc routines, such as the following code does:

   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;
    }

However, repeated calls to sinc are expensive, so we can use the identity $\sin(\theta - j\pi) = (-1)^{j}\sin(\theta)$ to write this as $$ 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} $$

which can be implemented as follows:

    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>();
    }

However, I have observed a vast decrease in accuracy using the fast method over the slow method. Can the speed of the fast method be preserved without a massive decrease in accuracy?

Working code, for those that care:

#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

Here's a test that reproduces the phenomenon:

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";
        }
    }
}

Answer 1

I was able to reproduce the behavior reported in the question, and traced the observed inaccuracies to the following line:

return y*sin(pi<Real>()*x)/pi<Real>();

The explicit multiplication with a floating-point approximation of π introduces a small error into the argument to sin, which comprises the representational error in the constant and the rounding error added by the multiply and is on the order of 1 ulp. This small error in the argument represents a phase error in sin which grows with the magnitude of x. In this case |x| is on the order of 1000, resulting in quite a bit of error magnification.

Because this is a relatively frequent scenario, various platforms offer a sinpi function that computes sin (π x) with the multiplication by π happening inside sinpi after argument reduction, minimizing the phase error. In the specific case of Boost the desired function is boost::math::sin_pi. Often the sinpi function is also more efficient than the regular sin function, since its internal argument reduction is simpler.

Replacing the original source line with the line below should fix the observed accuracy issue completely:

return y*boost::math::sin_pi(x)/pi<Real>();