# Finding root of strictly monotonic function by bisection

I need to find a root of strictly monotonic function f (only raising or only lowering). There is only a one root of a function. If given range from a to b doesn't contain root (fabs(f(a)*f(b) is not lower or equal than 0), we need to change range, since the assumption is that there is a one root, so we need to move on x-axis to left or right and check again, if f(a). If it is, then (f)*f(b) (e.g (-3)*2=-6, and that means there is a root between a and b).

I've made an solution and I would like to improve performance in terms of changing ranges (obviously adding +2 is not good for huge functions) vs halving ranges further, when finding f(c).

• What would be the best strategy?
• Also will there be any other problem to look on?
• What can I do to avoid user epsilon to be smaller than DBL_EPSILON?
• Any ideas to improve my code?
#include <iostream>
#include <cmath>

double f(double x)
{
//double y = std::pow(2, x) - 6;
double y = x - 10000;
return y;
}

bool isRaising(double a, double b)
{
if (f(a) > f(b)) return false;
return true;
}

bool areBothNumbersPositive(double a, double b)
{
if (f(a) > 0 && f(b) > 0) return true;
return false;
}

void moveLeft(double &a, double &b)
{
b = a;
a = a - 2;
}

void moveRight(double& a, double& b)
{
a = b;
b = b + 2;
}

void changeRanges(double &a, double &b)
{
if (!isRaising(a, b))
{
if (areBothNumbersPositive(a, b))
{
moveRight(a, b);
}
else
{
moveLeft(a, b);
}
}
else
{
if (areBothNumbersPositive(a, b))
{
moveLeft(a, b);
}
else
{
moveRight(a, b);
}
}
}

double findRoot(double a, double b, double epsilon)
{
if (epsilon < DBL_EPSILON) epsilon = DBL_EPSILON;
if (a == b) throw std::exception("Range must be higher than zero!");

double c;
if (std::fabs(f(a)) < epsilon) return a;
if (std::fabs(f(b)) < epsilon) return b;

while (f(a) * f(b) > 0)
{
changeRanges(a, b);
}
if (std::fabs(f(a)) < epsilon) return a;
if (std::fabs(f(b)) < epsilon) return b;
do
{
c = (a + b) / 2;
if (f(a) * f(c) > 0) a = c;
else b = c;
} while (!(std::fabs(f(c)) < epsilon));
return c;
}

int main()
{
//f(a) * f(b) > 0
std::cout << findRoot(-100,  -300,   0.00000000001) << std::endl;
std::cout << findRoot(100, 300, 0.00000000001) << std::endl;
std::cout << findRoot(-300, -100,  0.00000000001) << std::endl;
std::cout << findRoot(300, 100,  0.00000000001) << std::endl;

//f(a) * f(b) <= 0
std::cout << findRoot(-100, 300, 0.00000000001) << std::endl;
std::cout << findRoot(100, -300, 0.00000000001) << std::endl;
std::cout << findRoot(300, -100, 0.00000000001) << std::endl;
std::cout << findRoot(-300, 100, 0.00000000001) << std::endl;
}

• are you asking for f(x) = x - 10000 or any function double f(double) of which all you know is that it is monotonic? – slepic Nov 7 at 13:35
• Any function that is strictly monotonic (is always increasing or always decreasing), so it is guaranteed that it has any roots and only one root. – Jan Dycz Nov 7 at 23:03
• I know what it means for function to be monotonic. But in your code you Are only using the f function which Is defined as x - 10000. And you have no way to invoke your code over a differrent monotonic function other then changing the body of the f function. Thats why im asking if you only care about f the way it Is defined or you wanna be able to call it on any monotonic function taking double And returning double... – slepic Nov 8 at 6:31

### 1. Use established terms to make your code more readable to others

In this case, use "increasing" instead of "raising". In such a simple context it is quite obvious what you mean by "raising", but you would be surprised how much effort and how many misunderstandings you can save by using existing, well-known terms.

