Implementation of Brent's Algorithm to find roots of a polynomial

Question

I made a program that contains a root-finding algorithm for polynomials as a function and contains 3 test polynomials. The algorithm is Brent's method and is based entirely off the pseudocode from Wikipedia. Is there anything I should change to make the code faster/easier to understand/etc?

The code does run and prints correct results for all test cases I've tried.

I'm mostly interested in writing optimized code with the priority being faster runtime, less memory usage, and fewer redundant function calls respectively, but I'm also interested in picking up better coding practices.

/*******************************************************************************
*
* Grant Williams
*
* Version 1.0.0
* August 27, 2015
*
* Test Program for Brent's Method Function.
*
* Brent's method makes use of the bisection method, the secant method, and inverse quadratic interpolation in one algorithm.
*
* To Compile Please use icc -std=c++11 if using intel or g++ -std=c++11 if using GCC.
*
********************************************************************************/

#include <iostream>
#include <cmath>
#include <algorithm>
#include <functional>
#include <ctime>
//#include "brent_fun.h"

void brents_fun(std::function<double (double)> f, double lower_bound, double upper_bound, double TOL, double MAX_ITER);

int main()
{
    //clock stuff
    std::clock_t start;
    double duration;
    start = std::clock();
    //stop clock stuff 

    double a;               // lower bound
    double b;               // upper bound
    double TOL = 0.0001;    // tolerance
    double MAX_ITER = 1000; // maximum number of iterations

    int function;           // set polynomial to find roots of & boundary conditions for that polynomial

    std::cout << std::endl;

    for (function = 1; function <= 3; function++)
    {
        if (function == 1)
        {
            a = -1.5;
            b = 0;
            auto f = [](double x){ return (x+1) * (x+2) * (x+3); };
            brents_fun(f,a,b,TOL,MAX_ITER);
        }
        else if (function == 2)
        {
            a = -10;
            b = 10;
            auto f = [](double x){ return (x*x*x -4*x - 9); };
            brents_fun(f,a,b,TOL,MAX_ITER);
        }
        else if (function == 3)
        {
            a = -4;
            b = 3;
            auto f = [](double x){ return (x+3)*(x-1)*(x-1); };
            brents_fun(f,a,b,TOL,MAX_ITER);
        }
    }




    //clock stuff again
    duration = (std::clock() - start) / (double) CLOCKS_PER_SEC;
    std::cout << "Elapsed time: " << duration << " seconds" << std::endl;
    //finish clock stuff

    std::cout << std::endl;

    return 0;
}

