//本實驗求函數 f(x)=(x-3)2 的搜索區間[a,b]。并用黃金分割法或插值法求最優解
#include <math.h>
#include <stdio.h>
/*  函數 f(x)=(x-3)2 	*/
double f(double x) 
{  return (x-3)*(x-3);
}
/* 求搜索區間[a,b]的函數, x0---初始點; h0---初始步長 */
void find_ab(double x0,double h0,double *a,double *b)
{
/* 參照流程圖自己寫出外推法（進退法）程序的內容 */

}

/*
  黃金分割法
  a,b---搜索區間[a,b]; e---精度
  x,y---最優解X*,F*
*/
void search_gold(double a,double b,double e,
			  double *x,double *y)
{
	double x1,x2,y1,y2;
	x1=a+0.382*(b-a); y1=f(x1);
	x2=a+0.618*(b-a); y2=f(x2);
	do {
	  if (y1<y2) {
		 b=x2; x2=x1; y2=y1;
		 x1=a+0.382*(b-a); y1=f(x1);
	  } else {
		 a=x1; x1=x2; y1=y2;
		 x2=a+0.618*(b-a); y2=f(x2);
	  }
	} while (b-a>e);
	*x=0.5*(a+b); *y=f(*x);
}
/* 本實驗用鮑威爾方法求函數 f(x)=4(x1-5)2+(x2-6)2 的最優解*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#define MAXN 10
double xkk[MAXN],xk[MAXN],sk[MAXN];
int N;

double F(double *x) 
{
 return 4*pow(x[0]-5,2)+pow(x[1]-6,2);
}

double f(double x)
{   int i; 
	for (i=0; i<N; i++)
	{ xkk[i]=xk[i]+x*sk[i];
	}
	return f(xkk);
}
/*
	無約束坐標輪換法
	x0--初始點
	e1--一維搜索精度
	e2--求解精度
*/
double nc_trans(double *x0,double e1,double e2)
{
	int i,j,k=1;
	double a,b,ax,ay,d;
	for (;;) {
	  for (j=0; j<N; j++) xk[j]=x0[j];
	  for (i=0; i<N; i++) 
	  {
		 for (j=0; j<N; j++) if (j==i) sk[j]=1; 
         else sk[j]=0;
		 find_ab(0,1,&a,&b);
		 search_gold(a,b,e2,&ax,&ay);
		 for (j=0; j<N; j++) xk[j]=xkk[j];
	  }
	  d=0;
	  for (j=0; j<N; j++) 
      d+=(x0[j]-xkk[j])*(x0[j]-xkk[j]);
	  d=sqrt(d);
	  printf("k=%d;",k);
	  for (j=0; j<N; j++)
      printf(",x[%d]=%lf;",j+1,xkk[j]);
	  printf("d=%lf\n",d);
	  if (d<=e1) break;
	  for (j=0; j<N; j++) x0[j]=xkk[j];
	  k++;
	}
	for (j=0; j<N; j++) x0[j]=xkk[j];
	return F(xkk);
}