function [ t, s, errbd, err, tcpu, scpu ] = neqfft( n, c, m, x, k, q )
% Evaluate nonequidistant FFT of n coeffs c at m arbitrary points x:
% t(x_j) = sum_{k=1}^n c(k) e^{i(k-1)x_j} for j = 1, ..., m.
% First FFT onto an q-oversampled mesh, then use centered polynomial interpolation
% of degree 2k-1. Error bounded by (pi/q)^{2k}.
% Alternatively could use Aitken or Neville formula if m is small, 
% Alternatively could evaluate derivatives t'(x), t''(x), ... and
% use Hermite interpolation to decrease stencil size and thereby slightly decrease error.
%
tic;
errbd = ( pi / q )^(2*k) * sum( abs( c ) );
nq = q * n;
hq = 2*pi/nq;
nk = nq + 2*k+1;
d = zeros( 2*k, nk );
% Take standard FFT (of automatically zero-padded coeffs c) to get nq mesh values d.
% d( j, 1 ) = sum_{k=1}^n c(k) exp( -i 2 pi (j-1) (k-1) / nq ) for j = 1, ..., nq
d( 1, 1:nq ) = fft( c, nq ) ;
%
% Add 2k+1 periodic values at right end---shift x if near left end.
d( 1, 1:nk ) = [ d( 1, 1:nq ) d( 1, 1:2*k+1 ) ];
% Compute divided differences of d --- this is efficient if m is large.
for j = 1 : 2*k-1
   % Skip dividing by hq as we won't be including it in our Newton formula.
   d( j+1, 1:nk-j ) = ( d( j, 2:nk-j+1 ) - d( j, 1:nk-j ) ) / j;
end
%
% Locate each x_j in mesh (shift at left end) and interpolate by Newton formula.
t = zeros( 1, m );
for j = 1:m
   % Index of nearest grid point x_k = (kc-1)*hq left of x(j)=(kc-1+tc)*hq.
   kc = 1 + floor( x(j) / hq );
   % Index of left endpoint of stencil.
   kl = kc - k + 1;
   % Integer multiple plus fraction of mesh spacing from kl to x(j).
   tl = ( x(j) - ( kl-1 ) * hq ) / hq;
   % Wrap around if kl runs off left end of periodic array.
   if kl < 1 
      kl = kl + nq;
   end
   % Horner's rule for Newton interpolation from kl to t.
   t(j) = d( 2*k, kl );
   for s = 2*k-2:-1:0
      t(j) = d( s+1, kl ) + ( tl-s ) * t(j);
   end
end
tcpu = toc;
clear d
% Check against direct evaluation
tic;
s = zeros( 1, m );
for j = 1:m
   sum = 0.0;
   for k = 1:n
      sum = sum + c(k) * exp( -i * (k-1) * x(j) );
   end
   s(j) = sum;
end
err = norm( s - t );
scpu = toc;