using System;
using static System.Console;
public static partial class roots{
//public static int nsteps; /* already defined in newton */

public static vector broyden
(Func<vector,vector> f, vector x, double eps=1e-3){
	vector fx=f(x),z,fz;
	matrix J = jacobian(f,x,fx);
	var qrj = new givensQR(J);
	matrix B=qrj.inverse();
	nsteps=0;
	do{	nsteps++;
		vector Dx=-B*fx;
		double lambda=2;
		do{                  /* line search */
			lambda/=2;
			z=x+Dx*lambda;
			fz=f(z);
			if(lambda<1.0/32){ /* reset */
				J=jacobian(f,z,fz);
				qrj=new givensQR(J);
				B=qrj.inverse();
				break;
				}
		}while(fz.norm()>(1-lambda/2)*fx.norm());
		vector dx=z-x;
		vector df=fz-fx;

		vector c=(dx-B*df)/(dx % df); B.update(c,dx);
		//vector c=(dx-B*df)/(df%df); B.update(c,df);
		//vector c=(dx-B*df)/(dx%(B*df)); B.update(c,B*dx);
		
		x=z; fx=fz;
	}while(fx.norm()>eps && nsteps<nstepsmax);
	return x;
}//broyden

}//class