Aitken-Neville interpolation, version 2

If Cut/Copy and Paste fails, then click here for download.

void aitkenneville2 (unsigned int n, unsigned int m, double *xi, double *yi, double *xo, double *yo, double *w)
{
  /* xi[n], yi[n], xo[m], yo[m], w[((n + 1) * n >> 1) + 1] */
  unsigned int i, ii, j, k;
  double yok;
  for (i = 0; i < n; i++)
    w[i] = yi[i];
  for (i = 1; i < n; i++) {
    ii = (n + n - i) * (i - 1) >> 1;
    for (j = i; j < n; j++)
      w[ii + j + n - i] = (w[ii + j] - w[ii + j - 1]) / (xi[j] - xi[j - i]);
    }
  for (i = 0; i < n; i++)
    w[n - (i + 1)] = w[(n + n + 1 - i) * i >> 1];
  for (i = 0; i < n; i++)
    w[n + n - (i + 1)] = xi[i];
  if (n > 0)
    for (k = 0; k < m; k++) {
      yok = w[0];
      for (i = 1; i < n; i++)
        yok = yok * (xo[k] - w[n + i]) + w[i];
      yo[k] = yok;
      }
  return;
  }