Решение кубического уравнения O(logN)

a * x³ + b * x² + c * x + d = 0

    Пусть f(x) = a * x³ + b * x² + c * x + d = 0,
тогда f'(x) = 3 * a * x² + 2 * b * x + c, найдём критические точки.
    Функция в области между критическими точками монотонно возрастает или убывает. Поэтому уравнение на этих отрезках можно решить с помощью бинарного поиска.

Листинг C++

int sign (double x)
{
    return eq (x, 0) ? 0 : x < 0 ? - 1 : 1;
}

bool cubic_binfind (double l, double r, double a, double b, double c, double d, double & x)
{
    double fl = ((a * l + b) * l + c) * l + d;
    double fr = ((a * r + b) * r + c) * r + d;
    if (sign (fl) == 0)
    {
        x = l;
        return 1;
    }
    if (sign (fr) == 0)
    {
        x = r;
        return 1;
    }
    if (sign (fl) == sign (fr))
        return 0;


    while (r - l > eps / 1000.000)
    {
        double mid = 0.500 * (l + r);

        double fmid = ((a * mid + b) * mid + c) * mid + d;

        if (sign (fmid) == 0)
        {
            x = mid;
            return 1;
        }

        if (sign (fmid) == sign (fr))
            r = mid;
        else
            l = mid;
    }

    x = 0.5 * (r + l);

    return 1;
}

int cubic_equation (double a, double b, double c, double d, double & x1, double & x2, double & x3)
{
    if (eq (a, 0))
        return sqr_equation (b, c, d, x1, x2);

    double xk1, xk2;    // критические точки функции f(x) = ax^3 + bx^2 + cx + d
    int cnt = sqr_equation (3.000 * a, 2.000 * b, c, xk1, xk2);
    if (cnt == 1)
        xk2 = xk1;
    if (cnt == 0)
        xk2 = xk1 = 0;

    double left = - inf;
    double right = inf;

    double xp1, xp2, xp3;
    int c1 = cubic_binfind (left, xk1, a, b, c, d, xp1);
    int c2 = cubic_binfind (xk1, xk2, a, b, c, d, xp2);
    int c3 = cubic_binfind (xk2, right, a, b, c, d, xp3);

    if (c1)
        x1 = xp1;
    if (c2)
        if (c1)
            x2 = xp2;
        else
            x1 = xp2;

    if (c3)
        if (c2)
            if (c1)
                x3 = xp3;
            else
                x2 = xp3;
        else
            if (c1)
                x2 = xp3;
            else
                x1 = xp3;

    return c1 + c2 + c3;
}

14:57
06.08.2009

По всем вопросам обращаться: rumterg@gmail.com