#include <conio.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define XuYimin_Debug_1
#define N 3  /*N维空间*/

double xq[N]={0,0,0};

double fun(double x[N])
{
  double y;
 /*目标函数f(X)=
(x1-4)^2+(x2-3)^2+1.5*x3^2+(x1-x2+x3-1)^2+sin(x1+x2+x3-7)*sin(x1+x2+x3-7)
可根据需要修改，但需要同时修改维数N值*/
  y=pow((x[0]-4),2)+pow((x[1]-3),2)+1.5*x[2]*x[2]
    +pow(x[0]-x[1]+x[2]-1,2)+pow(sin(x[0]+x[1]+x[2]-7),2);
  return y;
}

double fd(double x[N],double d[N],double alf)
{
  int i;
  double y,t[N];
  for(i=0;i<N;i++)
    t[i]=x[i]+alf*d[i];
  y=fun(t);
  return y;
}

double alf(double x[N],double d[N])
{
  double h0=0.1,a0=0.0,ep1;
  double a1,a2,a3,y1,y2,y3,c1,c2,ap,yp,ax,yx,h;
  ep1=pow(0.1,12);
  a1=a0,h=h0;a2=h,y1=fd(x,d,a1),y2=fd(x,d,a2);
  if(y2>y1)
    {
      h=-1*h,a3=a1,y3=y1;
      loop1:a1=a2,y1=y2,a2=a3,y2=y3;
      loop2:a3=a2+h,y3=fd(x,d,a3); if(y3<y2) {h=2*h;goto loop1;}
    }
  else  goto loop2;
  if( fabs(y2)<ep1 )  goto loop5;
  loop3:
    c1=(y3-y1)/(a3-a1);
    c2=((y2-y1)/(a2-a1)-c1)/(a2-a3);
    if(c2==0) goto loop4;
    ap=0.5*(a1+a3-c1/c2);
    yp=fd(x,d,ap);
  while(!fabs((y2-yp)/y2)<ep1)
    {
      if((ap-a2)*h>0)
        {
          if(y2>=yp)  a1=a2,y1=y2,a2=ap,y2=yp;
          else        a3=ap,y3=yp;
        }
      else
        {
          if(y2>=yp)  a3=a2,y3=y2,a2=ap,y2=yp;
          else        a1=ap,y1=yp;
        }
      goto loop3;
    }
  if(y2<yp)  ax=a2,yx=y2;
  else       ax=ap,yx=yp;
  return ax;
  loop4:return a2;
  loop5:return 0;
}

double maxt(double del[N+1])
{
  int i;
  double y=del[0];
  for (i=1;i<N+1;i++)
    if(del[i]>y)
      y=del[i];
  return y;
}

