/* usage; polyfit (order ) (datafile) */
/* compile; gcc -o polyfit polyfit.c -lm */
#include<stdlib.h>
#include<stdio.h>
#include<math.h>

#define MAX 100   /* maximum array size */

float a[MAX][MAX], b[MAX],x[MAX];
int i,j,k,l, row, col,M;
float X[MAX],Y[MAX];   /* data arrays , (x_k, y_k) */

void solve(float (*A)[MAX], float *B, float *X,int N); 
float determinant(float (*A)[MAX], int N);

float (*a_ptr)[MAX], *b_ptr, *x_ptr;
int N,NN,k;
FILE *fpt;             /* the data file to be read */



main(int argc , char **argv){

 if(argc !=3){
     printf("polyfit (order ) (datafile)\n");
     exit(1);
 }

NN=atoi(argv[1]);      /* order of polynomial fit */
fpt=fopen(argv[2],"r");
if(fpt == NULL)
printf("\n Error opening file");

/* read in M lines of the datafile */
 M=0;
 do{
fscanf(fpt,"%f %f\n",&X[M],&Y[M]);
 M=M+1;

 }while(!feof(fpt));
 
 printf("Read %d lines\n", M);

  for (i = 0; i <= NN; i++)
    {
      for (j = 0; j <= NN; j++)
	{
	  a[i][j]=0.0;
	}
    }
  /* set up a and b for polynomial fit of order N */
  for (j = 0; j <= NN; j++)
    {
      b[j]=0.0;
      x[j]=0.0;
      for (k = 0; k < M; k++)
	{
	  b[j]=b[j]+Y[k]* pow (X[k], j);
	}
    }

  for (i = 0; i <= NN; i++)
    {
      for (j = 0; j <= NN; j++)
	{
	  for (k = 0; k < M; k++)
	    a[i][j]=a[i][j] + pow (X[k], j + i);
	}
    }

b_ptr=&b[0];      /* initialize the first addresses */
x_ptr=&x[0];
a_ptr=&a[0][0];

solve(a_ptr,b_ptr,x_ptr,NN);  /* solve the equations */

 printf("y=a_0+a_1 x+...+a_n x^n\n");
for(i=0;i<=NN;i++)            /* output the results */
    printf("a_%d=%f\n",i,*(x_ptr+i));
printf("The determinant is %f\n", determinant(a,NN));

fclose(fpt);
}


void solve(float (*A)[MAX], float *B, float *X, int N) 
{
  /* the routine for solving the simultaneous equations */
  float temp, factor, sum;
  int m,n,k,u,p;
  
  for (k = 0; k <= N; k++)
    {
      m = k;
      u = 0;
      do
	{
	  u + 1;
	}
      while (*(*(A+k)+u+k) == 0.0);

      m = k + u;
      for (l = m + 1; l <= N; l++)
	{
	  if (*(*(A+k)+m) != 0.0)
	    {
	      factor =*(*(A+k)+l) / *(*(A+k)+m);
              

	      /* perform row operation on [a] */

	      for (j = 0; j <= N; j++)
		{
		  *(*(A+j)+l) = *(*(A+j)+l) - factor * (*(*(A+j)+m));
		}

/* perform row operation on [b] */
	      *(B+l) = *(B+l) - factor * (*(B+m));
	    }
	}
      for (j = 0; j <= N; j++)
	{
	  temp = *(*(A+j)+k);
	  *(*(A+j)+k) = *(*(A+j)+m);
	  *(*(A+j)+m) = temp;
	}
      /* do the switch for [b] */
      temp = *(B+k);
      *(B+k) = *(B+m);
      *(B+m) = temp; 
  }

  /* compute the solutions x[i] */
  *(X+N) = *(B+N) / (*(*(A+N)+N));
  for (p = 0; p <= N; p++)
    {
      sum = 0.0;
      for (u = N - p + 1; u <= N; u++)
	sum = sum + (*(*(A+u)+N-p)) * (*(X+u));
      *(X+N-p) = (*(B+N-p) - sum) / (*(*(A+N-p)+N-p));
    }

}



float determinant(float (*A)[MAX], int N) 
{
  /* the routine for solving the simultaneous equations */
  float temp, factor, product;
  int m,n,k,u,p;
  
  for (k = 0; k <= N; k++)
    {
      m = k;
      u = 0;
      do
	{
	  u + 1;
	}
      while (*(*(A+k)+u+k) == 0.0);

      m = k + u;
      for (l = m + 1; l <= N; l++)
	{
	  if (*(*(A+k)+m) != 0.0)
	    {
	      factor =*(*(A+k)+l) / *(*(A+k)+m);
              

	      /* perform row operation on [a] */

	      for (j = 0; j <= N; j++)
		{
		  *(*(A+j)+l) = *(*(A+j)+l) - factor * (*(*(A+j)+m));
		}
	    }
	}
      for (j = 0; j <= N; j++)
	{
	  temp = *(*(A+j)+k);
	  *(*(A+j)+k) = *(*(A+j)+m);
	  *(*(A+j)+m) = temp;
	}
  }

  /* compute the determinant */
 product=1.0;
 for(i=0;i<=N;i++)
 product=product*(*(*(A+i)+i));
return(product);

}