#include "hjfgf.c"
double oneoptim(double x0[],double s[],double h0,double epsg,int n,double x[])
{double *a,*b,ff;
 a=(double *)malloc(n*sizeof(double));
 b=(double *)malloc(n*sizeof(double));
 jtf(x0,h0,s,n,a,b);
 ff=gold(a,b,epsg,n,x);
 free(a);
 free(b);
 return (ff);
}
double powell(double p[],double h0,double eps,double epsg,int n,double x[])
{int i,j,m;
 double *xx[4],*ss,*s;
 double f,f0,f1,f2,f3,fx,dlt,df,sdx,q,d;
 ss=(double *)malloc(n*(n+1)*sizeof(double));
 s=(double *)malloc(n*sizeof(double));
 for(i=0;i<n;i++)
  {for(j=0;j<=n;j++)
   *(ss+i*(n+1)+j)=0;
   *(ss+i*(n+1)+i)=1;
  }
  for(i=0;i<4;i++)
  xx[i]=(double *)malloc(n*sizeof(double));
  for(i=0;i<n;i++)
  *(xx[0]+i)=p[i];
  for(;;)
  {for(i=0;i<n;i++)
    {*(xx[1]+i)=*(xx[0]+i);
     x[i]=*(xx[1]+i);
    }
   f0=f1=objf(x);
   dlt=-1;
   for(j=0;j<n;j++)
    {for(i=0;i<n;i++)
      {*(xx[0]+i)=x[i];
       *(s+i)=*(ss+i*(n+1)+j);
      }
     f=oneoptim(xx[0],s,h0,epsg,n,x);
     df=f0-f;
     if(df>dlt)
      {dlt=df;
       m=j;
      }
    }
    sdx=0;
    for(i=0;i<n;i++)
     sdx=sdx+fabs(x[i]-(*(xx[1]+i)));
    if(sdx<eps)
     {free(ss);
      free(s);
      for(i=0;i<4;i++)
      free(xx[i]);
      return(f);
     }
    for(i=0;i<n;i++)
     *(xx[2]+i)=x[i];
    f2=f;
    for(i=0;i<n;i++)
    {*(xx[3]+i)=2*(*(xx[2]+i)-(*(xx[1]+i)));
     x[i]=*(xx[3]+i);
    }
    fx=objf(x);
    f3=fx;
    q=(f1-2*f2+f3)*(f1-f2-dlt)*(f1-f2-dlt);
    d=0.5*dlt*(f1-f3)*(f1-f3);
    if((f3<f1)||(q<d))
     {if(f2<=f3)
       for(i=0;i<n;i++)
        *(xx[0]+i)=*(xx[2]+i);
      else
       for(i=0;i<n;i++)
        *(xx[0]+i)=*(xx[3]+i);
     }
    else
    {for(i=0;i<n;i++)
     {*(ss+(i+1)*(n+1))=x[i]-(*(xx[1]+i));
      *(s+i)=*(ss+(i+1)*(n+1));
     }
   f=oneoptim(xx[0],s,h0,epsg,n,x);
   for(i=0;i<n;i++)
    *(xx[0]+i)=x[i];
   for(j=m+1;j<=n;j++)
   for(i=0;i<n;i++)
    *(ss+i*(n+1)+j-1)=*(ss+i*(n+1)+j);
    }
  }
}