aitkens_interpolation.py

def aitken_interpolation(x, y):
  if len(x) != len(y):
    raise ValueError()
  
  # Constants
  x_poly = np.polynomial.Polynomial([0, 1])
  
  # Initialize algorithm
  n = len(x)
  p = [np.polynomial.Polynomial([y[0]])]
  
  # Diagonal values
  p_diag = [0] * len(x) 
  p_diag[0] = y[0]

  for k in range(1, n):
    p = [np.polynomial.Polynomial([0.0])] * (k + 1)
    p[0] = np.polynomial.Polynomial([y[k]])
    
    xk_poly = x[k] - x_poly
    for i in range(1, k + 1):
      xi_poly = x[i-1] - x_poly
      
      p[i] = xk_poly * p_diag[i-1] - xi_poly * p[i-1]
      p[i] /= x[k] - x[i-1]
      
    p_diag[k] = p[-1]
          
  return p[-1]