using System;
using math=System.Math;

public static class cholesky_extension{

public static matrix cholesky(this matrix A, double tiny=1e-6)
{
	if(A.size1!=A.size2){return null;}
	int n=A.size1;
	var L = new matrix(n,n);
	for(int i=0;i<n;i++){
		double sum=0; for(int k=0;k<i;k++)sum+=L[i,k].pow(2);
		double aii=A[i,i]-sum;
		if(aii<tiny){return null;}
		L[i,i]=math.Sqrt(aii);
		for(int j=i+1;j<n;j++){
			sum=0; for(int k=0;k<i;k++)sum+=L[i,k]*L[j,k];
			L[j,i]=(A[j,i]-sum)/L[i,i];
		}
	}
	return L;
}

public static void L_solve(this matrix L, vector x){
	for(int i=0;i<x.size;i++){
		double sum=0; for(int k=0;k<i;k++)sum+=L[i,k]*x[k];
		x[i]=(x[i]-sum)/L[i,i];
		}
	}

public static void LT_solve(this matrix L, vector x){
	for(int i=x.size-1;i>=0;i--){
		double sum=0; for(int k=i+1;k<x.size;k++) sum+=L[k,i]*x[k];
		x[i]=(x[i]-sum)/L[i,i];
		}
	}

public static void solve(this matrix L, vector x){
	L.L_solve(x);
	L.LT_solve(x);
	}

public static matrix inverse(this matrix L){
	matrix M=new matrix(L.size1,L.size1);
	M.set_identity();
	for(int i=0;i<M.size2;i++) L.solve(M[i]);
	return M;
}

public static matrix L_inverse(this matrix L){
	matrix S=new matrix(L.size1,L.size1);
	S.set_identity();
	for(int i=0;i<L.size1;i++) L.L_solve(S[i]);
	return S;
}

}/* cholesky */