/**
拉格朗日插值計算
拉格朗日插值公式：
                n     n
  Pn(x(i))= ∑〔  ∏（ x-x(j)）/（x(k)-x(j)） 〕y(k)
               k=0   j=0,
                      j≠k
 
  屬性：插值計算法
                                                                           n
  精度(局部截斷誤差)：| f(x) - Pn(x) | = [f(ε)] / (n+1)!  ∏  ( x - x(k) ) 
 (注：其中[f(ε)]為f(ε)第n+1次求導的表達式)              k=0  
                                        
 《數值計算方法與算法》第二版 - 科學出版社 P19

  代碼維護：2007.04.18   pengkuny
**/

float Lagrange(float px[], float py[], int n, float x)
{//px,py:插值點(Xi,Yi) n:插值點個數  x:待計算的函數點
    float y = 0;
    for(int k=0; k<n; k++)//k控制Lagrange基函數序列
    {
        float tmp = 1;//tmp表示Lagrange基函數
        for(int i=0;i<n;i++)
        {
            if(i!=k)
            {
                tmp = tmp * (x-px[i])/(px[k]-px[i]);  //key step
            }
        }
        y = y + py[k]*tmp;
    }
    return y;
}

int main()
{
    float x;//插值
    float px[10];//已知(x0,y0),(x1,y1)
    float py[10];
    int n;//輸入已知插值組數

    cout<<"輸入插值組數:"<<endl;
    cin>>n;
    cout<<"輸入"<<n<<"組已知插值數（X,Y）"<<endl;
    for(int i=0; i<n; i++)
    {
        cin>>px[i]>>py[i];
    }
    cout<<"輸入插值:"<<endl;
    cin>>x;

    cout<<"Lagrange插值結果:"<<Lagrange(px, py, n, x)<<endl;

    system("pause");
    return 0;
}