using System;
using static System.Console;
using static System.Math;

public static class steep{

public static (vector,vector,double,double)
makestep(matrix H,matrix N,vector x,matrix iHess){
	vector Hx=H*x,Nx=N*x;
	double xHx=x%Hx;
	double xNx=x%Nx;
	double E=xHx/xNx;
	vector grad=(Hx-E*Nx)*(2/xNx);
	vector dx=iHess*grad;
//	vector dx=grad;
	double dxHx=dx%Hx,dxNx=dx%Nx,dxHdx=dx%(H*dx),dxNdx=dx%(N*dx);
	double A=xHx,B=2*dxHx,C=dxHdx;
	double a=xNx,b=2*dxNx,c=dxNdx;
	Func<double,double> f = z=>(A+B*z+C*z*z)/(a+b*z+c*z*z);
	double D=A*A*c*c+((B*B-2*A*C)*a-A*B*b)*c+A*C*b*b-B*C*a*b+C*C*a*a;
	double s1=-(Sqrt(D)+A*c-C*a)/(B*c-C*b);
	double s2= (Sqrt(D)-A*c+C*a)/(B*c-C*b);
	double f0=f(0),f1=f(s1),f2=f(s2);
	double s;
	if(f1<f0 && f1<f2) {
//Error.WriteLine("s1");
		s=s1;
		}
	else if(f2<f0 && f2<f1) {
Error.WriteLine("uups: s2");
		s=s2;
		}
	else {
Error.WriteLine("uups: too bad...");
		s=-1.0/128;
		}
	return (grad,s*dx,f(s),xNx);
	}

public static (double,vector) find
(matrix H, matrix N, vector x, double acc=1e-6){
	vector gradx,gradz,Dx;
	matrix iHes=matrix.id(H.size1);
	double E,xNx;
	int nsteps=0;
	(gradx,Dx,E,xNx)=makestep(H,N,x,iHes);
	do{
		nsteps++;
		x+=Dx;
		(gradz,Dx,E,xNx)=makestep(H,N,x,iHes);
		vector y=gradz-gradx;
		vector u=Dx-iHes*y;
		double uy=u%y,sy=Dx%y;
		double gamma=uy/sy/2;
		vector a=(u-gamma*Dx)/sy;
		iHes.update(a,Dx); iHes.update(Dx,a);
		gradx=gradz;
Error.WriteLine($"nsteps={nsteps} E={E}");
	}while(gradx.norm()>acc*x.size);
	return (E,x);
}//geneig

}//class


class main{
static void Main(string[] argv){
	int n=8;
	if(argv.Length>0)n=int.Parse(argv[0]);
	matrix A = new matrix(n,n);
	matrix N = matrix.id(n);
	var rndA=new System.Random(1);
	var rndN=new System.Random(1);
	for(int i=0;i<n;i++)
	for(int j=i;j<n;j++){
		A[i,j]=A[j,i]=rndA.NextDouble();
		N[i,j]=N[j,i]=rndN.NextDouble();
		}
	N=matrix.id(n)+N.T*N;
	//N.setid();
if(n<5)A.print("A=");
if(n<5)N.print("N=");
	vector v=new vector(n);
	for(int i=0;i<n;i++)v[i]=rndA.NextDouble()-0.5;
	v/=Sqrt(v%(N*v));
	(double e,vector x) = steep.find(A,N,v);
	WriteLine($"n={n} e={e}");
	var Ax=A*x;
	var eNx=e*(N*x);
	double miss=(Ax-eNx).norm();
	//(Ax-eNx).print("Ax-eNx");
	WriteLine($"|A*x-e*N*x|={miss}");
if(n<5)x        .print("x   =");
if(n<5)(A*x)    .print("Ax  =");
if(n<5)(e*(N*x)).print("eNx =");
}
}