return std::make_pair(b, root_is_valid);
return b;
double brent_zero    (double          a    ,double          b    ,double          t    ,FunctionalType& f    ){    double const epsilon = std::numeric_limits<double>::epsilon();    // Returns a zero of the function f in [a,b] to within a tolerance    //   6 * epsilon * |z| + 2 * t    // f(a) and f(b) must have different signs.    double fa = f(a);    double fb = f(b);    double fc = fb;    for(;;)        {        double c, d, e;        if((0.0 < fb) == (0.0 < fc))            {            c = a;            fc = fa;            d = e = b - a;            }        if(std::fabs(fc) < std::fabs(fb))            {             a =  b;  b =  c;  c =  a;            fa = fb; fb = fc; fc = fa;            }        double tol = 2.0 * epsilon * std::fabs(b) + t;        double m = 0.5 * (c - b);        if(0.0 == fb || std::fabs(m) <= tol)            {            return std::make_pair(b, root_is_valid);            }        if(std::fabs(e) < tol || std::fabs(fa) <= std::fabs(fb))            {            // Bisection.            d = e = m;            }        else            {            double p, q;            double s = fb / fa;            if(a == c)                {                // Linear interpolation.                p = 2.0 * m * s;                q = 1.0 - s;                }            else                {                // Inverse quadratic interpolation.                q = fa / fc;                double r = fb / fc;                p = s * (2.0 * m * q * (q - r) - (b - a) * (r - 1.0));                q = (q - 1.0) * (r - 1.0) * (s - 1.0);                }            if(0.0 < p)                {                q = -q;                }            else                {                p = -p;                }            s = e;            e = d;            if                (   p < 1.5 * m * q - std::fabs(tol * q)                &&  p < std::fabs(0.5 * s * q)                )                {                d = p / q;                }            else                {                d = e = m;                }            }        a = b;        fa = fb;        if(tol < std::fabs(d))            {            b += d;            }        else if(0.0 < m)            {            b += tol;            }        else            {            b -= tol;            }        fb = f(b);        }}