Вначале немного теории. Как обычно, будем рассматривать уравнение в котором коэффициент при x³ равен нулю:

x4 + px2 + qx + r = 0    (1)

Представим это уравнение в виде произведения двух квадратных:

x4 + px2 + qx + r = ( x2 + ax + b ) ( x2 - ax + c )    (2)

Тогда для коэффициентов a, b, c получим следующую систему уравнений:

b + c - a2 = p
a ( c - b ) = q    (3)
bc = r

Из первых двух уравнений системы (3) получаем:

2c = p + a2 + q/a    (4)
2b = p + a2 - q/a

Затем подставляем выражения для b и c в третье уравнение системы (3) и получаем уравнение шестой степени относительно a:

a6 + 2pa4 + (p2 - 4r)a2 - q2 = 0    (5)

Сделаем подстановку: y = a².

y3 + 2py2 + (p2 - 4r)y - q2 = 0    (6)

Таким образом получается следующий метод решения уравнения четвёртой степени (1). Вначале решаем кубическое уравнение (6). Берём произвольный корень и получаем a. Затем подставляем a в (4) и находим b и c. Остаётся только решить два квадратных уравнения из (2) и получить четыре корня.

Для нахождения действительных корней уравнения четвёртой степени предлагаю следующие две функции:

typedef unsigned int nat;

// x^4 + p * x^2 + q * x + r = 0

nat root4s ( double p, double q, double r, double * x )
{
    if ( r == 0 ) // это условие добавлено 21.10.2008
    {
        *x = 0;
        return 1 + root3s ( p, q, x+1 );
    }
    if ( q == 0 )
    {
m1:     double t[2];
        nat n = root2 ( p, r, t );
        if ( n == 0 ) return 0;
        if ( t[0] >= 0 )
        {
            x[0] = sqrt ( t[0] );
            x[1] = - x[0];
            n = 2;
        }
        else
        {
            n = 0;
        }
        if ( t[1] >= 0 )
        {
            x[n] = sqrt ( t[1] );
            x[n+1] = - x[n];
            n += 2;
        }
        return n;
    }
    nat n = root3 ( p + p, p * p - 4. * r, -q * q, x );
    double a = *x;
    if ( n == 3 )
    {
        if ( a < x[1] ) a = x[1];
        if ( a < x[2] ) a = x[2];
    }
    if ( a <= 0 ) goto m1; // этот переход сделан 09.12.2011
    p += a;
    a = sqrt ( a );
    const double b = q / a;
    n = root2 (  a, 0.5 * ( p - b ), x   );
    n+= root2 ( -a, 0.5 * ( p + b ), x+n );
    return n;
}

// x^4 + a * x^3 + b * x^2 + c * x + d = 0

nat root4 ( double a, double b, double c, double d, double * x )
{
    if ( a == 0 )
    {
        return root4s ( b, c, d, x );
    }
    if ( d == 0 )
    {
        *x = 0;
        return 1 + root3 ( a, b, c, x+1 );
    }
    const double e = a / 4;
    const double h = e * e;
    const double p = -6 * h + b;
    const double q =  8 * h * e - 2 * b * e + c;
    const double r = -3 * h * h + b * h - c * e + d;
    const nat n = root4s ( p, q, r, x );
    for ( nat i = 0; i < n; ++i ) x[i] -= e;
    return n;
}

Первая из них решает "неполное" уравнение, т.е. у которого коэффициент при x³ равен нулю, а вторая сводит общий случай к предыдущему.
Исходники этих функций находятся в файле mathem.cpp.

Следующие две функции, кроме вещественных корней, вычисляют также комплексные корни:

void root4s ( double p, double q, double r, 
              Complex & x1, Complex & x2, Complex & x3, Complex & x4 );

void root4 ( double a, double b, double c, double d, 
             Complex & x1, Complex & x2, Complex & x3, Complex & x4 );

void root4s ( const Complex & p, const Complex & q, const Complex & r, 
              Complex & x1, Complex & x2, Complex & x3, Complex & x4 );

void root4 ( const Complex & a, const Complex & b, const Complex & c, const Complex & d, 
             Complex & x1, Complex & x2, Complex & x3, Complex & x4 )