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

void print_matrix(int n, double a[][10])
{
   int x, y;
   for (x=0; x<n; x++)
   {
	 for (y=0; y<n; y++)
	 {
	    printf("%lf ", a[x][y]);
	 }
	 printf("\n");
   }
}

void print_column(int n, double b[])
{
	int i;
	for (i=0; i<n; i++)
	{
			printf("%lf\n", b[i]);
	}
}


void do_gauss(int n, double a[10][10], double x[10], double b[10])
{

int l[10] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

//real array (aij)nxn, (li)n;
//real array allocate (si)n;
int i=0;
int j=0;
int k=0;

double r=0.0;
double rmax=0.0;
double smax=0.0;
double xmult=0.0;
double tempval =0.0;
double s[10];

double sum=0.0;


for (i = 0; i<n; i++) 
{
	l[i] = i;
	smax = 0.0;
	
	for(j=0; j<n; j++)
	{
		//smax = max(smax, |aij|)
		tempval = (double)a[i][j];
				if(tempval < 0)
		{
			tempval = -tempval;
		}
		if(tempval > smax)
		{
			smax = tempval;
		}
		
	}

	s[i] = smax;

}


for (k=0; k<n-1; k++)
{
	rmax=0;

	for(i=k; i<n; i++)
	{
		//r=|ali,k/sli|
		r = fabs(a[l[i]][k] / s[l[i]]);
		
		
		if (r > rmax) 
		{
			rmax = r;
			j=i;
		}

	}

	//lj=lk;
    //lk=lj;
	tempval = l[j];
	
	l[j] = l[k];
	
	l[k] = tempval;
	
	for (i=(k+1); i<n; i++)
	{
		xmult = a[l[i]][k]/a[l[k]][k];
		a[l[i]][k] = xmult;
		
		for (j=(k+1); j<n; j++)
		{
		//ali,j = ali,j - (xmult)alk,j;
		
		a[l[i]][j] = a[l[i]][j] - xmult * a[l[k]][j];
		}
	}


}		

////second part of procedure...................

for(i = 0;  i < 10 ; i++)
{
	x[i] = 0;
}
i =0;
k =0; 
sum = 0;

for (k=0; k<n-1; k++)
	{
		for (i=(k+1); i<n; i++)
		{
			b[l[i]] = b[l[i]] - a[l[i]][k] * b[l[k]];
		}
	}



	x[n-1] = b[l[n-1]]/a[l[n-1]][n-1];

	for (i=(n-2); i >=0; i--)
	{
		sum = b[l[i]];
		
		for (j=(i+1); j<n; j++)
		{
		sum = sum - a[l[i]][j]*x[j];
		}
		

		x[i] = fabs(sum/a[l[i]][i]);
		if (x[i] == -0.00000)
		{
			x[i] = 0.00000;
		}

	}

}


void multiply_matrix(double a[][10], double b[][10], double result[][10], int n, int degree)
{
   

	int i = 0;
	int j = 0;
	int k = 0;
	for(i=0; i <= degree; i++)
	{
		for(j=0; j <= degree; j++)
		{
			for(k=0; k < n; k++)
			{
			result[i][j] = (a[k][j] * b[i][k]) + result[i][j];
			}
		}
	}



}

void multiply_column(double a[][10], double b[], double result[], int n, int degree)
{
	int i;
	int j;
	for(i=0; i<10; i++)
	{
		for(j=0; j<10; j++)
		{
			result[i] = (a[i][j]*b[j]) + result[i];
		}
	}
}


void transpose_matrix(double C[][10], double CT[][10])
{
	int i = 0;
	int j = 0;
	for(i=0; i<10; i++)
	{
		for(j=0; j<10; j++)
		{
			CT[i][j] = C[j][i];
		}
	}
}











int main()
{
	double A[10][10];
	double X[10];
	double B[10];
	double Y[10];
	double x_arr[10], y_arr[10];
	double C[10][10];
	double CT[10][10];
	double result[10][10];
	double resulty[10];
	double min, max, inc;	
	double a[10];
	int I[10];
	int iterations=0;
	
	double A_result[10][10];
	double B_result[10];

	int x = 0, y = 0, n=0;
	int degree=0;
	int i=0;
	int j=0;
	int k=0;
	double floor=0.0;
	double summation=0;
	
	for(i=0; i<10; i++)					//zero matrix
	{
		B[i] = 0;
		X[i] = i;
		Y[i] = 0;
		for(j=0; j<10; j++)
		{
			A[i][j] = 0;
			C[i][j] = 0;
			CT[i][j] = 0;
			result[i][j] = 0;
			resulty[i] =0;
			a[i] =0;
			I[i]=0;
		}
	}
	

	printf("Analysis #7 Devin Scott CS 407 \n");


	printf("please enter the degree\n");
	scanf("%d",&degree);
	
	printf("please enter number of points: \n");
	scanf("%d",&n);
	
	printf("\nplease enter data for %d points in the form X space Y <ENTER>\n",n);
	for (x = 0; x < n; x++)
	{
		printf("X%d Y%d: \n",x+1,x+1);
		scanf("%lf %lf",&A[x][0], &A[x][1]);
	}
	
	printf("what is the minimum x? \n");
	scanf("%lf",&min);
	
	printf("what is the increment? \n");
	scanf("%lf",&inc);
	
	printf("what is the max x? \n");
	scanf("%lf",&max);

	iterations = (fabs(min - max))/(inc);
	
	//make C matrix///////////////////////////////////

	
	for (x = 0; x < n; x++)
	{
		B[x]=A[x][1];
		C[x][degree]=1;
		for(y = degree-1; y > -1; y--)
		{
			C[x][y]=pow(A[x][0], degree-y);
		}
	}
	


	printf("\n printing the C matrix only\n");
	print_matrix(n,C);
	printf("\n the Y matrix: (B matrix)\n");
	print_column(n,B);
	printf("C Transpose is: \n");
	transpose_matrix(C,CT);
	print_matrix(n,CT);
	multiply_matrix(C, CT, result, n, degree);
	multiply_column(CT, B, resulty, n, degree);	
	printf("\n printing CT *C\n");
	print_matrix(n,result);
	printf("\n printing CT *Y\n");
	print_column(n,resulty);
	printf("\n\n");
	do_gauss(degree+1, result, a, resulty);

	printf("min X was: %lf \n",min);
	printf("max X was: %lf \n",max);\
	printf("increment was: %lf \n",inc);
	printf("iterations was: %d \n",iterations);
	


	printf("for n = %d , the result of gaussian elimination is: \n",n);
	//print the result!
	for (x = 1; x <n+1; x++)
	{
		
		printf("a [ %d ] = %lf\n",x,a[x-1]);
		
	}

	
	printf("X    Y\n");
	for( i = 0; i <= iterations; i++)
	{
		floor = min + i*inc;
		for(j = 0; j <= degree; j++)
		{
			summation = a[j]*pow(floor,degree-j) + summation;
		}
	
		printf("%lf    %lf\n", floor, summation);
		summation = 0;
	}

	
	return 0;
}