Polynomial interpolation - Newton divided difference one dimensional array

#include<stdio.h>
#include<stdlib.h>

double horner(int,double*,double,double* );
void invhorner(int,double*,double,double*);
void newton(int,double*,double*,double*);

int main()
{
    FILE* in;
    FILE* out;
    double *x, *y, *a;
    int i,n;
    if((in=fopen("in.txt","r"))==NULL)
    {
        printf("Nie mozna otworzyc pliku do czytania \n");
        exit(1);
    }
    if((out=fopen("out.txt","w"))==NULL)
    {
        printf("Nie mozna otworzyc pliku do pisania \n");
        exit(1);
    }
    while(!feof(in))
    {
        fscanf(in,"%d",&n);
        x=(double*)malloc((n+1)*sizeof(double));
        y=(double*)malloc((n+1)*sizeof(double));
        a=(double*)malloc((n+1)*sizeof(double));
        for(i=0; i<n; i++)
            fscanf(in,"%lf",&x[i]);
        for(i=0; i<n; i++)
            fscanf(in,"%lf",&y[i]);
        newton(n,x,y,a);
        for(i=0; i<n; i++)
            fprintf(out,"a[%d]=%.15lf\n",i,a[i]);
        fprintf(out,"\n");
        free(x);
        free(y);
        free(a);
    }
    fclose(in);
    fclose(out);
    return 0;
}

double horner(int n,double* a,double c,double* b)
{
    int i;
    double r;
    r=a[n];
    for(i=n-1; i>=0; i--)
    {
        b[i]=r;
        r=r*c+a[i];
    }
    return r;
}

void invhorner(int n,double* b,double c,double* a)
{
    int k;
    double* u;
    u=(double*)malloc((n+1)*sizeof(double));
    for(k=0; k<=n; k++)
        u[k]=b[k];
    a[n+1]=b[n];
    for(k=0; k<=n; k++)
        if((k-1)<0)
            a[k]=-c*u[k];
        else
            a[k]=u[k-1]-c*u[k];
    free(u);
}

void newton(int n,double* x,double* y,double* a)
{
    int i,j;
    double* b;
    double* dd;
    b=(double*)malloc((n+1)*sizeof(double));
    dd=(double*)malloc((n+1)*sizeof(double));
    for(i=0; i<=n; i++)
    {
        b[i]=0.0;
        a[i]=0.0;
    }
    for(j=0; j<=n; j++)
        dd[j]=y[j];
    for(j=1; j<=n; j++)
        for(i=n; i>=j; i--)
            dd[i]=(double)(dd[i]-dd[i-1])/(x[i]-x[i-j]);
    b[0]=1.0;
    for(i=0; i<=n; i++)
    {
        for(j=0; j<=i; j++)
            a[j]+=dd[i]*b[j];
        invhorner(i,b,x[i],b);
    }
    free(b);
    free(dd);
}