Чаще всего кубические уравнения решают "тригонометрическим" методом, но там надо вычислять специальные функции, которые по сути являются подпрограммами, хотя и весьма оптимизированными. Можно поступить по другому - вначале найти один корень методом Ньютона, а затем, используя полученный корень, свести уравнение к квадратному для поиска остальных корней.

// x^3 + a * x + b = 0

static void cubic ( double a, double b, double & x )
{
    double s = 1.;
    while ( b + a > -1. )
    {
        a *= 4.;
        b *= 8.;
        s *= 0.5;
    }
    while ( b + a + a < -8. )
    {
        a *= 0.25;
        b *= 0.125;
        s *= 2.;
    }
    x = 1.5;
    for ( int i = 9; --i >= 0; )
    {
        x -= ( x * ( x * x + a ) + b ) / ( 3. * x * x + a );
    }
    x *= s;
}

// x^3 + p * x + q = 0

int root3s ( double p, double q, double * x )
{
    if ( q < 0 )
    {
        cubic ( p, q, *x );
    }
    else
    if ( q > 0 )
    {
        cubic ( p, -q, *x );
        *x = - *x;
    }
    else
    {
        *x = 0;
    }
    return 1 + root2 ( *x, p + (*x)*(*x), x+1 );
}

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

int root3 ( double a, double b, double c, double * x )
{
    if ( c == 0 )
    {
        *x = 0;
    }
    else
    {
        const double a3 = a / 3.;
        double p = b - 3.* a3 * a3;
        double q = c - ( a3 * a3 + p ) * a3;
        if ( q < 0 )
        {
            cubic ( p, q, *x );
        }
        else
        if ( q > 0 )
        {
            cubic ( p, -q, *x );
            *x = - *x;
        }
        else
        {
            *x = 0;
        }
        *x -= a3;
        const double t = *x * ( *x * 3. + a * 2. ) + b;
        if ( fabs ( t ) > 1e-3 )
        {
            *x -= ( *x * ( *x * ( *x + a ) + b ) + c ) / t;
        }
        a += *x;
        b += *x * a;
    }
    return 1 + root2 ( a, b, x+1 );
}

Здесь рассматриваются два типа уравнений: общее (root3) и "неполное" (root3s). Обе функции возвращают количество найденных вещественных корней, а сами корни помещаются в массив по указателю x. Вначале отыскивается один корень ( а он всегда существует) методом Ньютона. Перед этим параметры уравнения путём масштабирования меняются таким образом, чтобы корень лежал между 1 и 2. В этом случае экспериментально (несколько миллионов испытаний) я установил, что достаточно 9 итераций. На самом деле нужно даже меньше, но на всякий случай пусть будет 9. После того как найден один корень кубическое уравнение сводится к квадратному.

В первоначальном варианте функции root3 уточнения корня методом Ньютона не было. Но затем выяснилось, что в некоторых случаях ( например, когда а3 и х после cubic очень близки ) уточнение необходимо. Константа 1е-3 выбрана произвольно.

Исходники находятся в mathem.zip.