bisection root finding segmentation fault

Question

I use following code to find roots of functions using a simple bisection algorithm

#include <stdio.h>
#include <stdlib.h>

typedef float (*continous_function)(float);

static const float epsilon = 0.0000001;

#if __STDC__
static __inline float fabsf(float x)
{
        return x >=0 ? x : -1*x;
}
#else
#include <math.h>
#endif

float poly_1(float x)
{
        return x*x*x+2*x*x+3;
}

float poly_2(float x)
{
        return 2*x + 3;
}

float poly_3(float x)
{
        return 5*x*x*x*x*x+3*x*x*x+2*x+7;
}

static __inline void swap_in_place(float *a, float *b)
{
        float tmp = *a;
        *a = *b;
        *b = tmp;
}

float bisection_root(continous_function f, float a, float b)
{

        float neg = f(a);
        float pos = f(b);
        float c = 0.5 * (a + b);
        float mid;

        if (neg * pos > 0) {
                /* neg and pos should have different sizes */
                abort();
        }

        if (neg > 0) {
                /* Ensure f(a)=neg is negative */
                swap_in_place(&neg, &pos);
                swap_in_place(&a, &b);
        }

        mid = f(c);

        if (fabsf(mid) < epsilon) {
                return c;
        } else {
                if (mid > 0)
                        return bisection_root(f, a, c);
                else
                        return bisection_root(f, c, b);
        }

}

int main()
{
        {
                float a = -5;
                float b = 2;
                printf("Root of x^3+2*x^2+3 is %f\n", bisection_root(poly_1,a,b));
        }

        {
                float a = -2;
                float b = 0;
                printf("Root of 2*x+3 is %f\n", bisection_root(poly_2,a,b));
        }

        {
                float a = -1;
                float b = 0;
                printf("Root of 5*x^5+3*x^3+2*x+7 is %f\n", bisection_root(poly_3,a,b));
        }

        return 0;
}

This program cause segmentation fault when compiled on windows xp (32-bit) using mingw gcc.

When the number of decimal digits of epsilon decreased, segmentation fault can be avoided. Therefore, I conclude that segmentation fault have something to do with an overflow.

I would like to know why and how exactly this error occurs so that I can find a reliable way to set epsilon or fix other errors that might cause the problem.

Answer 1

static const float epsilon = 0.0000001;

There does not need to be a float for which the value of your polynomial is at distance epsilon from 0. In particular if the root is for largish values of x, the polynomial may jump from being less than -epsilon to being more than +epsilon between one float and its successor.

In this case, and in many others, your algorithm will loop, Since you implemented it with recursion and C compilers usually do not guarantee tail-call optimization, this infinite loop can result in a segmentation fault (when the stack becomes full).

The solution, applicable whether you use recursion or a simple while loop, is to limit the number of iterations.

Note: you should read about Horner's scheme for computing polynomials.