void brents_fun(std::function<double (double)> f, double lower_bound, double upper_bound, double TOL, double MAX_ITER)
{
    double a = lower_bound;
    double b = upper_bound;
    double fa = f(a);   // calculated now to save function calls
    double fb = f(b);   // calculated now to save function calls
    double fs = 0;      // initialize 

    if (!(fa * fb < 0))
    {
        std::cout << "Signs of f(lower_bound) and f(upper_bound) must be opposites" << std::endl; // throws exception if root isn't bracketed
        return;
    }

    if (std::abs(fa) < std::abs(b)) // if magnitude of f(lower_bound) is less than magnitude of f(upper_bound)
    {
        std::swap(a,b);
        std::swap(fa,fb);
    }

    double c = a;           // c now equals the largest magnitude of the lower and upper bounds
    double fc = fa;         // precompute function evalutation for point c by assigning it the same value as fa
    bool mflag = true;      // boolean flag used to evaluate if statement later on
    double s = 0;           // Our Root that will be returned
    double d = 0;           // Only used if mflag is unset (mflag == false)

    for (unsigned int iter = 1; iter < MAX_ITER; ++iter)
    {
        // stop if converged on root or error is less than tolerance
        if (std::abs(b-a) < TOL)
        {
            std::cout << "After " << iter << " iterations the root is: " << s << std::endl;
            return;
        } // end if

        if (fa != fc && fb != fc)
        {
            // use inverse quadratic interopolation
            s =   ( a * fb * fc / ((fa - fb) * (fa - fc)) )
                + ( b * fa * fc / ((fb - fa) * (fb - fc)) )
                + ( c * fa * fb / ((fc - fa) * (fc - fb)) );
        }
        else
        {
            // secant method
            s = b - fb * (b - a) / (fb - fa);
        }

        /*
            Crazy condition statement!:
            -------------------------------------------------------
            (condition 1) s is not between  (3a+b)/4  and b or
            (condition 2) (mflag is true and |s−b| ≥ |b−c|/2) or
            (condition 3) (mflag is false and |s−b| ≥ |c−d|/2) or
            (condition 4) (mflag is set and |b−c| < |TOL|) or
            (condition 5) (mflag is false and |c−d| < |TOL|)
        */
        if (    ( (s < (3 * a + b) * 0.25) || (s > b) ) ||
                ( mflag && (std::abs(s-b) >= (std::abs(b-c) * 0.5)) ) ||
                ( !mflag && (std::abs(s-b) >= (std::abs(c-d) * 0.5)) ) ||
                ( mflag && (std::abs(b-c) < TOL) ) ||
                ( !mflag && (std::abs(c-d) < TOL))  )
        {
            // bisection method
            s = (a+b)*0.5;

            mflag = true;
        }
        else
        {
            mflag = false;
        }

        fs = f(s);  // calculate fs
        d = c;      // first time d is being used (wasnt used on first iteration because mflag was set)
        c = b;      // set c equal to upper bound
        fc = fb;    // set f(c) = f(b)

        if ( fa * fs < 0)   // fa and fs have opposite signs
        {
            b = s;
            fb = fs;    // set f(b) = f(s)
        }
        else
        {
            a = s;
            fa = fs;    // set f(a) = f(s)
        }

        if (std::abs(fa) < std::abs(fb)) // if magnitude of fa is less than magnitude of fb
        {
            std::swap(a,b);     // swap a and b
            std::swap(fa,fb);   // make sure f(a) and f(b) are correct after swap
        }

    } // end for

    std::cout<< "The solution does not converge or iterations are not sufficient" << std::endl;

} // end brent_fun

Answer 1

Comments

OK I like the first comment.

/*******************************************************************************
*
* Grant Williams
*
* Version 1.0.0
* August 27, 2015
*
* Test Program for Brent's Method Function.
*
* Brent's method makes use of the bisection method, the secant method, and inverse quadratic interpolation in one algorithm.
*
* To Compile Please use icc -std=c++11 if using intel or g++ -std=c++11 if using GCC.
*
********************************************************************************/

You may want to add a copyright notice.
It is over my head the description but I assume somebody with some maths skills could understand it. So it provides some good information.

But there are a lot of other comments in this application are useless. Comments should be used judiciously within code. Bad comments are worse than no comments. You have to maintain comments in the same way you maintain the code and the compiler will not help you with the comments.

Your comments should not explain what the code is doing (the code should be easy to read and by self explanatory). You should explain WHY the code is doing something in a comment.

Example:

    fs = f(s);  // calculate fs
    d = c;      // first time d is being used (wasnt used on first iteration because mflag was set)
    c = b;      // set c equal to upper bound
    fc = fb;    // set f(c) = f(b)

These comments are all should be removed. The code should also be fixed so that it makes the comments redundant.

    fs = f(s);
    c  = upper_bound;
    fc = fb;

You'll notice that I even removed a line. That's because of after much clicking I found it was never used. What?!

    double a;               // lower bound
    double b;               // upper bound
    double TOL = 0.0001;    // tolerance
    double MAX_ITER = 1000; // maximum number of iterations

Much easier to write if you use names that define what they mean.

    double lowerBound;
    double upperBound;
    double tolerance     = 0.0001;
    double maxIterations = 1000;

Identifier Names

