#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#define ep 0.0001
#define e 0.01
double f(double *p)
{
    double y;
    y=25/(p[0]*pow(p[1],3));
    return y;
}
int cons(double *q)
{
    int n;
    if((30/(q[0]*pow(q[1],3))-50<=0)&&(0.0004*q[0]*q[1]-0.001<=0)&&(q[0]>=2)&&(q[0]<=-4)&&(q[1]>=0.5)&&(q[1]<=1))
    n=1;
    else
    n=0;
    return n;
}
void order(double *p1,double *p2,double *p3,double *p4)
{
    double f1,f2,f3,f4,h[2];
    int i;
    f1=f(p1);
    f2=f(p2);
    f3=f(p3);
    f4=f(p4);
    if(f1<f2)
    for(i=0;i<=1;i++)
    {
        h[i]=p1[i];
        p1[i]=p2[i];
        p2[i]=h[i];
    }
    if(f1<f3)
    for(i=0;i<=1;i++)
    {
        h[i]=p1[i];
        p1[i]=p3[i];
        p3[i]=h[i];
    }
    if(f1<f4)
    for(i=0;i<=1;i++)
    {
        h[i]=p1[i];
        p1[i]=p4[i];
        p4[i]=h[i];
    }
    if(f2<f3)
    for(i=0;i<=1;i++)
    {
        h[i]=p2[i];
        p2[i]=p3[i];
        p3[i]=h[i];
    }
    if(f2<f4)
    for(i=0;i<=1;i++)
    {
        h[i]=p2[i];
        p2[i]=p4[i];
        p4[i]=h[i];
    }
    if(f3<f4)
    for(i=0;i<=1;i++)
    {
        h[i]=p3[i];
        p3[i]=p4[i];
        p4[i]=h[i];
    }
}
    void main()
{
    int i,j=0,tf1,tf2;
    double x1[2]={2.2,0.55},x2[2]={2.4,0.6},x3[2]={2.6,0.65},x4[2]={2.8,0.7},xc[2],xr[2],a[2],b[2],h=1.3;
    loop1: order(x1,x2,x3,x4);
    if((pow((f(x1)-f(x4)),2)+pow((f(x2)-f(x4)),2)+pow((f(x3)-f(x4)),2))/3<ep)
    {
        printf("  \n\n\n  输出最优解及目标函数值：\n");
        printf("\n  x1=%.5f   x2=%.5f    f(x1,x2)=%.5f\n\n ",x4[0],x4[1],f(x4));
    }
    else
    {
        loop2: for(i=0;i<=1;i++)
        xc[i]=(x2[i]+x3[i]+x4[i])/3;
        tf1=cons(xc);
        if(tf1==0)
        {
          for(i=0;i<=1;i++)
          {
            a[i]=x4[i];
            b[i]=xc[i];
          }
          for(i=0;i<=1;i++)
          {
            x1[i]=a[i]+rand()/RAND_MAX*(b[i]-a[i]);
            x2[i]=a[i]+rand()/RAND_MAX*(b[i]-a[i]);
            x3[i]=a[i]+rand()/RAND_MAX*(b[i]-a[i]);
            x4[i]=a[i]+rand()/RAND_MAX*(b[i]-a[i]);  
          }
          goto loop1;
        }
        else   
        loop3: for(i=0;i<=1;i++)
        xr[i]=xc[i]+h*(xc[i]-x1[i]);
        tf2=cons(xr);
        if(tf2==0)
        {
          loop4: h=h/2;
          goto loop3;
        }
        else
         if(f(xr)<f(x1))
         {
           for(i=0;i<=1;i++)
           x1[i]=xr[i];
           j++;
           printf("第%d次迭代\n",j);
           printf("xr1=%f   xr2=%f\n",xr[0],xr[i]);
           goto loop1;
         }
         else
           if(h<=pow(10,-10))
           {
             for(i=0;i<=1;i++)
             x1[i]=x2[i];
             goto loop2;
           }
           else
           goto loop4;
    }
}