#include <algorithm>
using namespace std;
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <time.h>

const int maxn = 500 + 10;
const int n = 500;

double mat[maxn][maxn], b[maxn], x[maxn]; //儲存系數矩陣mat、常數向量b以及解向量x

//高斯消元法
void GaussElimination(double _mat[maxn][maxn], double b[maxn], double x[maxn])
{
	double mat[maxn][maxn];
	int i,j,k; //將矩陣和b向量粘合，成為n*(n+1)的矩陣方便運算
	for (i=1; i<=n; i++)
	{
		for (j=1; j<=n; j++)
			mat[i][j] = _mat[i][j];
		mat[i][n+1] = b[i];
	}
	for (i=1; i<=n; i++) //對每一列消元
	{
		for (j=i+1; j<=n; j++) //使用第i行對第i+1..n行進行消元
		{
			double pp = -mat[j][i] / mat[i][i];
			for (k=i; k<=n+1; k++)
				mat[j][k] += mat[i][k] * pp;
		}
	}
	for (i=n; i>=1; i--) //倒推求出x的解
	{
		for (j=i+1; j<=n; j++)
			mat[i][n+1] -= mat[i][j] * x[j];
		x[i] = mat[i][n+1] / mat[i][i];
	}
}

//列主元的高斯消元法
void GaussEliminationImprove(double _mat[maxn][maxn], double b[maxn], double x[maxn])
{
	double mat[maxn][maxn];
	int i,j,k; //將矩陣和b向量粘合，成為n*(n+1)的矩陣方便運算
	for (i=1; i<=n; i++)
	{
		for (j=1; j<=n; j++)
			mat[i][j] = _mat[i][j];
		mat[i][n+1] = b[i];
	}
	for (i=1; i<=n; i++) //對每一列消元
	{
		double maxp = 0; int dir; //選擇合適的列主元（絕對值最大項）
		for (j=i; j<=n; j++)
			if (fabs(mat[j][i])>maxp)
				maxp = fabs(mat[j][i]), dir = j;
		for (j=i; j<=n+1; j++)
			swap(mat[i][j], mat[dir][j]);

		for (j=i+1; j<=n; j++) //使用第i行對第i+1..n行進行消元
		{
			double pp = -mat[j][i] / mat[i][i];
			for (k=i; k<=n+1; k++)
				mat[j][k] += mat[i][k] * pp;
		}
	}
	for (i=n; i>=1; i--) //倒推求出x的解
	{
		for (j=i+1; j<=n; j++)
			mat[i][n+1] -= mat[i][j] * x[j];
		x[i] = mat[i][n+1] / mat[i][i];
	}
}

// Cholesky分解法
double l[maxn][maxn];
void CholeskyDecomposition(double a[maxn][maxn], double b[maxn], double x[maxn])
{
	int i,j,k;
	for (j=1; j<=n; j++) //根據公式計算l矩陣
	{
		l[j][j] = a[j][j];
		for (k=1; k<j; k++)
			l[j][j] -= l[j][k]*l[j][k];
		l[j][j] = sqrt(l[j][j]);
		for (i=j+1; i<=n; i++)
		{
			l[i][j] = a[i][j];
			for (k=1; k<j; k++)
				l[i][j] -= l[i][k]*l[j][k];
			l[i][j] /= l[j][j];
		}
	}
	double y[maxn]; //根據Ly=b求解y
	for (i=1; i<=n; i++)
	{
		y[i] = b[i];
		for (k=1; k<i; k++)
			y[i] -= l[i][k]*y[k];
		y[i] /= l[i][i];
	}
	for (i=n; i>=1; i--)  //根據L^Tx=y求解x
	{
		x[i] = y[i];
		for (k=i+1; k<=n; k++)
			x[i] -= l[k][i]*x[k];
		x[i] /= l[i][i];
	}
}

double t[maxn][maxn], d[maxn];
void CholeskyDecompositionImprove(double a[maxn][maxn], double b[maxn], double x[maxn])
{
	int i,j,k;
	d[1] = a[1][1];
	for (i=2; i<=n; i++) //根據公式分解A=LDL^T
	{
		for (j=1; j<i; j++)
		{
			t[i][j] = a[i][j];
			for (k=1; k<j; k++)
				t[i][j] -= t[i][k]*l[j][k];
			l[i][j] = t[i][j]/d[j];
		}
		d[i] = a[i][i];
		for (k=1; k<i; k++)
			d[i] -= t[i][k]*l[i][k];
	}
	double y[maxn];
	for (i=1; i<=n; i++) //根據Ly=b求解y
	{
		y[i] = b[i];
		for (k=1; k<i; k++)
			y[i] -= l[i][k]*y[k];
	}
	for (i=n; i>=1; i--)  //根據L^Tx=y求解x
	{
		x[i] = y[i]/d[i];
		for (k=i+1; k<=n; k++)
			x[i] -= l[k][i]*x[k];
	}
}

int main()
{
	int i,j;
	for (i=1; i<=n; i++) //給矩陣初始值
		for (j=1; j<=n; j++)
			mat[i][j] = 1.0/(i+j-1);
	for (int i=1; i<=n; i++)
	{
		b[i] = 0;
		for (j=1; j<=n; j++)
			b[i] += mat[i][j];
	}
	clock_t st;
	
	st = clock();
	GaussElimination(mat, b, x);
	printf("%.8lf\n",(double)(clock()-st)/CLOCKS_PER_SEC);

	st = clock();
	GaussEliminationImprove(mat, b, x);
	printf("%.8lf\n",(double)(clock()-st)/CLOCKS_PER_SEC);

	st = clock();
	CholeskyDecomposition(mat, b, x);
	printf("%.8lf\n",(double)(clock()-st)/CLOCKS_PER_SEC);

	st = clock();
	CholeskyDecompositionImprove(mat, b, x);
	printf("%.8lf\n",(double)(clock()-st)/CLOCKS_PER_SEC);
	/*for (i=1; i<=n; i++)
		printf("%.10lf\n",fabs(x[i]-1));*/
	return 0;
}