using System;

public class gaussian{
	public matrix A;
	public gaussian(matrix a){A=a.copy();}
	}

public class nbody{

static Random rnd = new Random();

public static double rndexp(double x){
	double r = rnd.NextDouble();
	return -r.log()/x;
	}

public static matrix K = new matrix("0.5 0; 0 0.5");
public static int N=K.size1;
public static vector w1  = new vector(" 1;0");
public static vector w2  = new vector(" 0;2");
public static vector w12 = new vector("-1;1");

public static double[] matrix_elements(gaussian g1,gaussian g2){
	matrix A1=g1.A, A2=g2.A;
	matrix B=A1+A2;
	double detB=B.determinant();
	double overlap = (pi.pow(N)/detB).pow(1.5);
	double kinetic = 6*(B%A1*K*A2)*overlap;
	double beta;
	beta=1/(w1.T*Bi*w1);
	double E1=  (-2)*(2*beta.sqrt()/pi.sqrt())*overlap;
	beta=1/(w2.T*Bi*w2);
	double E1=  (-2)*(2*beta.sqrt()/pi.sqrt())*overlap;
	beta=1/(w12.T*Bi*w12);
	double E1=  (+1)*(2*beta.sqrt()/pi.sqrt())*overlap;
	double[] result = {kinetic+E1+E2+E3,overlap};
	return result;
	}

public static gaussian[] basis;
public static matrix N,H;
public static GEVD z;

public static double energy(gaussian[] basis){
	int n=basis.Length;
	for(int i=0;i<n;i++)
	for(int j=i;j<n;j++)
		{
		double[] r = matrix_elements(basis[i],basis[j]);
		H[i,j]=H[j,i]=r[0];
		N[i,j]=N[j,i]=r[1];
		}
	z=new GEVD(H,N);
	double e=z.e[0];
	return e;
	}	

public static int Main(){
	amoeba.function master = energy;
	//var x = new vector("0.2 1.0");
	//var x = new vector("0.2");
	//var x = new vector("0.2 1.0 5");
	var x = new vector("0.2 0.4 2 11");
	//var x = new vector("0.1 0.2 2 11 15");
	int K=x.size;
	var ss = new vector(K);
	for(int i=0;i<K;i++)ss[i]=0.05;
	basis = new gaussian[K];
	for(int i=0;i<K;i++)basis[i]=new gaussian(1);
	N = new matrix(K,K);
	H = new matrix(K,K);
	amoeba A = new amoeba(master,x,ss);
	double acc=0.001;
	A.downhill(acc);
	x=A.p[A.lo];
	x.print("x=","{0,9:g4}");
	}

}