### 2. Give clear, unambiguous names to functions

When I see a function areBothNumbersPositive(double a, double b), I assume it returns a > 0 && b > 0. A more suitable name would be areBothFuncValuesPositive(double a, double b), even though it's a bit longer.

### 3. Don't return true or false based on the value of a boolean expression, return the expression directly

bool isIncreasing(double a, double b)
{
return f(a) < f(b);
}

bool areBothFuncValuesPositive(double a, double b)
{
return f(a) > 0 && f(b) > 0;
}


### 4. Duplicate code in chageRanges can be reduced

void changeRanges(double &a, double &b)
{
if (isIncreasing(a, b) == areBothFuncValuesPositive(a, b))
{
moveLeft(a, b);
}
else
{
moveRight(a, b);
}
}


Although I must admit, this seems a bit confusing. Maybe it would be a good idea to define a function which would figure out whether to go left or right, and then use this function in the if. This might be an overkill in this simple case, but is quite commonly used with more complex boolean expressions.

bool shouldMoveLeft(double a, double b)
{
return (isIncreasing(a, b) && areBothFuncValuesPositive(a, b))
||(isDecreasing(a, b) && areBothFuncValuesNegative(a, b));
}

void changeRanges(double &a, double &b)
{
if (shouldMoveLeft(a, b))
{
moveLeft(a, b);
}
else
{
moveRight(a, b);
}
}


### 5. Use functions where appropriate

Functions serve two main purposes: 1) code reuse; 2) conceptual abstraction. Even if you don't reuse them, they can make code more readable. For example, when you calculate f(a) * f(b) > 0, what you really want to know is whether the function has the same sign at both values, and neither f(a) nor f(b) are zero.

bool areBothFuncValuesNonzeroSameSign(double a, double b)
{
return f(a) * f(b) > 0;
}


An even better idea would be to introduce the concept of a sign, and then build upon it (I copied the sign function from this SO answer).

template <typename T> int sign(T val) {
return (T(0) < val) - (val < T(0));
}

{
return sign(x) == sign(y);
}

bool areBothNumbersNonzero(double x, double y)
{
return x != 0 && y != 0;
}

bool areBothFuncValuesNonzeroSameSign(double a, double b)
{
const double
fa = f(a),
fb = f(b);

return areBothNumbersNonzero(fa, fb) && areBothNumbersSameSign(fa, fb);
}


This is not as succinct and elegantly clever as your original code. But imagine it wasn't a single line, but an entire function / class / project written in such a clever way. Then you would likely have a bit of trouble understanding it when you come back to it after some time. It would be even more trouble for someone else who has never seen your code before and now has to understand it, because you don't work at the company anymore and they inherited your code, or because they are extending your project, using your library, etc. By keeping your code simple and "stupid", you make it more readable and accessible to more people, and therefore, more maintainable.

Another two important aspects of software, besides readability and maintainability, are composability and reusability. Like readability and maintainability, these two are related. The functions sign, areBothNumbersSameSign and areBothNumbersNonzero are quite generic, so they are easily combined together (composed). This improves their reusability. Make your functions as generic as possible, so they can be reused later in different contexts, and composed to form larger abstractions. This includes the use of generic types (templates), as in the sign function. Making the function closely tailored for one specific use case prevents it from being used in other cases. You can observe that the function areBothFuncValuesNonzeroSameSign can only be used in the context of the function f, since it depends on it. We can decompose its functionality into two logically distinct parts: 1) calculating the values f(a), f(b); and 2) checking if the computed values are both non-zero and have the same sign.

bool areBothNumbersNonzeroSameSign(double x, double y)
{
return areBothNumbersNonzero(x, y) && areBothNumbersSameSign(x, y);
}

bool areBothFuncValuesNonzeroSameSign(double a, double b)
{
return areBothNumbersNonzeroSameSign(f(a), f(b));
}


