using System;
using static System.Console;
using static System.Math;
public static partial class rayleigh{

public static double pow (this double x,double y){return Pow(x,y);}
public static double pow (this double x,int n){return Pow(x,n);}

public static (double,vector,double) gengrad(matrix H,matrix N,vector x){
	vector Hx=H*x,Nx=N*x;
	double xNx=x%Nx;
	double r=x%Hx/xNx;
	vector grad=(Hx-r*Nx)*(2/xNx);
	return (r,grad,xNx);
	}

public static (double,vector) geneig
(matrix H, matrix N, vector x, double acc=1e-6){
	matrix B=matrix.id(x.size);
	x/=Sqrt(x%(N*x));
	(double fx,vector gx,double xNx)=gengrad(H,N,x);
	int nsteps=0;
	while(nsteps<3*x.size){
		if(gx.norm()<acc*x.size){
			Error.Write($"broyden: job done: |gx|<acc\n");
			break;
			}
		nsteps++;
		vector Dx=-B*gx;
		vector z,gz;
		double fz,zNz,lambda=1.0;
		while(true){// backtracking linesearch
			z=x+Dx*lambda;
			(fz,gz,zNz)=gengrad(H,N,z);
			if(fz<fx){
Error.WriteLine($"good step # {nsteps:d4} lambda={lambda:f5} e={fz}");
				break; // good step
				}
			if(lambda<1.0/128){
Error.WriteLine($"BAD  step # {nsteps} lambda={lambda} deltaf={fz-fx}");
				//B.setid();
				break; // accept anyway
				}
			lambda/=2;
		}
		vector s=z-x, y=gz-gx, u=s-B*y;
		double uy=u%y;

//		B.update(u,u,1/uy); //SR1

// symmetric Broyden update
		double sy=s%y, gamma=uy/sy/2;
		vector a=(u-gamma*s)/sy;
		B.update(a,s); B.update(s,a);

// BFGS update
//		vector By=B*y;
//		B.update(s,s,(s%y+y%(B*y))/(s%y).pow(2));
//		B.update(s,By,-1/(s%y));
//		B.update(By,s,-1/(s%y));
		

// DFP update
//		vector By=B*y;
//		B.update(s,s,1/(s%y));
//		B.update(By,By,-1/(y%By));

		x=z/Sqrt(zNz);
		(fx,gx,xNx)=gengrad(H,N,x);
	}
	Error.WriteLine($"H.size={H.size1} nstpes={nsteps}");
	return (fx,x);
}//SR1

}//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();
		}
if(n<5)A.print("A=");
	N=matrix.id(n)+N.T*N;
	//A.print("A=");
if(n<5)N.print("N=");
	//N.setid();
	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) = rayleigh.geneig(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=");
}
}