www.pudn.com > Face3DModel.zip > cubicequation.cpp

//========================================================== 
//  Solve the equation a0 x^3 + a1 x^2 + a2 x + a3 = 0 
//      reference web page -  
//      http://mathworld.wolfram.com/CubicFormula.html 
//  History: 
//      Created on 01/24/05 by Xu Jiang 
 
#define PI 3.14159265358979 
 
#include  
 
double cubic_root(double x) 
{ 
    if (x >= 0) return pow(x, 1./3); 
    else return -pow(-x, 1./3); 
}    
 
int solve_equation_3(double &a0, double &a1, double &a2, double &a3, 
                     double &x1, double &x2, double &x3)  
{ 
    a1 /= a0; 
    a2 /= a0; 
    a3 /= a0; 
 
    double Q, R, D; 
    double a1_div_3 = a1 / 3; 
     
    Q = (3 * a2 - a1 * a1) / 9; 
    R = (9 * a1 * a2 - 27 * a3 - 2 * pow(a1, 3)) / 54; 
    D = pow(Q, 3) + pow(R, 2); 
 
    if (D <= 0)  
    { 
        double sqrt_minus_Q; 
        double theta; 
        if (Q != 0) 
        { 
            sqrt_minus_Q = sqrt(-Q); 
            theta = acos(-R / (Q * sqrt_minus_Q)); 
        } 
        else 
        { 
            sqrt_minus_Q = 0; 
            theta = 0; 
        } 
        x1 = -a1_div_3 + 2 * sqrt_minus_Q * cos(theta / 3); 
        x2 = -a1_div_3 + 2 * sqrt_minus_Q * cos((theta + 2 * PI) / 3); 
        x3 = -a1_div_3 + 2 * sqrt_minus_Q * cos((theta + 4 * PI) / 3); 
        return 3; 
    } 
    else  
    { 
        double sqrt_D = sqrt(D); 
        double S = cubic_root(R + sqrt_D); 
        double T = cubic_root(R - sqrt_D); 
        double S_plus_T = S + T; 
//        double t1 = -S_plus_T / 2 - a1_div_3; 
//        double t2 = (S - T) * i * sqrt(3) / 2; 
        x1 = S_plus_T - a1_div_3; 
//        x2 = t1 + t2; 
//        x3 = t1 - t2; 
        return 1; 
    } 
} 
 
double distance(int x1, int y1, int z1, int x2, int y2, int z2)  
{ 
    return((double)(x1 - x2) * (x1 - x2)  
        + (double)(y1 - y2) * (y1 - y2)  
        + (double)(z1 - z2) * (z1 - z2)); 
}