一元三次方程的分治解法

/*  2011-09-21
    f(x)=ax^3+bx^2+cx+d 
    分析：设三个根为 x0，x1，x2 有 
            f(x)=a(x-x0)(x-x1)(x-x2)有(韦达定理)
            -b = a(x0+x1+x2)
             c = a(x0x1+x1x2+x2x3) 
            -d = ax0x1x2
            以 x0为已知量求解 x1，x2 
            x1 = 1/(2ax)(+sqrt((ax^2+bx)^2+4adx)-ax^2-bx) 
            x2 = 1/(2ax)(-sqrt((ax^2+bx)^2+4adx)-ax^2-bx) 
          
          接下来对 x在给定区间内进行二分即可 
         
         
    附韦达定理：
        ∑Xi=(-1)^1*A(n-1)/A(n) 　　
        ∑XiXj=(-1)^2*A(n-2)/A(n) 　　
            … 　　
        ΠXi=(-1)^n*A(0)/A(n)  
        {A(n)对应度数的系数}
*/ 

#include <iostream>
#include <cstdlib>
#include <cmath>

#define esp 1e-6

using namespace std;

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

double x1( double a, double b ,double c, double d, double x )
{
    return 1.0/(2.0*a*x)*(sqrt((a*x*x+b*x)*(a*x*x+b*x)+4.0*a*d*x)-a*x*x-b*x);
}

double x2( double a, double b ,double c, double d, double x )
{
    return 1.0/(2.0*a*x)*(-sqrt((a*x*x+b*x)*(a*x*x+b*x)+4.0*a*d*x)-a*x*x-b*x);
}

int main()
{
    double a,b,c,d;
    double A,B,M;
    while ( cin >> a >> b >> c >> d ) {
        cin >> A >> B;
        while ( B-A > esp ) {
            M = (A+B)/2.0;
            if ( f( a, b, c, d, M )*f( a, b, c ,d, A ) <= 0 )
                B = M;
            else A = M;
        }
        cout << M << " ";
        cout << x1( a, b, c, d, M ) << " ";
        cout << x2( a, b, c, d, M ) << " ";
        cout << endl;
    }
    return 0;
}