#include <math.h>
#include <iostream>
#include <iostream.h>
#define NRANSI
#include "nrutil.h"
#define ITMAX 900
#define EPS 1.0e-10
#define FREEALL free_nrvector(xi,0,n-1);free_nrvector(h,0,n-1);free_nrvector(g,0,n-1);

void frprmn(float p[], int n, float ftol, int *iter, float *fret,
	float (*func)(float []), void (*dfunc)(float [], float []))
{
	void linmin(float p[], float xi[], int n, float *fret,
		float (*func)(float []));
	int j,its;
	float gg,gam,fp,dgg;
	float *g,*h,*xi;

	// Debug 
	//cout << "frprmn 0 : p=";
	//for(int i=0;i<n;i++){cout << p[i] << " ";}cout << "\n";	
	
	g=nrvector(0,n-1);
	h=nrvector(0,n-1);
	xi=nrvector(0,n-1);
	fp=(*func)(p);
	(*dfunc)(p,xi);
	for (j=0;j<n;j++) {
	  g[j] = -xi[j];
	  xi[j]=h[j]=g[j];
	}

	////float prevret=+999;
	
	for (its=1;its<=ITMAX;its++) {
	  
	  *iter=its;
	  linmin(p,xi,n,fret,func);
	  
	  //cout << "its=" << its << " f=" << (*fret)<<"\n";	  
	  if (2.0*fabs(*fret-fp) <= ftol*(fabs(*fret)+fabs(fp)+EPS)) {
	    FREEALL
	    return;
	  }	  
	  fp= *fret;

	  ////cout << "iter " << its << "df = "<< prevret-(*fret) << "\n";
	  ////prevret = (*fret);
	  
	  (*dfunc)(p,xi);
	  dgg=gg=0.0;
	  for (j=0;j<n;j++) {
	    gg += g[j]*g[j];
	    dgg += (xi[j]+g[j])*xi[j];
	  }
	  if (gg == 0.0) {
	    FREEALL
	      return;
	  }
	  gam=dgg/gg;
	  for (j=0;j<n;j++) {
	    g[j] = -xi[j];
	    xi[j]=h[j]=g[j]+gam*h[j];
	  }
	}
	//cout<<("Too many iterations in frprmn\n");
	return;	
}
#undef ITMAX
#undef EPS
#undef FREEALL
#undef NRANSI