1
\$\begingroup\$

For those who aren't familiar with the bisection method for finding the root of a function (i.e. finding where \$f(x) = 0\$) the basic idea is:

  • Take a function \$f(x)\$ and an interval \$[a,b]\$
  • If \$[f(a)] \cdot [f(b)] < 0\$, there is a root in the interval, as one of those is a negative value
  • Repeatedly half the interval until you find an acceptable root within a maximum error of \$\varepsilon\$

I have the following code to emulate this process:

float bisect(float(*f)(float), float e, float a, float b){
    if (f(a)*f(b) > 0) { printf("Invalid interval - returning\n");  return 0; }
    float c;
    while (1){
        c = (a + b) / 2;
        if ( ((b-c) / c) <=e ) break;
        f(b)*f(c) <= 0 ? a = c : b = c;
    }
    return c;
}

Which I would pass a test function representing \$f(x) = x^{3}-x-2\$ with \$\varepsilon = 10^{-6}\$ and check for a root in the interval \$[a,b]\$.

typedef float (*func)(float);

float test(float x){
    return (x*x*x) - (x) - 2;
}

int main(){
    func f = &test;
    double start, end, ans;

    start = get_time();
    ans = bisect(f, .00001f, 1, 2);
    end = get_time();

    printf("The root is %f and took %e seconds", ans, end-start);
    return 0;
}

This code works, and returns 1.521385, which properly approximates the actual value of 1.5213928 (the actual error is 0.0000078).

I'm on Windows right now (using Visual Studio 2013), and I'm using this function for the benchmarking:

#include <windows.h>
double get_time()
{
    LARGE_INTEGER t, f;
    QueryPerformanceCounter(&t);
    QueryPerformanceFrequency(&f);
    return (double)t.QuadPart / (double)f.QuadPart;
}

I was interested to see if there were a way to optimize any of this to make it run as fast as possible.

\$\endgroup\$
  • \$\begingroup\$ So is it C or C++? \$\endgroup\$ – kraskevich Feb 21 '15 at 22:44
  • \$\begingroup\$ @kraskevich, Well it's C code, but I included the C++ tag because as far as I understand, the c++ compiler is being used (since I'm in Visual Studio) \$\endgroup\$ – galois Feb 21 '15 at 22:46
2
\$\begingroup\$
  • It is customary to use eps instead of e.

  • A ((b-c) / c) <=e condition doesn't look right (try \$f(x) = x\$ on the \$(-1, 1)\$ interval and observe a division by 0).

  • You seem to assume that b > a (otherwise the loop breaks immediately). Make sure to swap a and b if it is not so. Making sure that eps is positive also helps.

  • The code produces wrong results for a non-continuous \$f\$. I recommend to test that the function value at a suspected root is close enough to \$0\$.

  • Avoid break. while (b - a > eps) is much cleaner.

\$\endgroup\$
  • \$\begingroup\$ "A ((b-c) / c) <=e condition doesn't look right" - I think it would also be alright to do (b - c) < e, or f(b) < e ? \$\endgroup\$ – galois Feb 22 '15 at 3:37
  • \$\begingroup\$ These are two different stories (interval is narrow enough, and function value is close to 0 enough); both of them shall be taken into account. \$\endgroup\$ – vnp Feb 22 '15 at 5:44
0
\$\begingroup\$

Well you don't really need the multiplications to test if f(b) and f(c) have the same sign. And because you change only one endpoint of the interval on each iteration if you record the function values at a and b you can save one function call per iteration.

\$\endgroup\$

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.