Again, this is just a simple example. But as the code grows bigger and the functionality gets more complex, such refactoring is necessary, although very often neglected.

### 6. Changing range in moveLeft and moveRight

As you pointed out, it is not ideal to use a constant (2) when adjusting the range. The first naive idea that would come to my mind is to simply shift the range, without changing its size.

void moveLeft(double &a, double &b)
{
const double diff = b - a;
b = a;
a = a - diff;
}

void moveRight(double& a, double& b)
{
const double diff = b - a;
a = b;
b = b + diff;
}


But, as you said, this is not very efficient. What you can do is increase the size of the range exponentially. I think (but I'm not sure), that this would converge more quickly on average.

void moveLeft(double &a, double &b)
{
const double diff = b - a;
b = a;
a = a - diff*2;
}

void moveRight(double& a, double& b)
{
const double diff = b - a;
a = b;
b = b + diff*2;
}


Note that I have not tested any of this, so I can't guarantee it works.

First thing first. Your assumption that a strictly monotonic function always has exactly one root is wrong. For it to be true, the function must not only be monotonic but also continuous on its entire domain (or in other words, it cannot have jump discontinuities, https://en.wikipedia.org/wiki/Classification_of_discontinuities#Jump_discontinuity).

EDIT: Actually it must also be defined on all real numbers, or in other words, its domain must be R.
EDIT2: Actually it can be discontinuous, but must not have jump discontinuities.

For example, consider f(x) => x + (x<0?-1:1). This function is strictly monotonic and is defined for any R, but it has a jump discontinuity and it does not have a root.

f(x) => (x < 0 ? x-1 : nan) is strictly monotonic on its entire domain, it is continuous on its entire domain, but its entire domain is not all R and it does not have a root.

f(x) => (x<0 ? x : 2*x) is strictly monotonic on its entire domain, it is not continuous, its entire domain is all R, but it does not have a jump discontinuity, and it has a root.

Now, let's assume that the function being continuous is part of your definition, but you just forgot to mention that.

Another problem I see is that function double findRoot(double,double,double) is accessing the f() function globally. And since f() is defined as (x) => x - 10000, the most efficient implementation of the findRoot() function is to simply return constexpr 10000.

If you want findRoot() to work over any double (*f)(double) you have to pass the function pointer to the function:

double findRoot(double (*f)(double), double, double, double);


I'm not going to say much more about your implementation, as @kyrill already said a lot. There is one thing I would object to in that answer though:

I.e. here:

bool shouldMoveLeft(double a, double b)
{
return (isIncreasing(a, b) && areBothFuncValuesPositive(a, b))
||(isDecreasing(a, b) && areBothFuncValuesNegative(a, b));
}


Although it's nice and readable, it's terrible in terms of performance. isIncreasing, areBothFuncValuesPositive, isDecreasing and areBothFuncValuesNegative - all these functions invoke f(a) and f(b) with the same a and b. This means that f(a) is invoked 4 times and returns the same result on each call. Same with f(b). Because f is an arbitrary function, it can be arbitrarily complex. And invoking it four times to get one result is four times less performant than it could be. And shouldMoveLeft is just one instance of this problem. Other functions have the same problem.

And one more point to make, although I haven't really investigated the consequences, but I just feel like there might be some. There can be functions which although strictly monotonic, change so slowly that ẟy is incredibly small even for a big ẟx. For example, (x) => x/lot where lot is a very big number. I leave that for you to find out what consequences it might have.

EDIT3: actually the opposite (very small ẟx yields large ẟy) could also cause some problems, e.g. (x) => lot * x.

• Yes, it would make sense to compute the values of f beforehand and define the function isIncreasing and others to work on their inputs directly instead of invoking f. Also, good point about f being used globally. Passing a function pointer is what occurred to me first, but I didn't bother to mess with the funny business of function pointers in C++. BTW which consequences in the case of y = (x) => x/lot do you have in mind? – kyrill Nov 8 at 23:17