3
$\begingroup$

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?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.