Why segmentation fault in recursion to find roots of a function?

Question

Beginner in C.

I working on exercise to find roots of a function f(x) with newton method.

This is my code:

float bisection_root(float a, float b, float delta) {

    // TODO > pass functions ?
    // TODO check f(a) * f(b) < 0 o non finisce mai...
    // 

    float func(float x);

    printf("test %f \n",  func(a) * func(b) );


    // mean:
    float mean = (a + b) / 2.0;

    if ( fabs(func(mean)) < delta) {
        printf("Called function is: %s. Found delta. Result: %f\n",__func__, mean);
        return mean;
    } 
    else 
    // check if a, mean have same signs:
    if ((func(a) * func(mean)) < 0) {
        printf("Called function is: %s. Taking first half. Result: %f %f\n",__func__, delta, func(a) * func(mean));
        return bisection_root(a, mean, delta);
    }

    else {
        printf("Called function is: %s. Taking second half. Result: %f %f\n",__func__, delta, func(mean) * func(b));
        return bisection_root(mean, b, delta);
    }
    

}

float func(float x) {
    return 2*pow(x,3) - 4*x + 1;
}

and called in main like:

int main(int argc, char *argv[]) {
    
    bisection_root(0, 1, 0.00000000002);
    return 0;
}

The result is that recursion does not stop, and I have a segmentation fault: 11.

Could you help in understanding why recursion fail and where is that "program writes out of memory" ?

What does negative zero means in C ?

This is my output:

Called function is: bisection_root. Taking first half. Result: 0.000000 -0.750000
test -0.750000 
Called function is: bisection_root. Taking second half. Result: 0.000000 -0.023438
test -0.023438 
Called function is: bisection_root. Taking first half. Result: 0.000000 -0.012329
...

Called function is: bisection_root. Taking first half. Result: 0.000000 -0.000000 0.258652
test -0.000000 
Called function is: bisection_root. Taking first half. Result: 0.000000 -0.000000 0.258652
test -0.000000 
Called function is: bisection_root. Taking first half. Result: 0.000000 -0.000000 0.258652
Segmentation fault: 11

You are going in a stack overflow because your call recursively to many time you function. — OznOg, Sep 02 '21 at 14:54
Using text rather than picures of text makes for a better post. — chux - Reinstate Monica, Sep 02 '21 at 15:00
Note that your code checks if `func(mean)` is close enough to zero, but it doesn't check if `mean` is numerically different from either `a` or `b`. When implementing a bisection method it's ususally a good idea to stop when the interval becomes smaller than a target value. — Bob__, Sep 02 '21 at 15:20
Also, it would be better to put the declaration of `func` *outside* the function `bisection_root`. — Bob__, Sep 02 '21 at 15:29
Change your `%f` conversions to `%g`. This will adapt to print small numbers with exponential notation, so you will see better representations of their values than “0.000000”. — Eric Postpischil, Sep 02 '21 at 15:36
Instrument the `bisection_root` to print the ongoing values of `a` and `b` and the values of `mean` and `func(mean)` as it calculates them. This would make it easier for you to follow what the function is doing than printing the more complicated value `func(mean) * func(b)`. You would see the function is getting stuck because `func(mean)` cannot go below a certain threshold due to arithmetic rounding errors. — Eric Postpischil, Sep 02 '21 at 15:38
*"I['m] working on exercise to find roots of a function f(x) with newton method."* The posted code tries to implement the bisection method, is it a fallback in case the Newton method (which is quite different) doesn't converge? — Bob__, Sep 02 '21 at 16:09
@Bob__ why would it be better? Coming from python, i thought - if I declare within the function, I know it 's only used there, kinda belongs to it. Advice of best practice are welcome — user305883, Sep 02 '21 at 21:19
I imagined that, but it doesn't really work the same way in C (see e.g. https://stackoverflow.com/questions/2608158/nested-function-in-c ). You may consider a different design pattern, decoupling the algorithm from a specific target function: https://godbolt.org/z/jjW486qhW — Bob__, Sep 02 '21 at 22:07
@user305883 `float func(float x)` is not defined as `static`. Its name is in the global name space. Declaring it within `bisection_root()` does not make it more local. — chux - Reinstate Monica, Sep 02 '21 at 22:08
BTW, about the *"What does negative zero means in C?"* question, that's an artifact of the internal representation of floating-point values (see e.g. https://stackoverflow.com/questions/5095968/does-float-have-a-negative-zero-0f or even https://en.wikipedia.org/wiki/IEEE_754). — Bob__, Sep 02 '21 at 22:16

chux - Reinstate Monica · Answer 1 · 2021-09-02T22:45:07.000

2

[Update]

why recursion fail

delta as 0.00000000002 is too small an error limit for this task.

In this case, float mean = (a + b) / 2.0; keeps generating a mean were fabs(func(mean)) < delta is never true.

Consider when a=0.258652002f and b=0.258652031f, b is the next possible float after a. a, b straddles the mathematical root of 0.25865202250415276.... Their func() values are +7.5155981e-08 and -3.20905009e-08, both magnitudes exceeding 2e-09. There is no solution within that error bound.

The change from pow() to powf() (in below earlier answer) simply changed func() a little, enough to sometimes return a small enough result.

With the bisection method, recursion could simply stop when the mean is equal to a or b. delta not really needed.

