cubic function root bisection search time limit exceeded

Question

I have implemented this solution for finding a root of a cubic function

f(x) = ax³ + bx² + cx + d

given a, b, c, and d, ensuring it's being monotonic.

After submitting the solution to an online judge without being shown the test cases, I am being faced by a time limit error. a, b, c, and d guarantee that the function is monotonic and we know it is being continuous. The code first finds the interval [A, B] such that f(A) * f(B) < 0; then the code moves to implement the bisection search.

What I want to know is if there is some possibility to minimize the time complexity of my code so it passes the online judge. The input is a, b, c, d, and the output should be the root with an error 0.000001.

Code:

#include <iostream>
#include <algorithm>
//#include <cmath>
//#include <string>

using namespace std;

int f(double a, double b, double c, double d, double x) {
    return x*(x*(a*x + b) + c) + d;
}

int main() {

    freopen("input.txt", "r", stdin);
    freopen("output.txt", "w", stdout);

    double a, b, c, d, A, B, x = 1, res;
    cin >> a >> b >> c >> d; 

    //determinning the interval
    double f_x = f(a, b, c, d, x);
    if (a > 0) { // strictly increasing
        if (f_x > 0) { B = 0;
            while (f(a, b, c, d, x) >= 0) { x -= x; }
            A = x; }
        else { A = 0;
            while (f(a, b, c, d, x) <= 0) { x += x; }
            B = x; }
    }

    else { //strictly decreasing
        if (f_x > 0) { A = 0;
            while (f(a, b, c, d, x) >= 0) { x += x; }
            B = x; }
        else { B = 0;
            while (f(a, b, c, d, x) <= 0) { x -= x; }
            A = x; }    
    }
    // Bisection Search
    double l = A;
    while ((B - A) >= 0.000001)
    {
        // Find middle point 
        l = (A + B) / 2;

        // Check if middle point is root 
        if (f(a, b, c, d, l) == 0.0)
            break;

        // Decide the side to repeat the steps 
        else if (f(a, b, c, d, l)*f(a, b, c, d, A) < 0)
            B = l;
        else
            A = l;
    }
    res = l;
    cout.precision(6);
    cout << fixed << " " << res;

    return 0;
}

`if (f(a, b, c, d, l) == 0.0)` -- Checking for exactly 0.0 is not a good idea. [See this](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) — PaulMcKenzie, Nov 15 '19 at 19:10
i saw the code already earlier today. `x -= x;` is still wrong ;) — 463035818_is_not_an_ai, Nov 15 '19 at 19:11
@formerlyknownas_463035818 It has been corrected, please tell me where, it prints correct results now, you can check if you have a b c d to make the function monotonic. — v_head, Nov 15 '19 at 19:16
`x -= x;` is the same as `x = 0;` and if that wont make you leave the while loop you never will leave it — 463035818_is_not_an_ai, Nov 15 '19 at 19:17
`x -= x;` is exactly the same as `x = x - x;` which is... `x = 0;` — 463035818_is_not_an_ai, Nov 15 '19 at 19:19
formerlyknownas_463035818 is telling you. Change to `x = -x`. — Ripi2, Nov 15 '19 at 19:20
when writing an iteration with maths involved I would always implement a maximum loop count. Any mistake you make can make the iteration not converge. This can be caught easily if you limit the number of iterations — 463035818_is_not_an_ai, Nov 15 '19 at 19:21
@Ripi2 no I am not! I am just saying that `x = 0;` cannot be correct — 463035818_is_not_an_ai, Nov 15 '19 at 19:21
@formerlyknownas_463035818 Strictly speaking, yes, you told what you told. But you also pointed to the bug. — Ripi2, Nov 15 '19 at 19:24
@Ripi2 `x = -x;` would iterate forever as well. I really dont know what OP wanted to actually choose as next step — 463035818_is_not_an_ai, Nov 15 '19 at 19:25
@formerlyknownas_463035818 The problem with these online judges is that we don't know the input. For test, I suggest to take a never null 2nd order polynomial and to integrate it to get a monotonic one — Damien, Nov 15 '19 at 20:02
@Damien the problem is not really the online judges tests, but the absense of any other known reproducible tests, just saying ;) — 463035818_is_not_an_ai, Nov 15 '19 at 20:05

Evg · Answer 1 · 2019-11-15T20:18:38.173

There is no need to determine the initial interval, just take [-DBL_MAX, +DBL_MAX]. The tolerance can be chosen to be 1 ULP.

The following code implements these ideas:

// This function will be available in C++20 as std::midpoint
double midpoint(double x, double y) {
    if (std::isnormal(x) && std::isnormal(y))
        return x / 2 + y / 2;
    else
        return (x + y) / 2;
}

int main() {
    ...
    const auto fn = [=](double x) { return x * (x * (x * a + b) + c) + d; };

    auto left  = -std::numeric_limits<double>::max();
    auto right =  std::numeric_limits<double>::max();

    while (true) {
        const auto mid = midpoint(left, right);
        if (mid <= left || mid >= right)
            break;

        if (std::signbit(fn(left)) == std::signbit(fn(mid)))
            left = mid;
        else
            right = mid;
    }

    const double answer = left;
    ...
}

Initially, fn(x) can overflow and return inf. No special handling of this case is needed.

hard to beat, any better start interval will just save only a couple of iterations — 463035818_is_not_an_ai, Nov 15 '19 at 20:09
Very impressive. A little room for speed improvment: here function is called twice per iteration. Could be handled easily if needed — Damien, Nov 15 '19 at 20:51

cubic function root bisection search time limit exceeded

1 Answers1