#include "math.h"
#include "stdio.h"
#define iDim 2
#define iNum iDim+1
#define ep 0.001
#define v1 0.75              /*壓縮因子*/
#define v2 0.25              /*壓縮因子*/
#define u  2                 /*加長(zhǎng)因子*/

void select(float *p,int arraynum,int *imax,int *ismax,int *imin)
/*比較諸函數(shù)值的大小，選出最好點(diǎn)x(L)、最差點(diǎn)x(H)（對(duì)目標(biāo)函數(shù)而言）及次差點(diǎn) x(G)*/
{
 int i;
 *imax=*ismax=*imin=0;
 for(i=1;i<arraynum;i++)
 {
  if(*(p+i)>*(p+*ismax))
   {
    if(*(p+i)>*(p+*imax))
     {
      *ismax=*imax;
      *imax=i;
     }
    else
      *ismax=i;
   }
  else
    { if(*(p+i)<*(p+*imin))
      *imin=i;
    }
 }
}

float ffunction(float x,float y)
/*求目標(biāo)函數(shù)值*/
{
 float f;
 f=60-10*x-4*y+x*x+y*y-x*y;   
 return f;
}

initx(float a[iNum][iDim],int ifL,float h)
/*初始單純形產(chǎn)生*/
{int i,j;
 if(ifL!=0)
  {
   for(j=0;j<iNum-1;j++)
     a[0][j]=a[ifL][j];
  }
 for(i=1;i<iNum;i++)
   for(j=0;j<iDim;j++)
     {a[i][j]=a[0][j];}
 for(i=1;i<iNum;i++)
    a[i][i-1]=a[i][i-1]+h;
 return 0;
}


main()
{
 float x[iNum][iDim],fx[iNum];
 float xS[iDim],xE[iDim],xR[iDim],fS,fE,fR;
 int ifH,ifL,ifG;
 int i,j,count=0;
 float v,wucha;
 FILE *fp;

 if((fp=fopen("datef1.txt","w"))==NULL)
   {
    printf("File cannot be opened.\n");
    exit(1);
   }

for(i=0;i<iDim;i++)                  /*初始點(diǎn)*/
  x[0][i]=0;
initx(x,0,2.0);
/*x[0][0]=0;x[0][1]=0; x[1][0]=2;x[1][1]=0; x[2][0]=0;x[2][1]=2;   */
fprintf(fp,"x(0)=[%f,%f]  ",x[0][0],x[0][1]);
fprintf(fp,"x(1)=[%f,%f]  ",x[1][0],x[1][1]);
fprintf(fp,"x(2)=[%f,%f]\n",x[2][0],x[2][1]);

 fS=fE=fR=0;

 do{
    count++;
    for(i=0;i<iDim;i++)
      { xS[i]=xE[i]=xR[i]=0;}
    for(i=0;i<iNum;i++)
      {
       fx[i]=ffunction(x[i][0],x[i][1]);
	xR[0]=xR[0]+x[i][0];
	xR[1]=xR[1]+x[i][1];
       }

    select(fx,iNum,&ifH,&ifG,&ifL);
    if(fx[ifL]!=0)
	wucha=fabs(fx[ifH]-fx[ifL])/fabs(fx[ifL]);

    for(i=0;i<iDim;i++)
	xR[i]=xR[i]-x[ifH][i];
    for(i=0;i<iDim;i++)
	xR[i]=2*xR[i]/(iNum-1)-x[ifH][i];     /*求出反射點(diǎn)x(R)*/
    fR=ffunction(xR[0],xR[1]);

    if(fR>=fx[ifG])                          /*進(jìn)行步長(zhǎng)壓縮*/
      {
       v=fR<fx[ifH]?v1:v2;
	for(i=0;i<iDim;i++)
	  {xS[i]=(1-v)*fx[ifH]+v*xR[i];}
	fS=ffunction(xS[0],xS[1]);
      }//end if
    else    /*fR<fx[ifG]    進(jìn)行加大步長(zhǎng)*/
	 {
	 for(i=0;i<iDim;i++)
	    {xE[i]=(1-u)*x[ifH][i]+u*xR[i];}
	  fE=ffunction(xE[0],xE[1]);
	  if(fR>=fE)
	     {
	     for(i=0;i<iDim;i++)
		{xS[i]=xE[i];}
	      fS=fE;
	      }//end if
	  else   /*fR<fE    壓縮和加大步長(zhǎng)都不成立*/
	     {
	     for(i=0;i<iDim;i++)
	       {xS[i]=xR[i];}
		fS=fR;
	     }//end else
	  }//end else

    
    if(fS<fx[ifG])
       {
	   fprintf(fp,"\n第%d輪迭代:\n" , count);
	   fprintf(fp,"fL=%f,   x(L)=[%f,%f]  \n",fx[ifL],x[ifL][0],x[ifL][1]);
	   fprintf(fp,"fH=%f,   x(H)=[%f,%f]  ",fx[ifH],x[ifH][0],x[ifH][1]);
	   fprintf(fp,"ep=%f\n",wucha);
	   fprintf(fp,"fR=%f,   xR=[%f,%f]  ",fR,xR[0],xR[1]);
	   fprintf(fp,"fS=%f,   x(S)=[%f,%f]  ",fS,xS[0],xS[1]);
	   for(i=0;i<iDim;i++)
	     {x[ifH][i]=xS[i];}
	   fx[ifH]=fS;
       }//end if

    else      /*fS>=fx[ifG]*/
	{  if(fabs(fx[ifH]-fx[ifL])<=ep*fabs(fx[ifL]))
	    {
	      fprintf(fp,"\n最后結(jié)果:f=%f,   x=[%f,%f]  ",fx[ifL],x[ifL][0],x[ifL][1]);
             fprintf(fp,"fR=%f,   xR=[%f,%f]  ",fR,xR[0],xR[1]);
	      fprintf(fp,"fH=%f,   x(H)=[%f,%f]  ",fx[ifH],x[ifH][0],x[ifH][1]);
             fprintf(fp,"ep=%f\n",wucha);
	      break;
	     }//end if

	   else         /*fabs(fx[ifH]-fx[ifL])>ep*fabs(fx[ifL])*/
	  {
	      for(i=0;i<iNum;i++)        /*進(jìn)行單純形收縮*/
	        for(j=0;j<iDim;j++)
	          x[i][j]=(x[i][j]+x[ifL][j])/2;
	      fprintf(fp,"\n第%d輪迭代,進(jìn)行單純形收縮:\n",count);
             fprintf(fp,"fR=%f,   xR=[%f,%f]  ",fR,xR[0],xR[1]);
             fprintf(fp,"fH=%f,   x(H)=[%f,%f]\n",fx[ifH],x[ifH][0],x[ifH][1]);
	  }//end else

       }//end else

    if(count>50)
    break;
  }while(1);//end do_while

 fclose(fp);
 return 0;
 }