void powell()
{
  int i,j,id,pd[3],xym=1,m;
  double al,esp,ALF[N+1],x[2],pa[N+1],X[N][N+1],del[N+1],d[N][N+1],
      F[3],fu[N+1],te[N],tem[N],dk[N],xk[N],delm;
  esp=pow(0.1,10);x[0]=xq[0];x[1]=xq[1];
  for(i=0;i<N;i++)
    for(j=0;j<N+1;j++)
      d[i][j]=0;
  for(i=0;i<N;i++)  d[i][i+1]=1;

  loop:
  #ifdef XuYimin_Debug_1
    printf("fangxiang juzhen:d=\n");
    for(i=0;i<N;i++)
      for(j=0;j<N+1;j++)
    {
      printf("%15g",d[i][j]);
      if(j==N)  printf("\n");
    }
    printf("\n");
  #endif

    for(i=0;i<N;i++)  X[i][0]=x[i];
    ALF[0]=0;del[0]=0;
    fu[0]=fun(x);
    for(id=1;id<N+1;id++)
      {
    for(i=0;i<N;i++)  tem[i]=d[i][id];
    ALF[id]=alf(x,tem);
    for(i=0;i<N;i++)  te[i]=x[i]+ALF[id]*d[i][id];
    fu[id]=fun(te);
    for(i=0;i<N;i++)  X[i][id]=te[i];
    for(i=0;i<N;i++)  x[i]=te[i];
      }
    for(id=1;id<N+1;id++)  del[id]=fu[id-1]-fu[id];
  #ifdef XuYimin_Debug_1
    printf("zuijia buchang juzhen:ALF:\n");
    for(i=0;i<N+1;i++)  printf("%15g",ALF[i]);
    printf("\n\ndiedaidian juzhen:X=\n");
    for(i=0;i<N;i++)
      for(id=0;id<N+1;id++)
        {
      printf("%15g",X[i][id]);
      if(id==N)  printf("\n");
        }
    printf("\nhanshuzhi juzhen:fu=\n");
    for(id=0;id<N+1;id++)  printf("%15g",fu[id]);
    printf("\n\nhanshuzhi chazhi juzhen:del=\n");
    for(id=0;id<N+1;id++)  printf("%15g",del[id]);
    printf("\n");
  #endif

  for(i=0;i<N;i++)
    {
      dk[i]=X[i][N]-X[i][0];
      xk[i]=2*X[i][N]-X[i][0];
    }
  for(i=0;i<N;i++)  x[i]=X[i][0];
  F[0]=fun(x);
  for(i=0;i<N;i++)  x[i]=X[i][N];
  F[1]=fun(x);
  F[2]=fun(xk);
  delm=maxt(del);
  for(i=0;i<N+1;i++)
    if(del[i]==delm)
      m=i;

  #ifdef XuYimin_Debug_1
    printf("x0_k=\n");
    for(i=0;i<N;i++)  printf("%15g\n",X[i][0]);
    printf("\nxn_k=\n");
    for(i=0;i<N;i++)  printf("%15g\n",X[i][N]);
    printf("\nxn+1_k=\n");
    for(i=0;i<N;i++)  printf("%15g\n",xk[i]);
    printf("\nhanshuzhi zengliang juzhen F=\n");
    for(i=0;i<3;i++)  printf("%15g",F[i]);
    printf("\n\nzuida chazhi:m=%d    delm=%g\n\n",m,delm);

  #endif
  pd[0]=(F[2]<F[0]);

  pd[1]=((F[0]-2*F[1]+F[2])*(F[0]-F[1]-delm)*(F[0]-F[1]-delm)
            <0.5*delm*(F[0]-F[2])*(F[0]-F[2]));
  pd[2]=(F[1]<F[2]);
  #ifdef XuYimin_Debug_1
    printf("pd[0]=%d   pd[1]=%d   pd[2]=%d\n\n",pd[0],pd[1],pd[2]);
  #endif

  if(pd[0]*pd[1]==1)
    {
      for(i=0;i<N;i++)  te[i]=X[i][N];
      al=alf(te,dk);
      for(i=0;i<N;i++)  x[i]=X[i][N]+al*dk[i];
      for(j=m;j<N;j++)
        for(i=0;i<N;i++)
          d[i][j]=d[i][j+1];
      for(i=0;i<N;i++)  d[i][N]=dk[i];
      #ifdef XuYimin_Debug_1
    printf("xin diedaidian:x=\n");
        for(i=0;i<N;i++)  printf("%15g\n",x[i]);
        printf("\nxin fangxiang juzhen:d=\n");
        for(i=0;i<N;i++)
          for(j=0;j<N+1;j++)
        {
          printf("%15g",d[i][j]);
          if(j==N)  printf("\n");
        }
         printf("\n");
      #endif
    }
  else
    if(pd[2]==1)
      for(i=0;i<N;i++)  x[i]=X[i][N];
    else
      for(i=0;i<N;i++)  x[i]=xk[i];

  pa[N]=0;
  for(i=0;i<N;i++)  pa[i]=X[i][N]-X[i][0];
  for(i=0;i<N;i++)  pa[N]+=pa[i]*pa[i];
  pa[N]=sqrt(pa[N]);
  #ifdef XuYimin_Debug_1
    printf("piancha juzhen:PA=\n");
    for(i=0;i<N;i++)
      printf("%15g\n",pa[i]);
    printf("\npiancha juli:pa=%g\n\n",pa[N]);
    printf("************Di %d ci diedai jieshu!************\n\n\n\n\n",xym++);
  #endif

  if(pa[N]<esp)
    {
      printf("The Result:\n");
      printf("jizhi dian:X=\n");
      for(i=0;i<N;i++)
    printf("%15g\n",x[i]);
      printf("\njizhidian:f(X)=%f",fun(x));
    }
  else goto loop;
}

int main()
{
  powell();
  getch();
  return 0;
}


程序刚编完，正在修改，希望大家提点意见！