Using short identifier names is fine for maths. But when you write real applications that need to be maintained for decades you need people to understand your code. You can not rely on comments (because code changes over time and comments and code will drift apart). So you must use good names that explain their meaning for both variables and functions.

Single character names should be banned. A visual inspection of your code did not reveal any use of d. But I don't trust my eyes. So I had to search your code for occurrences using tools. There were so many false positives in the search because the single character d appears a lot in the middle of other words.

Pick names that are easy to see and spot and that provide meaning to the context that they are used in.

Reserved Words.

Certain styles of identifier are reserved for different Jobs. You should understand these reservations.

Identifiers that are all upper case are traditionally reserved for use by the pre-processor.

MAX_ITER

This will get clobbered by the pre-processors if you are unlucky. Pick names that use a combination of upper and lower case characters.

Note it is allso traditional to have type names (User Defined) start with an uppercase letter and all method/variable names start with a lower case letter.

Whether you prefer camel-casing your words or using _ between words is still a very divided subject. The _ is the minority group (but it is a big minority). I am in the group of camel-casing because using _ is error prone and hard to spot when you get it wrong.

Declare and initialize in one statement.

    std::clock_t start;
    start = std::clock();

    // It is best to initialize your variables on declaration.
    // It makes sure that value are always defined.

    std::clock_t start = std::clock();

Also don't declare variables until you need them. Its not much of a problem with POD types that don't run code. But objects that have a type and therefore a constructor will run the constructor at declaration.

Also for readability you want the variable declared just before use so that you can see they type of variable easily and not have to go searching to the top of the function just to find the type.

Prefer '\n' to std::endl

The std::endl places a '\n' on the stream then flushes the stream. Manually flushing the stream is usually a bad idea. The code will flush it for you at the best times. You should only manually flush the streams if there is something critical that needs to happen that the code does not understand directly (like the application could be crashing and will not flush the stream).

If you flush streams manually you are likely to find your IO very inefficient.

Declare loop variables so they are scoped correctly.

    for (function = 1; function <= 3; function++)

    // More usual to see this:
    for (int function = 0; function < 3; ++function)

1. Scope the loop variable to the for loop.
   This prevent you from leaking scope informaiton.
2. Most loops in C/C++ are < loops
   (Usually because we use zero indexed arrays.)
3. Prefer prefix increment
   For POD it does not matter. But for user defined types it does. If you change the type of the iterator later then you only need to change the type not change the type and then change all the postfix increments into prefix increments.

Your loop code seems to be an obtuse way of writing:

brents_fun(bradsFunction, -1.5, 0, tolerance, maxIteration);
brents_fun(alexsFunction, -10, 10, tolerance, maxIteration);
brents_fun(johnsFunction, -4,   3, tolerance, maxIteration);

Self documenting code

