5
\$\begingroup\$

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 <functional>
double bisection(double x1, 
                 double x2, 
                 double e, 
                 std::function<double(double)>& f);

In order to get plain C code, the std::function<double(double)> 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<double(double)>& f, 
                     double fx1, 
                     double fx2){
        double mid = x1*0.5 + x2*0.5;
        if (x2 - x1 < e)
            return mid;
        double fmid = f(mid);
        if ((fmid>0.0) == (fx1 > 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<double(double)>& 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 < x2. Moreover the user should be aware of the fact that there will be a log_2((x_2-x_1)/tolerance) calls stack.

\$\endgroup\$
  • \$\begingroup\$ Although this is probably well known, could you please provide us usage example? Small .cpp file with all includes and main() would be great. Usually include directives might be reviewed too. \$\endgroup\$ – Incomputable Jun 6 '18 at 12:51
4
\$\begingroup\$
  1. 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.
  2. 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.
  3. Not only must \$x_1 <x_2 \$, but also \$f(x_1) f(x_2) < 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<double(double)>& f, 
                 double fx1, 
                 double fx2){

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

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


    if ((fmid > 0.0) == (fx1 > 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.

\$\endgroup\$
  • \$\begingroup\$ I was aware of point 3, I didn't list it because this is a "bug" of the method itself while the fact that they are expected ordered is an implementation detail. In any case thanks for pointing it out. And thanks for points 1 and 2 too. \$\endgroup\$ – jimifiki Jun 6 '18 at 14:28
3
\$\begingroup\$
  • Temporary can't bind to a non-const reference in standard C++ (MSVC++ will compile it because of compiler extension): make f a const reference
  • Const correctness: use const wherever possible
  • Some temporary variables can be avoided. E.g., fmid
  • Is there a reason you do x1 * 0.5 + x2 * 0.5 instead of (x1 + x2) * 0.5? If for some numerical reasons, it's fine. Otherwise, I find the latter more intuitive
  • You don't really need fx2
  • Recursion can indeed be more both more costly (if not optimized by the compiler) and more difficult to read.

Given all of the above considerations, here's a simple non-recursive version:

double bisection(
    double x1,
    double x2,
    const double e,
    const std::function<double(double)>& f)
{
    while(x2 - x1 >= e)
    {
        const double mid = (x1 + x2) * 0.5;
        if((f(mid) > 0.0) == (f(x1) > 0.0))
            x1 = mid;
        else
            x2 = mid;
    }
    return (x1 + x2) * 0.5;
}

And here's the version that caches f(x1):

double bisection(
    double x1,
    double x2,
    const double e,
    const std::function<double(double)>& f)
{
    double fx1 = f(x1);
    while(x2 - x1 >= e)
    {
        const double mid = (x1 + x2) * 0.5;
        const double fmid = f(mid);
        if((fmid > 0.0) == (fx1 > 0.0))
        {
            x1 = mid;
            fx1 = fmid;
        }
        else
        {
            x2 = mid;
        }
    }
    return (x1 + x2) * 0.5;
}
\$\endgroup\$
  • \$\begingroup\$ I really appreciate your effort in producing a version caching f(x1). Offering this version is of capital importance IMO. \$\endgroup\$ – jimifiki Jun 6 '18 at 14:55
  • \$\begingroup\$ Would you consider if(fmid==0) break; ? \$\endgroup\$ – jimifiki Jun 6 '18 at 14:57
  • 1
    \$\begingroup\$ About the strange x1 * 0.5 + x2 * 0.5 I've used it in order to avoid any risk of overflow. I am not sure this is of help in 99.99999 percent of use cases. \$\endgroup\$ – jimifiki Jun 6 '18 at 15:01
  • \$\begingroup\$ @jimifiki Ah, I see. Thanks for the explanation. As for the if(fmid==0), it depends on the usage patterns. IMO it's not worth to have an extra comparison for this very rare case. In other words it's a trade-off: optimized function in 0.001% vs. slightly slower in 99.999%. \$\endgroup\$ – AMA Jun 6 '18 at 15:08
1
\$\begingroup\$

In addition to the suggestions from other answers, idiomatic C++ would use a template parameter for the function object type, rather than std::function. std::function is great when you need late binding of functions at runtime for whatever reason, but it's also a tremendous optimization barrier. I suspect this will be much faster:

template <typename Func>
double bisection(double x1, double x2, double e, Func f)
{
    double fx1 = f(x1);
    double mid = x1 * 0.5 + x2 * 0.5;
    while (x2 - x1 >= e) {
        double fmid = f(mid);
        if ((fx1 > 0.0) == (fmid > 0.0)) {
            x1 = mid;
            fx1 = fmid;
        } else {
            x2 = mid;
        }
        mid = x1 * 0.5 + x2 * 0.5;
    }
    return mid;
}

(Credit to AMA's answer for most of that implementation.)

You could even make it slightly more general by parameterizing the value type instead of assuming double.

\$\endgroup\$
  • 1
    \$\begingroup\$ You would also use an Func && f usually. However, std::function<…>::operator()(…) can get optimized by the compiler in some circumstances, so it's not "much faster" in general. \$\endgroup\$ – Zeta Jun 7 '18 at 7:22
  • \$\begingroup\$ A naive question: why is std::function<…>::operator()(…) likely to be slower than a pointer to function? I would be grateful to anyone sharing a link to any relevant content. \$\endgroup\$ – jimifiki Jun 8 '18 at 6:04
  • \$\begingroup\$ My understanding is that commonstd::function implementations will need two indirections, not one: first to handle the various types of function that can be wrapped such as function objects, then another to actually call the pointed-to function. But even one function pointer is much slower than zero. \$\endgroup\$ – Tavian Barnes Jun 8 '18 at 14:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.