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$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.