// http://en.wikipedia.org/wiki/Brent%27s_method (6/17/06)
double
Brent(double func(double, char *), double a, double b, double tol,
	       double ftol, char *params)
{
    const double eps = 1e-16;
    double x, c, d, fa, fb, fc, s, fs, tmp;
    fa = func(a, params);
    fb = func(b, params);
    if (fabs(fa)<fabs(fb)){
        tmp= a; a= b; b=tmp;
        tmp=fa;fa=fb;fb=tmp;
    }
    c =  a; fc = fa;
    bool mflag = true;
    while (fabs(a-b) > tol && (fabs(b-c) > tol || fabs(fb) > ftol)) {
//        cout << "c,b - " << c << ", " << b << endl;
        if (    fabs(fa-fc) > eps 
            &&  fabs(fb-fc) > eps 
            &&  fabs((fa-fc)*(fb-fc)) > eps) {
//            cout << "quadratic s ";
            s = a*fb*fc/((fa-fb)*(fa-fc))
              + b*fc*fa/((fb-fc)*(fb-fa))
              + c*fa*fb/((fc-fa)*(fc-fb));
        } else {
//            cout << "linear s ";
            s = b - fb*(b-a)/(fb-fa);
        }
        if ( fabs(s-(3.*a+5.*b)/8.) >  fabs(3./8.*(a-b))
            || ( mflag && fabs(s-b) >= fabs(b-c)/2)
            || (!mflag && fabs(s-b) >= fabs(c-d)/2)){
//            cout << "rejected" << endl;
            s = (a+b)/2;
            mflag = true;
        } else {
//            cout << "accepted" << endl;
            mflag = false;
        }
        fs = func(s, params);
        d = c; c = b; fc = fb;
        if (fa*fs < 0){
            b = s ; fb = fs;
        } else {
            a = s ; fa = fs;
        }
        if (fabs(fa) < fabs(fb)){
            tmp= a; a= b; b=tmp;
            tmp=fa;fa=fb;fb=tmp;
        }    
    }
    return b;
}