python計算對角線有理函數插值的方法

''' p = rational(xData,yData,x)
 Evaluates the diagonal rational function interpolant p(x)
 that passes through he data points
''' 
from numpy import zeros
def rational(xData,yData,x):
 m = len(xData)
 r = yData.copy()
 rOld = zeros(m)
 for k in range(m-1):
 for i in range(m-k-1):
 if abs(x - xData[i+k+1]) < 1.0e-9:
 return yData[i+k+1]
 else:
 c1 = r[i+1] - r[i]
 c2 = r[i+1] - rOld[i+1]
 c3 = (x - xData[i])/(x - xData[i+k+1])
 r[i] = r[i+1] + c1/(c3*(1.0 - c1/c2) - 1.0)
 rOld[i+1] = r[i+1]
 return r[0]