Bisection to find roots in C++ I have written a short C/C++ code finding root by bisection. (This is a simple iterative numerical method allowing to find the root of an equation i.e. x such that f(x) = 0). Bisection Method

The header simply consists of guards and of the following lines: 

#include &lt;functional&gt;
double bisection(double x1, 
                 double x2, 
                 double e, 
                 std::function&lt;double(double)&gt;&amp; f);


In order to get plain C code, the std::function&lt;double(double)&gt; in the signature of the function has to be turned into a double(*f)(double);

I have two questions. 

The first is if this code can be improved (minor improvements). 
The second is if recursion is a good choice. 

namespace { 
    double bisection(double x1, 
                     double x2, 
                     double e, 
                     std::function&lt;double(double)&gt;&amp; f, 
                     double fx1, 
                     double fx2){
        double mid = x1*0.5 + x2*0.5;
        if (x2 - x1 &lt; e)
            return mid;
        double fmid = f(mid);
        if ((fmid&gt;0.0) == (fx1 &gt; 0.0))
            return bisection(mid, x2, e, f, fmid, fx2);
        return bisection(x1, mid, e, f, fx1, fmid);
    }
}

double bisection(double x1, 
                 double x2, 
                 double e, 
                 std::function&lt;double(double)&gt;&amp; f){
    double fx1 = f(x1);
    double fx2 = f(x2);
    return bisection(x1, x2, e, f, fx1, fx2);
}


An example of main using my code is: 

#include "bisection.hpp"
int main(){
    bisection(0.0, 1.0, 1.0E-6, [](double x){return x*x-0.5;})

    return 0;
}       


please notice that x1 and x2 are expected to be ordered: x1 &lt; x2. Moreover the user should be aware of the fact that there will be a log_2((x_2-x_1)/tolerance) calls stack.
If any of fmid, fx1 or fx2 are zero (or close to zero), you should stop, since you already found an answer. That's also in Wikipedia's example algorithm. Even for the constant function f = [](double){return 0;}, your code will still continue and take a total of \$ \log_2 \left(\frac{x_2-x_1}{\varepsilon}\right)\$ calls.
Your recursion cannot get transformed into a loop by GCC. The recursive algorithm looks nice, but since we end up with \$ \log_2 \left(\frac{x_2-x_1}{\varepsilon}\right)\$ calls. An iterative variant is easy to write, so I'd use that one.
Not only must \$x_1 &lt;x_2 \$, but also \$f(x_1) f(x_2) &lt; 0\$, that is, the signs of \$f(x_1)\$ and \$f(x_2)\$ must differ. Otherwise your algorithm might not find the root. I only mention this because you wrote that \$x_1\$ and \$x_2\$ are expected to be ordered, but that's not enough.


Other than that, I strongly suggest you to use braces even for single statement ifs, as missing braces have been a source of errors too often. 

Alternatively a little bit more whitespace in your function can make it a lot easier to read:

namespace { 
double bisection(double x1, 
                 double x2, 
                 double e, 
                 std::function&lt;double(double)&gt;&amp; f, 
                 double fx1, 
                 double fx2){

    const double mid = x1*0.5 + x2*0.5;
    const double fmid = f(mid);

    if (x2 - x1 &lt; e || fmid == 0.0) {
        return mid;
    }


    if ((fmid &gt; 0.0) == (fx1 &gt; 0.0)) {
        return bisection(mid, x2, e, f, fmid, fx2);
    } else {
        return bisection(x1, mid, e, f, fx1, fmid);
    }
}
}


I would use the ternary operator though instead of if/else, but that's only a style question. But then again, I would write this in an iterative way, since we don't need an additional helper.