/* 以下為復合形法法子程序 */
#include "stdlib.h"
#include "math.h"
#include "stdio.h"
float objfx(float x[]);
void constraint(float x[],float g[]);

int gau(float x[],float g[],int kg)
{
       int i;
       constraint(x,g);
       for(i=0;i<kg;i++)
       {
         if(g[i]<0)
         goto s333;
       }
       return 1;
s333:return 0;
}

void xcent(int n,int k,int ll,int lh,float x0[],float xcom[][100])
{
       int i,l;
       float xs;
       for(i=0;i<n;i++)
       {
         xs=0;
         for(l=0;l<ll;l++)
         {
    if(l!=lh)
    xs=xs+xcom[i][l];
         }
         if(lh>-1)
    x0[i]=xs/(ll-1);
         else
    x0[i]=xs/ll;
       }
}

void fxse(int n,int k,float x[],float xcom[][100],float fxk[])
{
       int l,lp,lp1,i;
       float w;
       for(l=0;l<k-1;l++)
       for(lp=0;lp<k-l;lp++)
       {
         lp1=lp+1;
         if(fxk[lp]<=fxk[lp1])
         {
    w=fxk[lp];
    fxk[lp]=fxk[lp1];
    fxk[lp1]=w;
    for(i=0;i<n;i++)
    {
      x[i]=xcom[i][lp];
      xcom[i][lp]=xcom[i][lp1];
      xcom[i][lp1]=x[i];
    }
         }
       }
}


void complex(int n,int k,int kg,float ep,float x[],float bl[],float bu[],
        float xcom[][100],float *f,int *nf,int *ng)
{
       int i,iw,l,ll,lh,it;
       float fx,fx0,sdx,fxh,fxr,alp;
       float *x0=(float*)calloc(n,sizeof(float));
       float *xh=(float*)calloc(n,sizeof(float));
       float *xr=(float*)calloc(n,sizeof(float));
       float *fxk=(float*)calloc(k,sizeof(float));
       float *g=(float*)calloc(kg,sizeof(float));
s5:    for(i=0;i<n;i++)
       x[i]=bl[i]+rand()/40000.0*(bu[i]-bl[i]);
       iw=gau(x,g,kg);
       *ng=*ng+1;
       if(iw==0)goto s5;
       for(i=0;i<n;i++)
       xcom[i][0]=x[i];
       for(l=1;l<k;l++)
       for(i=0;i<n;i++)
       xcom[i][l]=bl[i]+rand()/50000.0*(bu[i]-bl[i]);
       lh=-1;
       for(ll=1;ll<k;ll++)
       {
         xcent(n,k,ll,lh,x0,xcom);
         iw=gau(x0,g,kg);
         *ng=*ng+1;
         if(iw==0)goto s5;
         for(i=0;i<n;i++)
         x[i]=xcom[i][ll+1];
s24:     iw=gau(x,g,kg);
         *ng=*ng+1;
         if(iw==0)
         {
    for(i=0;i<n;i++)
    x[i]=x0[i]+0.5*(x[i]-x0[i]);
    goto s24;
         }
         else
         {
    for(i=0;i<n;i++)
    xcom[i][ll+1]=x[i];
         }
       }
       for(l=0;l<k;l++)
       {
         for(i=0;i<n;i++)
         x[i]=xcom[i][l];
         fx=objfx(x);
         *nf=*nf+1;
         fxk[l]=fx;
       }
       it=0;
s14: it=it+1;
       printf("\n\n    +++ ITERATION COMPUTE +++");
       printf("\n       ITER    = %d",it);
       lh=-1;
       xcent(n,k,k,lh,x0,xcom);
       fx0=objfx(x0);
       *nf=*nf+1;
       iw=gau(x0,g,kg);
       *ng=*ng+1;
       printf("\n       Fmid    = %f",fx0);
       for(i=0;i<n;i++)
         printf("\n       X(%d)mid=%f",i,x0[i]);
       for(i=0;i<kg;i++)
         printf("\n       G(%d)mid=%f",i,g[i]);
       sdx=0;
       for(l=0;l<k;l++)
         sdx=sdx+(fx0-fxk[l])*(fx0-fxk[l]);
       sdx=sqrt(sdx/(float)k);
       if(sdx<ep)goto s38;
       fxse(n,k,x,xcom,fxk);
       lh=0;
s22: fxh=fxk[lh];
       for(i=0;i<n;i++)
         xh[i]=xcom[i][lh];
       xcent(n,k,k,lh,x0,xcom);
       iw=gau(x0,g,kg);
       *ng=*ng+1;
       if(iw==0)goto s36;
       alp=1.3;
s12: for(i=0;i<n;i++)
       xr[i]=x0[i]+alp*(x0[i]-xh[i]);
       iw=gau(xr,g,kg);
       *ng=*ng+1;
       if(iw==0)
       {
         alp=alp*0.5;
         goto s12;
       }
       fxr=objfx(xr);
       *nf=*nf+1;
       if(fxr>=fxh)
       {
         if(alp>1.0e-4)
         {
    alp=alp*0.5;
    goto s12;
         }
         lh=lh+1;
         if(lh<3)goto s22;
       }
       for(i=0;i<n;i++)
       xcom[i][lh]=xr[i];
       fxk[lh]=fxr;
       goto s14;
s36: for(i=0;i<n;i++)
       {
         bl[i]=xcom[i][k];
         bu[i]=x0[i];
       }
       goto s5;
s38: for(i=0;i<n;i++)
       x[i]=x0[i];
       *f=objfx(x);
       *nf=*nf+1;
       free(x0);
       free(xh);
       free(xr);
       free(g);
       free(fxk);
}