File Metadata

brent.c
View Options

	#include <stdio.h>
	#include <math.h>

	#define ITMAX 100
	#define CGOLD 0.3819660
	#define ZEPS 1.0e-10
	#define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
	#define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);

	double brent(double ax, double bx, double cx, double (f)(double), double tol, double xmin)
	{
	int iter;
	double a, b, d, etemp, fu, fv, fw, fx, p, q, r, tol1, tol2, u, v, w, x, xm;
	double e = 0.0;

	d = -1e10;

	a = ((ax < cx) ? ax : cx);
	b = ((ax > cx) ? ax : cx);
	x = w = v = bx;
	fw = fv = fx = (*f)(x);
	for (iter = 0; iter < ITMAX; iter++)
	{
	xm = 0.5 * (a + b);
	tol2 = 2.0 * (tol1 = tol * fabs(x) + ZEPS);
	if (fabs(x - xm) <= (tol2 - 0.5*(b - a)))
	{
	*xmin = x;
	return fx;
	}
	if (fabs(e) > tol1)
	{
	r = (x - w) * (fx - fv);
	q = (x - v) * (fx - fw);
	p = (x - v) * q - (x - w) * r;
	q = 2.0 * (q - r);
	if (q > 0.0) p = -p;
	q = fabs(q);
	etemp = e;
	e = d;
	if (fabs(p) >= fabs(0.5qetemp) \|\| p <= q(a - x) \|\| p >= q(b - x))
	d = CGOLD * (e = (x >= xm ? a - x : b - x));
	else
	{
	d = p / q;
	u = x + d;
	if (u - a < tol2 \|\| b - u < tol2)
	d = SIGN(tol1, xm - x);
	}
	}
	else
	{
	d = CGOLD * (e = (x >= xm ? a - x : b - x));
	}
	u = (fabs(d) >= tol1 ? x + d : x + SIGN(tol1, d));
	fu = (*f)(u);
	if (fu <= fx)
	{
	if (u >= x) a = x;
	else b = x;
	SHFT(v, w, x, u)
	SHFT(fv, fw, fx, fu)
	}
	else
	{
	if (u < x) a = u;
	else b = u;
	if (fu <= fw \|\| w == x)
	{
	v = w;
	w = u;
	fv = fw;
	fw = fu;
	}
	else if (fu <= fv \|\| v == x \|\| v == w)
	{
	v = u;
	fv = fu;
	}
	}
	}
	fprintf(stderr, "Too many iterations in BRENT\n");
	*xmin = x;
	return fx;
	}

	#undef ITMAX
	#undef CGOLD
	#undef ZEPS
	#undef SIGN