#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "nrutil.h"
#define FTOL 1.0e-5
#define TINY 1.0e-10
#define NMAX 5000
#define NDIM 5

#define GET_PSUM \
		   for (j=1;j<=ndim;j++) {\
		   for (sum=0.0,i=1;i<=mpts;i++) sum += p[i][j];\
		   psum[j]=sum;}
#define SWAP(a,b) {swap=(a);(a)=(b);(b)=swap;}

void amoeba(double **p, double y[], int ndim, double ftol,
    double funk(double []), int *nfunk)
{
    double amotry(double **p, double y[], double psum[], int ndim,
        double (funk)(double []), int ihi, double fac);
    int i,ihi,ilo,inhi,j,mpts=ndim+1;
    double rtol,sum,swap,ysave,ytry,*psum;

    psum=dvector(1,ndim);
    *nfunk=0;
    GET_PSUM
    for (;;) {
        ilo=1;
        ihi = y[1]>y[2] ? (inhi=2,1) : (inhi=1,2);
        for (i=1;i<=mpts;i++) {
           if (y[i] <= y[ilo]) ilo=i;
           if (y[i] > y[ihi]) {
              inhi=ihi;
              ihi=i;
	    } else if (y[i] > y[inhi] && i != ihi) inhi=i;
           }
	rtol=2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo])+TINY);
	if (rtol < ftol) {
	   SWAP(y[1],y[ilo])
		for (i=1;i<=ndim;i++) SWAP(p[1][i],p[ilo][i])
		 break;
	}
	if (*nfunk >= NMAX) nrerror("NMAX exceeded");
	*nfunk +=2;
	ytry=amotry(p,y,psum,ndim,funk,ihi,-1.0);
	if (ytry<=y[ilo])
	   ytry=amotry(p,y,psum,ndim,funk,ihi,2.0);///////////////////
	else if (ytry >= y[inhi]) {
	     ysave=y[ihi];
	     ytry=amotry(p,y,psum,ndim,funk,ihi,0.5);
	     if (ytry >= ysave) {
		for (i=1;i<=mpts;i++) {
		    if (i != ilo) {
			for (j=1;j<=ndim;j++)
			    p[i][j]=psum[j]=0.5*(p[i][j]+p[ilo][j]);
			y[i]=funk(psum);
		    }
		}
		*nfunk += ndim;
		GET_PSUM
	      }
	    } else --(*nfunk);
	}

	free_dvector(psum,1,ndim);
}
// #include "nrutil.h"

double amotry(double **p,double y[], double psum[], int ndim,
    double funk(double []), int ihi, double fac)
{
    int j;
    double fac1,fac2,ytry,*ptry;

    ptry=dvector(1,ndim);
    fac1=(1.0-fac)/ndim;
    fac2=fac1-fac;
    for (j=1;j<=ndim;j++) ptry[j]=psum[j]*fac1-p[ihi][j]*fac2;
    ytry=funk(ptry);
    if (ytry < y[ihi]) {
	y[ihi]=ytry;
	for (j=1;j<=ndim;j++) {
		psum[j] +=ptry[j]-p[ihi][j];
		p[ihi][j]=ptry[j];
	    }
	}
	free_dvector(ptry,1,ndim);
//printf("f=%lf,",ytry);
	return ytry;
}



double myFun(double x[])
{
	// x=dvector(1, NDIM);
double f;
	f = (x[1]-2.54)*(x[1]-2.54) + (x[2]+7.5)*(x[2]+7.5) + 
(x[3]+1.22)*(x[3]+1.22) + (x[4]-6.34)*(x[4]-6.34)+ 
(x[5]+8.78)*(x[5]+8.78) + 239.78984346;

printf("f=%lf\n",f);
	return f;
//	return (x[1]*x[1] + x[2]*x[2]);

}



void  main()
{
	double **p;
	double *y;
//	double ftol;
	int nfunk;
	int i,j;

	y=dvector(1, NDIM + 1);

	p=dmatrix(1, NDIM + 1, 1, NDIM);
	for(i=1; i<= NDIM+1 ; i++)
	{

		for(j=1; j<=NDIM; j++)
			p[i][j]=rand()/1000;

	}
	//p[1][1] = 1.0;
	//p[1][2] = 8.0;
	//p[2][1] = 1.0;
	//p[2][2] = 3.0;
	//p[3][1] = 4.0;
	//p[3][2] = -2.0;

	for(i=1; i<= NDIM+1 ; i++)
	{
		y[i] = myFun(p[i]);
printf("y=%lf\n",y[i]);
	}

	amoeba(p, y, NDIM, FTOL, myFun, &nfunk);
	
	printf("Thank God it's over! nfunk = %d\n", nfunk);
	
	for(i=1; i<= NDIM+1 ; i++)
	{
		printf("y[%d] = %f\n", i, y[i]);
//printf("ran=%d",rand());
		//for(j=1; j<=NDIM; j++)
		//	printf("p[%d][%d] = %f\n", i, j, p[i][j]);

	}


	free_dvector(y, 1, NDIM + 1);
	free_dmatrix(p, 1, NDIM + 1, 1, NDIM);
}