/*
************************************************************************
*
* Numerical Math Package
*
* Brent's root finder
* obtains a zero of a function of one variable
*
* Synopsis
* double zeroin(ax,bx,f,tol=EPSILON)
* const double ax The root is to be sought within
* const double bx the interval [ax,bx]
* UnivariateFunctor& f The function under consideration
* const double tol Acceptable tolerance for the root
* position. It is an optional parameter
* with default value DBL_EPSILON
*
* Zeroin returns an approximate location for the root with accuracy
* 4*DBL_EPSILON*abs(x) + tol
*
* Algorithm
* G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
* computations. M., Mir, 1980, p.180 of the Russian edition
*
* The function makes use of a bissection procedure combined with
* a linear or quadratic inverse interpolation.
* At each step the code operates three abscissae - a, b, and c:
* b - the last and the best approximation to the root
* a - the last but one approximation
* c - the last but one or even an earlier approximation such that
* 1) |f(b)| <= |f(c)|
* 2) f(b) and f(c) have opposite signs, i.e. b and c encompass
* the root
* Given these abscissae, the code computes two new approximations, one by the
* bissection procedure and the other one from interpolation (if a,b, and c
* are all different the quadratic interpolation is used, linear otherwise).
* If the approximation obtained by the interpolation looks
* reasonable (i.e. falls within the current interval [b,c], not too close
* to the end points of the interval), the point is accepted as a new
* approximation to the root. Otherwise, the result of the bissection is used.
* Therefore, the range of uncertainty is guaranteed to tighten at
* least by a factor of 1.6
*
*
************************************************************************
*/
#define DBL_EPSILON 2.22045e-16
#include "assert.h"
template
inline double zeroin( // An estimate to the root
const double ax, // Specify the interval the root
const double bx, // to be sought in
const double cx, // root between a and c, best approximation b
functor& f, // Function under investigation
const double tol) // Acceptable tolerance
{
Assert( tol > 0, "Tolerance must be positive");
Assert( cx > ax,
"Left end point of the interval should be strictly less than the "
"right one" );
double b = bx; // the last and the best approx to the root
double fb = f(b);
double a = ax; // the last but one approximation
double fa = f(a);
double c = cx; // the last but one or even an earlier approx
double fc = f(c);
if (fc*fb>0){
// swap a and c such that the root is between b and c
double t = a, ft = fa;
a = c; fa = fc;
c = t; fc = ft;
}
for(;;) // Main iteration loop
{
const double prev_step = b-a; // Step from the previous iteration
if( abs(fc) < abs(fb) )
{ // Swap data so that b would be the
a = b; b = c; c = a; // best approximation found so far
fa=fb; fb=fc; fc=fa;
}
// Estimate the effective tolerance
const double tol_act = 2*DBL_EPSILON*abs(b) + tol/2;
double new_step = (c-b)/2; // Bissection step for this iteration
if( abs(new_step) <= tol_act || fb == 0 )
return b; // Acceptable approximation is found
// Figuring out if the interpolation can be tried
if( abs(prev_step) >= tol_act // If prev_step was large enough
&& abs(fa) > abs(fb) ) // and was in true direction,
{ // Interpolatiom may be tried
double p; // Interpolation step is calcu-
double q; // lated in the form p/q; divi-
// sion operations is delayed
// until the last moment
const double cb = c-b;
if( a==c ) // If we've got only two distinct
{ // points linear interpolation
register const double t1 = fb/fa; // can only be applied
p = cb*t1;
q = 1.0 - t1;
}
else // Quadratic inverse interpolation
{
register const double t1=fb/fc, t2=fb/fa;
q = fa/fc;
p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
q = (q-1.0) * (t1-1.0) * (t2-1.0);
}
if( p > 0 ) // Formulas above computed new_step
q = -q; // = p/q with the wrong sign (on purpose).
else // Correct this, but in such a way so
p = -p; // that p would be positive
if( 2*p < (1.5*cb*q-abs(tol_act*q)) // If b+p/q falls in [b,c]
&& 2*p < abs(prev_step*q) ) // and isn't too large
new_step = p/q; // it is accepted
// If p/q is too large then the
// bissection procedure can
// reduce [b,c] to a larger
// extent
}
if( abs(new_step) < tol_act ) // Adjust the step to be not less
new_step = new_step > 0 ? // than the tolerance
tol_act : -tol_act;
a = b; fa = fb; // Save the previous approximation
b += new_step; fb = f(b); // Do step to a new approximation
if( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) )
{ // Adjust c for it to have the sign
c = a; fc = fa; // opposite to that of b
}
}
}