C語言 例題1-8 已知下面5點座標 建立牛頓差除表 並計算 P(1.5)之值

C語言 例題1-8  已知下面5點座標

x       f(x)
===========
0        0
1        -3
2        0
3        15
4        48
===========   
建立牛頓差除表 並計算 P(1.5)之值

程式
#include<stdio.h>
#define MAXN 100
#define ORDER 4

main()
{
    float ax[MAXN+1], ay [MAXN+1], diff[MAXN+1][ORDER+1], nr=1.0, dr=1.0,x,p,h,yp;
    int n,i,j,k;
    //printf("\nEnter the value of n:\n");
    scanf("n=%d",&n);

    //printf("\nEnter the values in form x,y:\n");
    for (i=0;i<=n;i++)
        scanf("%f %f",&ax[i],&ay[i]);
    //printf("\nEnter the value of x for which the value of y is wanted: \n");
    scanf("%f",&x);
   
    h=ax[1]-ax[0];

    //now making the difference table
    //calculating the 1st order of differences
    for (i=0;i<=n-1;i++)
    {   
        diff[i][1] = ay[i+1]-ay[i];
    }
    //now calculating the second and higher order differences
    for (j=2;j<=ORDER;j++)
    {
        for(i=0;i<=n-j;i++)
        {
            diff[i][j] = diff[i+1][j-1] - diff[i][j-1];

        }
    }   
    //now finding x0
    i=0;
    while (!(ax[i]>x))
        i++;

    //now ax[i] is x0 and ay[i] is y0
    i--;
    p = (x-ax[i])/h;
    yp = ay[i];

    //now carrying out interpolation
    for (k=1;k<=ORDER;k++)
    {
        nr *=p-k+1;
        dr *=k;
        yp +=(nr/dr)*diff[i][k];
    }
    printf("\nWhen x = %6.1f, corresponding y = %0.4f\n",x,yp);
}


輸入資料
n=4
0.0    0.0
1.0    -3.0
2.0    0.0
3.0    15.0
4.0    48.0
1.5


輸出畫面
When x =    1.5, corresponding y = -2.6250