golden.c

golden method in c, please enjoy it

/* 搏垛尸充恕≤跋い哈み恕 */
/* Coded by T.Takahama, 1998 */

#define SQRT5	2.2360679
#define F1	(3.0-SQRT5)/2.0
#define F2	(SQRT5-1.0)/2.0
#define EPSILON	1.0e-3

int NEval=0;

double evaluate(double x)
{
	++NEval;
	return((x-1.0)*(x-1.0));
}

main()
{
	double x[3], y[3], x_min, y_min;

	bracketing(evaluate, 1.0, 1.0, x, y);
	golden(evaluate, x[0], x[2], EPSILON, &x_min, &y_min);
	printf("f(%lf)=%lf\n", x_min, y_min);
	printf("%d evaluations\n", NEval);
}

/* golden(f, x_min, x_max, epsilon, x, y)
  簇眶 f を惰粗 x_min から x_max の粗で玫瑚し、呵井猛を艰る x と y を手す
  惰粗の升が epsilon 笆布の箕匿贿する
*/
golden(double (*f)(), double x_min, double x_max, double epsilon,
	double *x, double *y)
{
	double x1, y1, x2, y2;

	x1=x_min+F1*(x_max-x_min); y1=(*f)(x1);
	x2=x_min+F2*(x_max-x_min); y2=(*f)(x2);
	while(x_max-x_min>epsilon) {
		if(y1<y2) {
			x_max=x2;
			x2=x1; y2=y1;
			x1=x_min+F1*(x_max-x_min); y1=(*f)(x1);
		}
		else {
			x_min=x1;
			x1=x2; y1=y2;
			x2=x_min+F2*(x_max-x_min); y2=(*f)(x2);
		}
	}
	if(y1<=y2) {
		*x=x1; *y=y1;
	}
	else {
		*x=x2; *y=y2;
	}
}

/*
    bracketing(f, x0, h, 手猛 x[3], 手猛 y[3])
	簇眶 f に簇して介袋爬 x0 から 癸み h で跋い哈みを悸乖し、
	y[0]>=y[1]<=y[2] (霹规は办数のみ) となる x[0]<x[1]<x[2] を手す
*/
bracketing(double (*f)(), double x0, double h, double x[], double y[])
{
	double y0, x1, y1, x2, y2;

	y0=(*f)(x0);
	x1=x0+h;
	y1=(*f)(x1);
	if(y1>y0) {
		y2=y1; y1=y0; y0=y2; /* swap y0 y1 */
		x2=x1; x1=x0; x0=x2; /* swap x0 x1 */
		h=-h;
	}
	for(;;) {
		h*=2.0;
		x2=x1+h;
		y2=(*f)(x2);
		if(y2>y1) break;
		y0=y1; x0=x1;
		y1=y2; x1=x2;
	}
	if(h>=0.0) {
		x[0]=x0; x[1]=x1; x[2]=x2;
		y[0]=y0; y[1]=y1; y[2]=y2;
	}
	else {
		x[0]=x2; x[1]=x1; x[2]=x0;
		y[0]=y2; y[1]=y1; y[2]=y0;
	}
}