#include "stdio.h"
#include "conio.h"
#define N 4
static double x[N+1]={0.4,0.55,0.8,0.9,1};
static double y[N+1]={0.41075,0.57815,0.88811,1.02652,1.17520};
void main()
{
    double varx,f[N+1][N+1];
    int checkvalid ( double x[],int n);
    void Print(double x[], double y[],double f[][N+1]);
    void chashang ( double x[],double y[],double f[][N+1]);
    double compvalue (double f[][N+1],double x[],double varx);
    varx = 0.5;
    if(checkvalid(x,N)==1)
    {

        chashang(x,y,f);
        Print( x, y, f) ;
        varx = 0.5;
        printf("\n\n牛顿插值结果： P(%f)=%f\n",varx, compvalue(f,x,varx));
    }
    else printf("输入的插值节点的x值必须互异! ");
    getch();
}

void Print(double x[], double y[],double f[][N+1])
{
    int j,k;
    int n = 4;
    printf("i\t Xi\t  F(Xi)         1阶\t   2阶\t\t3阶\t    4阶\t\n");
    printf("-------------------------------------------------------------------------------\n");
    for(k = 0; k <= n; k++)
    {
       printf("x%d    %f   %f",k,x[k],y[k]);
       for(j=1;j<=k;j++)
           printf("    %f",f[0][j]);

       printf("\n");
    }
    printf("-------------------------------------------------------------------------------\n");
    printf("Ln(x)  =\n");
    printf("\t %f",y[0]);
    for(k = 1; k <= n; k++)
    {

       printf("\n\t+%f",f[0][k]);
       for( j = 1; j <= k; j++ )
       {
         printf("*(x-%f)",x[j-1]);
       }
    }
      
}

int checkvalid ( double x[],int n)
{
    int i , j;
    n = 4;

    for( i = 0;i <=n; i++)
      for( j = 0; j <= n; j++ )
        if((i!=j) && (x[i] == x[j]) )
           return -1;

    return 1;
}

void chashang ( double x[],double y[],double f[][N+1])
{
   int i,j,k;
  
   for( i=0;i<=N;i++)
   {
      k=0;
      for( j=i+1;k<=N;j++)
      {
        if(j==k+1)
           f[k][j]=(y[j]-y[k])/(x[j]-x[k]);
         
        else
           f[k][j]=(f[k+1][j]-f[k][j-1])/(x[j]-x[k]);
          
        k++;
      }
   }
}

double compvalue (double f[][N+1],double x[] ,double varx)
{
    double sum = y[0] , p = 1.0 , fx = 0.0;
    int j , n = 4 ;
   
    for( j = 1;  j <= n; j++ )
    {
       p = p*(varx - x[j-1]);
       fx = f[0][j]*p;
       sum = sum + fx;
    }
    return (sum);
}