Sum over ghosts or use a diffrent formula for periodic cardinal sinus.

Hi, working for the variation of the Quintic kernel #1265 along the line of Gary & Daniel 2014 paper, I mention that the code below in Interpolant.cpp

    double Interpolant::xvalWrapped(double x, int N) const
    {
        // sum over all arguments x+jN that are within range.
        // Start by finding x+jN closest to zero
        double xdown = x - N*std::floor(x/N + 0.5);
        xassert(std::abs(xdown) <= N);
        if (xrange() <= N) {
            // This is the usual case.
            return xval(xdown);
        } else {
            double xup = xdown+N;
            double sum = 0.;
            while (std::abs(xdown) <= xrange()) {
                sum += xval(xdown);
                xdown -= N;
            }
            while (xup <= xrange()) {
                sum += xval(xup);
                xup += N;
            }
            return sum;
        }
    }

implements the following equation of the 2014 paper

But, discussing with Gary, I've mebtioned that in fact the true kernel (aka periodic cardinal sinus)

K^{\mathrm{ideal}}_u(\nu) =  \frac{1}{N}\sum_{j=-N/2}^{N/2-1} e^{-2\pi i j \nu} = e^{i\pi \nu} \frac{\mathrm{sinc}(N \nu)}{\mathrm{sinc}(\nu)}$$

can be directly approximated by any "sinc" approximated kernel (noted ker) with the following generic python code
(of course this is not an optimized code)

def k_p(x,N,ker,xmax):
  x = np.abs(x)
  th=xmax/N
  return np.piecewise(x, [x<th,x>=th],[lambda x: ker(N*x)/ker(x), lambda x:0])

def rkWrapped(u,N,ker,xmax):
  return np.sign(0.5-np.floor(u+0.5)%2)*k_p(u-np.floor(u+0.5),N,ker,xmax)

def kernelWrapped(u,N,ker,umax):
  return np.exp(1j * np.pi * u) * kWrapped(u,N,ker,xmax)

and then one perform a single convolution (Eq. of 11 in the paper) with the wrapped kernel above

\tilde{F}(u) =  K_x(u) \times  \sum_{k=-N/2}^{N/2-1} \hat{a}_k\times   \tilde{K}_u(u-k/N)

One can reproduce the figure of the paper using such convolution.

Gary says that : "it would be interesting to see numerically how accurate this is. I suppose it’s a speedup of a few since you don’t have to sum over the aliased frequencies any more."

Start in Interpolant.cpp and read xvalWrapped alongside the linked 2014 paper and the proposed Python kernel functions. Numerically compare the existing aliased-sum approach with the wrapped-kernel convolution, measuring accuracy and speed; the work is done when the preferred formulation and its tradeoffs are established.