LagrangePolynomial

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

void viete(int,double*,double*);
double horner(int,double*,double,double*);
void lagrange(int,double*,double*,double*);

int main()
{
    FILE* in;
    FILE* out;
    double *x,*y,*l;
    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));
        l=(double*)malloc(n*sizeof(double));
        for(i=1; i<=n; i++)
            fscanf(in,"%lf",&x[i]);
        for(i=1; i<=n; i++)
            fscanf(in,"%lf",&y[i]);
        lagrange(n,x,y,l);
        for(i=0; i<=n-1; i++)
            fprintf(out,"l[%d]=%.15lf\n",i,l[i]);
        fprintf(out,"\n");
        free(x);
        free(y);
        free(l);
    }
    fclose(in);
    fclose(out);
    return 0;
}

void viete(int n,double* x,double* a)
{
    int i,j,e;
    double t;
    a[0]=1.0;
    for(i=1; i<=n; i++)
        a[i]=0.0;
    for(i=1; i<=n; i++)
        for(j=i; j>=1; j--)
            a[j]+=a[j-1]*x[i];
    e=n;
    for(i=0; i<(n+1)/2; i++)
    {
        t=a[i];
        a[i]=a[e];
        a[e]=t;
        e--;
    }
    for(i=0; i<=n; i++)
        if((n-i)%2==1)
            a[i]=-a[i];
}

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 lagrange(int n,double* x,double* y,double* l)
{
    int i,j;
    double* a;
    double* b;
    double d;
    a=(double*)malloc((n+1)*sizeof(double));
    b=(double*)malloc(n*sizeof(double));
    for(i=0; i<=n-1; i++)
        l[i]=0.0;
    viete(n,x,a);
    for(i=1; i<=n; i++)
    {
        horner(n,a,x[i],b);
        d=1.0;
        for(j=1; j<=n; j++)
            if(j!=i)
                d*=(x[i]-x[j]);
        for(j=0; j<=n-1; j++)
            l[j]+=(double)(y[i]/d)*b[j];
    }
    free(a);
    free(b);
}