自己编的一个求无约束最优化解的单纯形

#include<iostream.h>
#include<math.h>
#define n 2
double F(double p[n])
{
	double f;
	f=60-10*p[0]-4*p[1]+p[0]*p[0]+p[1]*p[1]-p[0]*p[1];    //此函数为课本53页中的例2
	return f;
}
void main(void)
{
	double x[n][n+1];
	double a,b,c,l,esp;
	cout<<"请输入反射率:\t";
	cin>>a;
	cout<<"请输入收缩率:\t";
	cin>>b;
	cout<<"请输入扩张率:\t";
	cin>>c;
	cout<<"请输入正规单纯形边长:\t";
	cin>>l;
	cout<<"请输入所需精度:\t";
	cin>>esp;
	int i,j;
	for(i=0;i<n;i++)
	{
		cout<<"请输入初始点的第"<<i+1<<"维坐标:\t";
		cin>>x[i][0];
	}
	double p,q;
	p=(sqrt(n+1)+n-1)/(n*sqrt(2))*l;
	q=(sqrt(n+1)-1)/(n*sqrt(2))*l;
	for(i=0;i<n;i++)
	{
		for(j=1;j<n+1;j++)
		{
			if(i==j-1)
				x[i][j]=x[i][0]+p;
			else
				x[i][j]=x[i][0]+q;
		}
	}
	double y[n+1],t[n],xm[n],xr[n],xc[n],xe[n];
	double fl,fh,fr,fe,fc;
    int cl,ch;
    cout.setf(ios::fixed,ios::floatfield);      //按定点数格式输出，即不以科学计数法表示
da1:do
	{
		for(i=0;i<n;i++)
		{
			xm[i]=0;
			xr[i]=0;
			xc[i]=0;
			xe[i]=0;
		}
		double *p1,*q1;
	    cl=0;
	    ch=0;
	    for(j=0;j<n+1;j++)
		{
		   for(i=0;i<n;i++)
			  t[i]=x[i][j];
		   y[j]=F(t);
		}
		for(i=0;i<n;i++)
			t[i]=0;
		p1=&y[0];
		q1=&y[n];
		while(*p1!=*q1)
		{
			if(*p1>*q1)
			{
				p1++;
				cl++;
			}
			else
				q1--;
		}
		p1=&y[0];
		q1=&y[n];
		while(*p1!=*q1)
		{
			if(*p1<*q1)
			{
				p1++;
				ch++;
			}
			else
				q1--;
		}
		for(i=0;i<n;i++)
			t[i]=x[i][ch];
		fh=F(t);
		for(i=0;i<n;i++)
			t[i]=x[i][cl];
		fl=F(t);
		if(fabs(fh-fl)<esp*fabs(fl))
			break;
		else
		{
			for(i=0;i<n;i++)
			{
               for(j=0;j<n+1;j++)
				  xm[i]+=x[i][j];
			}
            for(i=0;i<n;i++)
			{ xm[i]=(xm[i]-x[i][ch])/n;}
			for(i=0;i<n;i++)
			{ xr[i]=xm[i]+(xm[i]-x[i][ch])*a;}
			fr=F(xr);
			if(fr<fh)
			{
				if(fr>fl)
				{
					for(i=0;i<n;i++)
						x[i][ch]=xr[i];
					goto da1;
				}
				else
				{
					for(i=0;i<n;i++)
						xe[i]=xm[i]+(xr[i]-xm[i])*c;
					fe=F(xe);
					if(fe<fr)
					{
						for(i=0;i<n;i++)
							x[i][ch]=xe[i];
						goto da1;
					}
					else
					{
						for(i=0;i<n;i++)
							x[i][ch]=xr[i];
					    goto da1;
					}
				}
			}
			else
			{
				for(i=0;i<n;i++)
					xc[i]=xm[i]+(x[i][ch]-xm[i])*b;
				fc=F(xc);
				if(fc<fh)
				{
					for(i=0;i<n;i++)
						x[i][ch]=xc[i];
				    goto da1;
				}
				else
				{
					for(j=0;j<n+1;j++)
						for(i=0;i<n;i++)
							x[i][j]=1/2*(x[i][j]+x[i][cl]);
					goto da1;
				}
			}
		}
	}
	while(fabs(fh-fl)<esp*fabs(fl));
    cout<<"最优点为:(";
 	for(i=0;i<n-1;i++)
		 cout<<x[i][cl]<<",";
	cout<<x[n-1][cl]<<")"<<'\n';
	cout<<"其值为:"<<fl<<'\n';
}