void erase(double pbar[], double prr[], double pr[])
{
	for (int i=0; i<=20; i++)
	{
		pbar[i]=0;
		prr[i]=0;
		pr[i]=0;
	}
}

void amoeba(double p[], double y[], int mp, int np, int ndim,
			double ftol, int& iter)
{
    int i,j,ihi,inhi,mpts,nmax = 20;
    double ypr,yprr,rtol,alpha = 1.0;
    double beta = 0.5;
    double gamma = 2.0;
    int itmax = 500;
    double pr[21], prr[21], pbar[21];
    mpts = ndim + 1;
    iter = 0;
    do
	{
		int ilo = 1;
        if (y[1] > y[2])
		{
            ihi = 1;
            inhi = 2;
		}
        else
		{
            ihi = 2;
            inhi = 1;
		}
        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])
				{
				    if (i != ihi)
					{
					    inhi = i;
					}
				}
			}
		}
        rtol = 2.0*fabs(y[ihi]-y[ilo])/(fabs(y[ihi])+fabs(y[ilo]));
        if (rtol < ftol)
		{
		    erase(pbar, prr, pr);
		    return;
		}
        if (iter == itmax)
		{
            cout<<" amoeba exceeding maximum iterations."<<endl;
            return;
		}
        iter = iter + 1;
        for (j = 1; j<=ndim; j++)
		{
            pbar[j] = 0.0;
		}
        for (i = 1; i<=mpts; i++)
		{
            if (i != ihi)
			{
                for (j = 1; j<=ndim; j++)
				{
                    pbar[j] = pbar[j] + p[(i-1)*np+j];
				}
			}
		}
        for (j = 1; j<=ndim; j++)
		{
            pbar[j] = pbar[j] / ndim;
            pr[j] = (1.0 + alpha) * pbar[j] - alpha * p[(ihi-1)*np+j];
		}
        ypr = famoeb(pr);
        if (ypr <= y[ilo])
		{
            for (j = 1; j<=ndim; j++)
			{
                prr[j] = gamma * pr[j] + (1.0 - gamma) * pbar[j];
			}
            yprr = famoeb(prr);
            if (yprr < y[ilo])
			{
                for (j = 1; j<=ndim; j++)
				{
                    p[(ihi-1)*np+j] = prr[j];
				}
                y[ihi] = yprr;
			}
            else
			{
                for (j = 1; j<=ndim; j++)
				{
                    p[(ihi-1)*np+j] = pr[j];
				}
                y[ihi] = ypr;
			}
		}
        else
		{
		    if (ypr >= y[inhi])
			{
			    if (ypr < y[ihi])
				{
				    for (j = 1; j<=ndim; j++)
					{
						p[(ihi-1)*np+j] = pr[j];
					}
				}
				y[ihi] = ypr;
			
				for (j = 1; j<=ndim; j++)
				{
					prr[j]=beta*p[(ihi-1)*np+j]+(1.0-beta)*pbar[j];
				}
				yprr = famoeb(prr);
				if (yprr < y[ihi])
				{
					for (j = 1; j<=ndim; j++)
					{
						p[(ihi-1)*np+j] = prr[j];
					}
					y[ihi] = yprr;
				}
				else
				{
					for (i = 1; i<=mpts; i++)
					{
						if (i != ilo)
						{
							for (j = 1; j<=ndim; j++)
							{
							  pr[j]=0.5*(p[(i-1)*np+j]+p[(ilo-1)*np+j]);
							  p[(i-1)*np+j] = pr[j];
							}
							y[i] = famoeb(pr);
						}
					}
				}
			}
			else
			{
				for (j = 1; j<=ndim; j++)
				{
					p[(ihi-1)*np+j] = pr[j];
				}
				y[ihi] = ypr;
			}
		}
    }while(1);
}