The body of this loop is way too long.

    for (unsigned int iter = 1; iter < MAX_ITER; ++iter)
    {

You should put the body in a function (with a descriptive name). Then your code becomes more self documenting.

Yep Crazy

       /*
            Crazy condition statement!:
            -------------------------------------------------------
            (condition 1) s is not between  (3a+b)/4  and b or
            (condition 2) (mflag is true and |s−b| ≥ |b−c|/2) or
            (condition 3) (mflag is false and |s−b| ≥ |c−d|/2) or
            (condition 4) (mflag is set and |b−c| < |TOL|) or
            (condition 5) (mflag is false and |c−d| < |TOL|)
        */

Also the comment does not help. The comment basically restates the code, which is absolutely no help as I can read the code and see what it says. Now when I read the comment and the code I have to verify that they match (which is a pain).

What would be better is to describe WHY there is a conditional tests here. What are you trying to achieve with the test?

In the code rather than have this obtuse expression.

        if (    ( (s < (3 * a + b) * 0.25) || (s > b) ) ||
                ( mflag && (std::abs(s-b) >= (std::abs(b-c) * 0.5)) ) ||
                ( !mflag && (std::abs(s-b) >= (std::abs(c-d) * 0.5)) ) ||
                ( mflag && (std::abs(b-c) < TOL) ) ||
                ( !mflag && (std::abs(c-d) < TOL))  )

Make the expression a function with a name that explains the test.

      if (falloutDegredationHasDroppedBelowThreshold(sign, falloutCoeffecient, dropLevel, plank)

Answer 2

A couple remarks on the benchmark itself:

• It is much too short. On my machine, it ends in approx 30 microseconds, with large variations between runs. You need to run your function more times to eliminate noise and to make sure you're well over whatever clock resolution your timer has.
• It is timing I/O, which you really don't want. I/O is slow and irregular. Unless that's specifically what you want to measure, make sure your tests do no I/O. e.g. I ran your code with 1 million iterations:
  • With the output in the tested function: ~2.4s
  • Without the output: ~0.6s

So increase the number of iterations, and remove all output from your function - it should have a way of returning what it found or signaling an error if it failed.

For the performance aspect, std::function has a drawback: there's pretty much no way for the compiler to do inlining, and the call can be more expensive than a plain function call (I'm not familiar with the details of std::function, but you'll find a lot of relevant info but searching for "std::function performance" or "overhead"). Whether that matters or not depends on the specific piece of code of course, but it does appear to matter here.

What you can try is using a template instead. The modification is very simple: change the signature from

void brents_fun(std::function<double (double)> f, ...)

to

template <typename Func>
void brents_fun(Func f, ...)

(And move the whole definition into a header while you're at it.)

On my machine, after having removed the output from the function and again for one million runs, this reduces the execution time from ~0.61s to ~0.48s - not huge, but not just line-noise either.

That's all I have on that front, hopefully someone else will chime in with other improvements.

---

Coding style looks good to me, applied consistently. Variable names are a bit cryptic but that's due to following the pseudo-code in the article. I wouldn't change them.

One thing I would change is the spacing in your function calls, e.g.

std::swap(a,b);
brents_fun(f,a,b,TOL,MAX_ITER);

to

std::swap(a, b);
brents_fun(f, a, b, TOL, MAX_ITER);

---

The "loop" in main is very very strange. If you need to do three different things one after the other, using a loop is a strange idea. Just call your three tests, the loop and if are just noise. Go back to a loop version if you put your test "parameters" in a vector.
You can use anonymous block if you want to isolate each test though.

{
  double a = -1.5;
  double b = 0;
  auto f = [](double x){ return (x+1) * (x+2) * (x+3); };
  brents_fun(f,a,b,TOL,MAX_ITER);
}
{
  double a = -10;
  double b = 10;
  auto f = [](double x){ return (x*x*x -4*x - 9); };
  brents_fun(f,a,b,TOL,MAX_ITER);
}

If you do keep a loop, declare the loop variable inside the for itself, and take the habit of starting at 0. (Array indexes are 0-based in C++ and quite a few other languages).

for (int test = 0; test < 3; ++test) { ... }

(Or use the for (element: container) form if you get the parameters in a container.)

Answer 3

    double s = 0;           // Our Root that will be returned

    for (unsigned int iter = 1; iter < MAX_ITER; ++iter)
    {
        // stop if converged on root or error is less than tolerance
        if (std::abs(b-a) < TOL)
        {
            std::cout << "After " << iter << " iterations the root is: " << s << std::endl;
            return;
        } // end if
        ...

looks like a bug to me. If lower_bound and upper_bound are within a tolerance of each other, the function will claim that 0 is root after 1 iteration.

• A crazy conditional statement could and should be simplified. Join both mflag branches, both !mflag branches and realize a common expression into a function:

  (mflag && somefunc(s-b, b-c, TOL) || !mflag && somefunc(s-b, c-d, tol)
  

  with

  static inline bool somefunc(double x, double y, double tol) {
      return std::abs(x) < std::abs(y)*0.5 || std::abs(y) < TOL;
  }
  

  Of course you need to come up with a descriptive name for somefunc.