#include <stdio.h>
#include <math.h>

/*
  plot the output with the command

  plot "ddif.dat" index 0 with lines,"ddif.dat" index 1 ps 4 pt 5

 */

#define N 4
double x[N]= { 1,2,4,5 };
double fx[N]= { 0,4,-2,6 };
double dd[N][N];

void matprint(double A[N][N]){
    for(int i=0;i<N;i++){
        for(int j=0;j<N;j++){
            printf("%g ",A[i][j]);
        }
        printf("\n");
    }
}

double P(double t){
    double s=dd[0][N-1];
    for(int j=N-2;j>=0;j--){
        s=dd[0][j]+(t-x[j])*s;
    }
    return s;
}
int main(){
    for(int i=0;i<N;i++) dd[i][0]=fx[i];
    for(int j=1;j<N;j++){
        for(int i=0;i<N-j;i++) {
            dd[i][j]=(dd[i+1][j-1]-dd[i][j-1])
                     /(x[i+j]-x[i]);
        }
    }
//    matprint(dd);
    double xmin=0,xmax=6;
    int points=10;
    double h=(xmax-xmin)/points;
    for(int i=0;i<=points;i++){
        double xi=xmin+h*i;
        printf("%20.13e %20.13e\n",xi,P(xi));
    }
    printf("\n\n");
    for(int i=0;i<N;i++){
        printf("%20.13e %20.13e\n",x[i],P(x[i]));
    }
    return 0;
}