#include<stdio.h>
#include<math.h>
#define new(T,N) ((T*)malloc(sizeof(T)*N))

void inverse(const int n, double A[][n], double InvA[][n]);

double lstsqr(int n,const char * filename)
{
	FILE * f;
	char ch;
	int i,j,k,q;

	double InvX[n][n];
	
	double X[n][n],Y[n],p[n],detX,dydx;
	double *x,*y,*xlo,*xhi,*ylo,*yhi,*dy;

	q = 0;
	f = fopen(filename,"r");
	while((ch=fgetc(f))!=EOF){
		if(ch == '\n')
		{
		  q++;
		}
	}
	fclose(f);

	x   = new(double,q);	
	y   = new(double,q);	
	xhi = new(double,q);	
	xlo = new(double,q);	
	yhi = new(double,q);	
	ylo = new(double,q);	
	dy  = new(double,q);	
	
	for(i=0;i<n;i++){
		Y[i] = 0.0;
		p[i] = 0.0;
	 for(j=0;j<n;j++){
		 X[i][j]    = 0.0;
	 	 InvX[i][j] = 0.0;
	 }
	}

	f = fopen(filename,"r");

	for(i=0;i<q;i++){
	fscanf(f,"%lf %lf %lf",&x[i], &xhi[i], &xlo[i]);
	fscanf(f,"%lf %lf %lf",&y[i], &yhi[i], &ylo[i]);
	}

	for(i=0;i<q;i++){
		if(i!=q){
		dydx = (y[i]-y[i+1])/(x[i]-x[i+1]);
		}
		else{	
		dydx = (y[i-1]-y[i])/(x[i-1]-x[i]);
		}
		dy[i] = sqrt((pow(yhi[i],2)+pow(ylo[i],2))*0.5 + pow(dydx,2)*(pow(xhi[i],2)+pow(xlo[i],2))*0.5);
	}

         for(i=0;i<n;i++){
	  for(k=0;k<n;k++){
	   for(j=0;j<q;j++){
		  X[i][k] = X[i][k] + dy[j]*pow(x[j],i+k);
		  if(k==0){
		   Y[i]    = Y[i] + dy[j]*y[j]*pow(x[j],i);
		  }
	   }	
	  }
	 } 
	
	 free(x); 
	 free(y); 
	 free(xhi); 
	 free(xlo); 
	 free(yhi); 
	 free(ylo);
	 free(dy);
	 
	 inverse(n,X,InvX);
	  
	 for(i=0;i<n;i++)
	  for(j=0;j<n;j++){
	   	p[i] = p[i] + Y[j]*InvX[j][i];
	  }
	 
	for(i=0;i<n;i++){
		printf("p[%d] = %lf\n",i,p[i]);
	}
	
	return 3.14;
}

int main(int argc,char **argv)
{
	double out;

	out = lstsqr(3,"datafile");

}