using System;

public class qnewton
{

public delegate double function(vector x);

public static vector gradient(function f, vector x, double dx=1e-6){
	double fx=f(x);
	vector dfdx = new vector(x.size);
	for(int i=0;i<x.size;i++){
		x[i]+=dx;
		dfdx[i]=(f(x)-fx)/dx;
		x[i]-=dx;
		}
	return dfdx;
	}

public static int minimize
(function f, vector p, double dx=1e-6, double eps=1e-3){
	vector x=p, s, z, Dx, dfdx=gradient(f,x,dx);
	double fz, fx=f(x);
	matrix H1=new matrix(x.size,x.size);
	H1.set_unity();
	int nsteps=0;
	do{
		nsteps++;
		Dx=-H1*dfdx; s=Dx*2;
		do{
			s/=2; z=x+s; fz=f(z); /* backtracking */
			if( Math.Abs(fz) < Math.Abs(fx)+.1*s.dot(dfdx) ) break;
			if(s.norm() < dx) {H1.set_unity(); break;}
		}while(true);
		vector dfdx_z=gradient(f,z,dx), y=dfdx_z-dfdx, u=s-H1*y;
		H1+=matrix.outer(u,u)/u.dot(y);
		x=z; fx=fz; dfdx=dfdx_z;
		} while(Dx.norm()>dx && dfdx.norm()>eps);
	vector.copy(p,x);
	return nsteps;
	}
}