//if (fabs(func(mean)) < delta) {
if (mean == a || mean == b) {

What does negative zero means in C ?

OP had output like "-0.000000". This occurs when the value was between -0.0000005 and zero and printed with "%f". (C supports 2 zeros, +0.0 and -0.0, so "-0.000000" also comes up when the value is -0.0. This is seen when a calculation rounds to 0.0 from the negative side.)

For debug, better to use "%g" than "%f" to see more information when the value is small.

Other thoughts

If sticking to float math, float mean = (a+b) / 2.0f; better as float mean = a/2 + b/2 to avoid overflow.

Rather than 2*pow(x,3) - 4*x + 1, consider (2*x*x - 4)*x + 1. Mathematically more stable.

See @Bob__ for 2 good links concerning these ideas.

[Earlier answer]

It seemed more sensible to use float math rather than double math for a float problem.

pow() --> powf()
2.0 --> 2.0f
fabs() --> fabsf()
0.00000000002 --> 0.00000000002f

Curiously, using

pow() --> powf()

lead to completion.

Called function is: bisection_root. Found delta. Result: 0.258652`

edited Sep 02 '21 at 22:45

answered Sep 02 '21 at 15:09

chux - Reinstate Monica

143,097
13
135
256

1

yes, at some point f(a) *f(b) underflows floats – OznOg Sep 02 '21 at 15:13
That changing to `float` arithmetic in `func` (the rest is irrelevant; e.g., `fabs` and `fabsf` return the same value for `float` arguments) causes termination in this case is largely coincidence. 0.00000000002 is a small value to expect to achieve. We should do some analysis to show how small we might expect `func` to get near roots that are not representable. Since 2x^3−4x+1 is zero, 2x^3-4x is near −1 when we are near a root, so we are likely to have errors around the ULP of 1. – Eric Postpischil Sep 02 '21 at 15:18
1

@OznOg: There is no underflow here; the numbers involved are far from where `float` underflows, below 2^−149. `a`, `b`, and `mean` go to about .258652, and the function value goes to about −3.20905•10^−8. – Eric Postpischil Sep 02 '21 at 15:24
@EricPostpischil Agree'd. Answer up dated. As I see it, rather than using an error bound on the `|func(root)|`, easier and better to bisect to the ULP of `root`. – chux - Reinstate Monica Sep 02 '21 at 16:14
Could you please explain what "more stable" means, about 2*x*x vs 2*pow(x,3) ? Why float mean = (a+b) / 2.0f; is better than float mean = (a+b) / 2.0; ? Could you elaborate why in my solution the problem goes to stack overflow? IS it because the resolution of the error is to big or to small ? – user305883 Sep 02 '21 at 21:27
1

@user305883 "to avoid overflow" measn "to avoid numeric overflow". If `a + b` exceeds `FLT_MAX` or `-FLT_MAX`, then `(a+b) / 2.0f` is +/- infinfity. – chux - Reinstate Monica Sep 02 '21 at 21:40
1

@user305883 "(a+b) / 2.0f; is better than float mean = (a+b) / 2.0;" --> `(a+b)` is a `float` addition. `/ 2.0` the takes that sum, converts to `double` and then does a `double` division. The quotient is then converted back to `float` to assign to `mean`. `(a+b) / 2.0f` avoid all the `double` to/from `float` conversison. – chux - Reinstate Monica Sep 02 '21 at 21:42
1

@user305883 `pow(x,3)` vs `powf(x,3)` --> Has similar `float/double` issues. `pow()` is a `double` function, when `powf()` is likely sufficent here. – chux - Reinstate Monica Sep 02 '21 at 21:51
1

@user305883 Consider `a*a – b*b` done that way or `(a-b)*(a+b)`. Mathematically equivalent. When `a` is near `b`, `a*a – b*b` suffer precision loss forming each product and the subtraction amplifies that loss. `(a-b)*(a+b)` incurs little precision loss. `pow(x,3)` is not as tightly defined as `x*x*x`. These 2 factors together make `(2*x*x - 4)*x + 1` superior to `2*pow(x,3) - 4*x + 1`. It is not easy to pack courses in [numerical analysis](https://en.wikipedia.org/wiki/Numerical_analysis) in a short paragraph, but hope this helps. – chux - Reinstate Monica Sep 02 '21 at 22:02
2

@user305883 If I may, that's the [Horner's method](https://en.wikipedia.org/wiki/Horner%27s_method) for evaluating polynomials. About the evaluation of `mean`, you may find this talk (even if it's about a C++ feature) ["std::midpoint? How Hard Could it Be?"](https://www.youtube.com/watch?v=sBtAGxBh-XI) enlighting. – Bob__ Sep 02 '21 at 22:25
@chux-ReinstateMonica wow thank you so much Monica for your informative comments and tips. I'm now able to complete: the thing I changed was the condition `(mean == a || mean == b)` to exit recursion. But still not clear to me these things: 1. exercise shows exit condition on `f(y) < d`. If I used a `d = 0.02`, still I would have got the same error `Segmentation fault: 11`. So why overflow, if `d` was great enough this time? 2. Double is "greater" than float. `Pow()` return a double, which is a number with a slot greater than number returned by `Powf()`. So why `pow()` would not work ? – user305883 Sep 03 '21 at 12:56
1

@user305883 "exercise shows exit condition on f(y) < d. If I used a d = 0.02, still I would have got the same error Segmentation fault: 11" is unclear. `bisection_root(0, 1, 0.02);` leads to `Called function is: bisection_root. Found delta. Result: 0.257812`. – chux - Reinstate Monica Sep 07 '21 at 01:50

score -1 · Answer 2 · edited Sep 02 '21 at 15:08

A recursive function works on a bounded array.

First thing to look at is "the basis step correct and is the recursion stopping when it's supposed to" or, to set the base case correctly.
Second most common error while using recursion is, to set the base case for one variable and traverse the array using another reference variable. Thus making the base case unreachable.

Check for these 2 issues in your code.

Why segmentation fault in recursion to find roots of a function?

2 Answers2