Решение уравнения

Даны коэффициенты уравнения a₀*xⁿ + a₁*x^n-1 + ... + a_n = 0
Найти его корни

    Будем рассматривать уравнение как функцию. Исследуем её на монотонность: найдём производную и критические точки. Теперь просто найдём решение уравнения на отрезках между этими критическими точками за O(NlogL), так как на них функция монотонна.
    Производная будет на порядок меньше исходной функции количеством коэффициентов и степенью. Поэтому её нули находим рекурсивно. Как оптимизацию можно решение квадратных и линейных уравнений рассмотреть отдельно.
    Но есть один нюанс в этом алгоритме - это его точность. Экспериментальным путём было установлено, что при N = 40 алгоритм находит корни с точностью 2 знака после запятой. Для большей точности необходима поддержка длинных вещественных чисел.
    Асимптотика O((NlogL)^N-2).

Листинг C++

// решение уравнения
int solve_equation (list < double > & p, vector < double > & x)
{
    if (p.size () == 0)    // 0 = 0
        return - 1;
    if (eq (p.front (), 0))    // a = 0
    {
        p.pop_front ();
        return solve_equation (p, x);
    }
    if (p.size () == 1)
        return 0;
    if (eq (p.back (), 0))    // pn = 0
    {
        p.pop_back ();
        int flag = solve_equation (p, x);
        p.push_back (0.000);
        int i;
        forn (i, x.size ())
            if (eq (x[i], 0.000))
                return 1;
            else
                if (x[i] > 0)
                    break;
        x.insert (x.begin () + i, 0.000);
        return 1;
    }

    list < double > :: iterator it;
    if (p.size () == 2)    // ax + b = 0
    {
        it = p.begin ();
        double a = (* it);
        ++ it;
        double b = (* it);
        x.push_back (- b / a);
        return 1;
    }

    if (p.size () == 3)    // ax^2 + bx + c = 0
    {
        it = p.begin ();
        double a = (* it);
        ++ it;
        double b = (* it);
        ++ it;
        double c = (* it);
        double D = b * b - 4.000 * a * c;
        double a2 = 2.000 * a;
        if (eq (D, 0))
        {
            x.push_back (- b / a2);
            return 1;
        }
        if (D < 0)
            return 0;
        double sD = sqrt (D);
        x.push_back ((- b - sD) / a2);
        x.push_back ((- b + sD) / a2);
        if (x[0] > x[1])
            swap (x[0], x[1]);
        return 1;
    }

    int i, j;
    int n = p.size ();
    list < double > f; // f(x) = p'(x)
    it = p.begin ();
    forn (i, n - 1)
    {
        f.push_back ((double)(n - 1 - i) * (* it));
        ++ it;
    }
    vector < double > kx;
    int flag = solve_equation (f, kx);
    if (flag == - 1) return 0;

    kx.insert (kx.begin (), - 1e9);
    kx.insert (kx.end (), 1e9);
    forn (i, (int)kx.size () - 1)
    {
        double tx;
        if (bin_equation (p, kx[i], kx[i + 1], tx) == 1)
            x.push_back (tx);
    }
    return 1;
}

18:53
14.08.2009

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