/**
  牛頓插值法

  牛頓遞推公式：N_k(x) = N_k-1(x) + t_k(x) / f [x₀,x₁,..,x_k]

  基函數遞推公式： t_i (x)  = (x - x _i ) * t_i-1 (x) 

  屬性：插值計算法
  
   《數值計算方法與算法》第二版 - 科學出版社 P25

  代碼維護：2007.04.18   pengkuny
**/

#include <iostream>
using namespace std;

float NewtonInterpolate(float px[], float py[], int n, float x)
{//px,py;插值點(Xi,Yi) n:插值點個數  x:待計算的點
    float *g = new float[n+1];

    //計算差商表個g[k],k=0,1,,n
    for (int i=0; i<n+1; i++)
    {
        g[i] = py[i];
    }
    for (int k=1; k<n+1; k++)
    {
        for (int i=n; i>=k; i--)
        {
            g[i] = (g[i] - g[i-1]) / (px[i] - px[i-k]);//g[i]用來暫時存放f[X(i-k), X(i)]
        }
    }

    float t = 1;//基函數
    float newton = py[0];
    for (int i=1; i<n; i++)
    {
        t = (x - px[i-1]) * t;//基函數的遞推式
        newton = newton +  t * g[i];
    }

    delete g;
    return newton;
}

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<<"牛頓插值結果:"<<NewtonInterpolate(px, py, n, x)<<endl;

    system("pause");
    return 0;
}