Last time we looked at how to improve the accuracy of Whittaker-Shannon interpolation, where user njuffa demonstrated that judicious use of sin_pi
could greatly improve the accuracy of computing
$$
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)}
$$
Now I am struggling with accurate evaluation of
$$
f'(t) = \frac{\pi}{h} \sum_{j = 0}^{n-1} y_{j} \mathrm{sinc}'\left( \pi \frac{t-t_0}{h} -\pi j\right)
$$
My code for evaluating this derivative is here and for brevity the relevant portion is presented here:
Real prime(Real t) const {
using boost::math::constants::pi;
using std::isfinite;
using std::floor;
Real x = (t - m_t0)/m_h;
// The integer branch is not a problem:
if (ceil(x) == x) {
Real s = 0;
long j = static_cast<long>(x);
long n = m_y.size();
for (long i = 0; i < n; ++i)
{
if (j - i != 0)
{
s += m_y[i]/(j-i);
}
// else derivative of sinc at zero is zero.
}
if (j & 1) {
s /= -m_h;
} else {
s /= m_h;
}
return s;
}
Real z = x;
auto it = m_y.begin();
Real cospix = boost::math::cos_pi(x);
Real sinpix = boost::math::sin_pi(x);
Real s = 0;
auto end = m_y.end();
while(it != end)
{
s += (*it++)*(pi<Real>()*z*cospix - sinpix)/(z*z);
z -= 1;
}
return s/(pi<Real>()*m_h);
}
I'm looking to achieve errors of roughly 10 ULPs, but this is giving errors of roughly ~1000 ULPs, which is a sizable fraction of the bit budget. I believe the source of the error comes from the (pi<Real>()*z*cospix - sinpix)/(z*z)
term, but I haven't been able to rearrange it to get better accuracy.
Is there a way to get this more accurate?