SP_Maths Namespace Reference

Functions
void	CrossProduct (const double A[3], const double B[3], double C[3])
double	DotProduct (const double A[3], const double B[3])
void	eigenValsAndVecsOf3by3SymMatrix (double m11, double m12, double m13, double m22, double m23, double m33, double &l1, double &l2, double &l3, double v1[3], double v2[3], double v3[3])
void	eigenValsAndVecsOf3by3SymMatrixUsingQRFactor (double m11, double m12, double m13, double m22, double m23, double m33, double &l1, double &l2, double &l3, double v1[3], double v2[3], double v3[3])
double	eigenValsOf3by3SymMatrix (double m11, double m12, double m13, double m22, double m23, double m33, double &lamda1, double &lamda2, double &lamda3)
double	eigenValsOf3by3SymMatrixSmith1961 (double m11, double m12, double m13, double m22, double m23, double m33, double &lamda1, double &lamda2, double &lamda3)
void	FindAnEigenVector (double &m11, double &m12, double &m13, double &m22, double &m23, double &m33, double l1, double v1[3])
void	FindOrthogonalVectors (const double iA[3], double oB[3], double oC[3])
double	RootFindingByNewtonRaphsonMethod (void(*funcd)(double x, double *fx, double *dfx), double x1, double x2, double xacc)
double	rootOfCubicPolynomial (double a, double b, double c, double d, int rootNo)
double	rootsOfCubicPolynomial (double a, double b, double c, double d, double &r1, double &r2, double &r3)
double	rootsOfCubicPolynomial (double a, double b, double c, double d, std::complex< double > &r1, std::complex< double > &r2, std::complex< double > &r3)
void	rootsOfCubicPolynomialCodeGeneratedByMaple (double a, double b, double c, double d, std::complex< double > &r1, std::complex< double > &r2, std::complex< double > &r3)
Variables
const double	PI = 3.141592654

---

Function Documentation

void SP_Maths::CrossProduct	(	const double	A[3],
		const double	B[3],
		double	C[3]
	)

Definition at line 21 of file sp_mathematics.cpp.

    {
        C[0] = A[1]*B[2] - A[2]*B[1];
        C[1] = A[2]*B[0] - A[0]*B[2];
        C[2] = A[0]*B[1] - A[1]*B[0];
    } // END FN:    CrossProduct()

double SP_Maths::DotProduct	(	const double	A[3],
		const double	B[3]
	)

Definition at line 31 of file sp_mathematics.cpp.

    {
        return  A[0]*B[0] + A[1]*B[1] + A[2]*B[2];
    } // END FN:    DotProduct()

void SP_Maths::eigenValsAndVecsOf3by3SymMatrix	(	double	m11,
		double	m12,
		double	m13,
		double	m22,
		double	m23,
		double	m33,
		double &	l1,
		double &	l2,
		double &	l3,
		double	v1[3],
		double	v2[3],
		double	v3[3]
	)

Definition at line 526 of file sp_mathematics.cpp.

    { 
        double errThreshold = 1E-05;    // error tolerant value

        double discriminant = eigenValsOf3by3SymMatrix( m11, m12, m13, m22, m23, m33, l1, l2, l3 );

        // Smith's Method (1961)
        //double discriminant = eigenValsOf3by3SymMatrixSmith1961( m11, m12, m13, m22, m23, m33, ll1, ll2, ll3 );
        
        //cout << "discriminant = " << discriminant << endl;
        //cout << "l1: " << l1 << "  l2: " << l2 << "  l3: " << l3 << endl;

        // ********************************************************************
        // ASSUME all eigen values are real!

        // If discriminant > 0, then one root is real and two complex conjugate.
        /*
        if ( discriminant > errThreshold ) {
            cerr << "ERR:  Discriminant > 0, i.e. one eigen value is real and two complex conjugate!" << "\n";
            cerr << "CONT: The function is not intended for this kind of matrix!" << endl;
            return;
        }
        */



        // If all eigen values are equal, then eigen vectors can be
        // any vectors that are orthogonal to one another.
        // Also in this case, discriminant == 0.
        //if ( fabs(m12) < errThreshold && fabs(m13) < errThreshold && fabs(m23) < errThreshold ) {
        const double ratio = 100000.0;
        if ( fabs(l1 - l2) < fabs(l1/ratio) ) {
            if ( fabs(l1 - l3) < fabs(l1/ratio) ) {         // CASE: l1 == l2 == l3
                v1[0] = 1;  v1[1] = 0;  v1[2] = 0;
                v2[0] = 0;  v2[1] = 1;  v2[2] = 0;
                v3[0] = 0;  v3[1] = 0;  v3[2] = 1;
                cout << "CASE: l1 == l2 == l3\n";
            }
            else {                                          // CASE: l1 == l2 != l3
                FindAnEigenVector( m11, m12, m13, m22, m23, m33, l3, v3 );
                FindOrthogonalVectors( v3, v1, v2 );
                cout << "CASE: l1 == l2 != l3\n";
            }
        }
        else if ( fabs(l2 - l3) < fabs(l2/ratio) ) {        // CASE: l1 != l2 == l3
            FindAnEigenVector( m11, m12, m13, m22, m23, m33, l1, v1 );
            FindOrthogonalVectors( v1, v2, v3 );
            cout << "CASE: l1 != l2 == l3\n";
        }
        else if ( fabs(l1 - l3) < fabs(l1/ratio) ) {        // CASE: l1 == l3 != l2
            FindAnEigenVector( m11, m12, m13, m22, m23, m33, l2, v2 );
            FindOrthogonalVectors( v2, v1, v3 );
            cout << "CASE: l1 == l3 != l2\n";
        }
        else {                                              // CASE: l1 != l2 != l3
            //if ( fabs(l1) > fabs(l2) && fabs(l1) > fabs(l3) ) {}
            FindAnEigenVector( m11, m12, m13, m22, m23, m33, l1, v1 );
            //FindOrthogonalVectors( v1, v2, v3 );
            FindAnEigenVector( m11, m12, m13, m22, m23, m33, l2, v2 );
            FindAnEigenVector( m11, m12, m13, m22, m23, m33, l3, v3 );
            cout << "CASE: l1 != l2 != l3\n";
        }
        
        cout << "v1*v2 = " << DotProduct(v1,v2) << "\n";
        cout << "v2*v3 = " << DotProduct(v2,v3) << "\n";
        cout << "v1*v3 = " << DotProduct(v1,v3) << "\n";
    } // END FN:    eigenValsAndVecsOf3by3SymMatrix()

void SP_Maths::eigenValsAndVecsOf3by3SymMatrixUsingQRFactor	(	double	m11,
		double	m12,
		double	m13,
		double	m22,
		double	m23,
		double	m33,
		double &	l1,
		double &	l2,
		double &	l3,
		double	v1[3],
		double	v2[3],
		double	v3[3]
	)

Definition at line 603 of file sp_mathematics.cpp.

    { 
        //double errThreshold = 5E-3;   // error tolerant value
        static const int TIMES_LIMIT = 256;
        int iterationTimes = 1;

        SPtMatrix< double > IT( 3, 3 );
        SPtMatrix< double > EG( 3, 3 );
        IT.Set( 0, 0, m11 );    IT.Set( 0, 1, m12 );    IT.Set( 0, 2, m13 );
        IT.Set( 1, 0, m12 );    IT.Set( 1, 1, m22 );    IT.Set( 1, 2, m23 );
        IT.Set( 2, 0, m13 );    IT.Set( 2, 1, m23 );    IT.Set( 2, 2, m33 );
        
        //cout << "IT is \n";
        //IT.DisplayElements();
        do {
            iterationTimes *= 2;
        
            MatrixOperations::EigenValsByQRFactorization( IT, iterationTimes, &EG );

            cout << "QR Factorization iterates " << iterationTimes << " times.\n";

            //cout << "IT is \n";
            //IT.DisplayElements();
            //cout << "EG is \n";
            //EG.DisplayElements();

            l1 = EG.Get( 0, 0 );
            l2 = EG.Get( 1, 1 );
            l3 = EG.Get( 2, 2 );

            // If all eigen values are equal, then eigen vectors can be
            // any vectors that are orthogonal to one another.
            // Also in this case, discriminant == 0.
            //if ( fabs(m12) < errThreshold && fabs(m13) < errThreshold && fabs(m23) < errThreshold ) {
            const double ratio = 10000.0;
            if ( fabs(l1 - l2) < fabs(l1/ratio) ) {
                if ( fabs(l1 - l3) < fabs(l1/ratio) ) {         // CASE: l1 == l2 == l3
                    v1[0] = 1;  v1[1] = 0;  v1[2] = 0;
                    v2[0] = 0;  v2[1] = 1;  v2[2] = 0;
                    v3[0] = 0;  v3[1] = 0;  v3[2] = 1;
                    cout << "CASE: l1 == l2 == l3\n";
                }
                else {                                          // CASE: l1 == l2 != l3
                    FindAnEigenVector( m11, m12, m13, m22, m23, m33, l3, v3 );
                    FindOrthogonalVectors( v3, v1, v2 );
                    cout << "CASE: l1 == l2 != l3\n";
                }
            }
            else if ( fabs(l2 - l3) < fabs(l2/ratio) ) {        // CASE: l1 != l2 == l3
                FindAnEigenVector( m11, m12, m13, m22, m23, m33, l1, v1 );
                FindOrthogonalVectors( v1, v2, v3 );
                cout << "CASE: l1 != l2 == l3\n";
            }
            else if ( fabs(l1 - l3) < fabs(l1/ratio) ) {        // CASE: l1 == l3 != l2
                FindAnEigenVector( m11, m12, m13, m22, m23, m33, l2, v2 );
                FindOrthogonalVectors( v2, v1, v3 );
                cout << "CASE: l1 == l3 != l2\n";
            }
            else {                                              // CASE: l1 != l2 != l3
                //if ( fabs(l1) > fabs(l2) && fabs(l1) > fabs(l3) ) {}
                FindAnEigenVector( m11, m12, m13, m22, m23, m33, l1, v1 );
                //FindOrthogonalVectors( v1, v2, v3 );
                FindAnEigenVector( m11, m12, m13, m22, m23, m33, l2, v2 );
                FindAnEigenVector( m11, m12, m13, m22, m23, m33, l3, v3 );
                cout << "CASE: l1 != l2 != l3\n";
            }
            
            cout << "v1*v2 = " << DotProduct(v1,v2) << "\n";
            cout << "v2*v3 = " << DotProduct(v2,v3) << "\n";
            cout << "v1*v3 = " << DotProduct(v1,v3) << "\n";
        } while ( fabs(DotProduct(v1,v2) + DotProduct(v2,v3) + DotProduct(v1,v3)) > 1E-06 && iterationTimes <= TIMES_LIMIT );
    } // END FN:    eigenValsAndVecsOf3by3SymMatrixUsingQRFactor()

double SP_Maths::eigenValsOf3by3SymMatrix	(	double	m11,
		double	m12,
		double	m13,
		double	m22,
		double	m23,
		double	m33,
		double &	lamda1,
		double &	lamda2,
		double &	lamda3
	)

Definition at line 397 of file sp_mathematics.cpp.

    { 
        double errThreshold = 1E-6; // error tolerant value
        double b, c, d = 0;
        double f[5];

        // coefficients of the characteristic polynomial of the 3x3 symmetric matrix
        b = - m11 - m22 - m33;
        c = (m11*m22 + m22*m33 + m33*m11) - (m12*m12 + m23*m23 + m13*m13);
        //d = m11*m23*m23 + m22*m13*m13 + m33*m12*m12 - m11*m22*m33 - 2*m12*m23*m13;
        
        f[0] = m11*m23*m23;
        f[1] = m22*m13*m13;
        f[2] = m33*m12*m12;
        f[3] = -m11*m22*m33;
        f[4] = -2*m12*m23*m13;
        double dNeg = 0;
        for ( int i = 0; i < 5; i++ ) {
            if ( f[i] >= 0 )    d += f[i];
            else                dNeg += f[i];
        }
        d += dNeg;
        
        // the roots of the polynomial function are the eigen values of the matrix
        if ( fabs(m12) < errThreshold && fabs(m13) < errThreshold && fabs(m23) < errThreshold ) {
            double difThreshold = fabs( m11 * 1E-06 );
            if ( fabs(m11-m22) < difThreshold && fabs(m11-m33) < difThreshold ) {   // all roots are real and at least two are equal if discriminant = 0.
                lamda1 = lamda2 = lamda3 = ( m11 + m22 + m33 ) / 3.0;
                return  0;  // return zero discriminant value
            }
            else {                              // all roots are real and unequal if discriminant < 0.
                lamda1 = m11;
                lamda2 = m22;
                lamda3 = m33;
                return   -1;    // return negative discriminant value
            }
        }

        // Find roots using formulas from "Handbook of Mathematic"
        double discriminant = rootsOfCubicPolynomial( 1, b, c, d, lamda1, lamda2, lamda3 );

        return static_cast<double>(discriminant);
    } // END FN:    eigenValsOf3by3SymMatrix()

double SP_Maths::eigenValsOf3by3SymMatrixSmith1961	(	double	m11,
		double	m12,
		double	m13,
		double	m22,
		double	m23,
		double	m33,
		double &	lamda1,
		double &	lamda2,
		double &	lamda3
	)

Definition at line 449 of file sp_mathematics.cpp.

    { 

        double errThreshold = 1E-10;

        // 3m = trace(A)
        // 2q = det(A - mI)
        // 6p = the sum of squares of elements of (A - mI)
        // discriminant = p^3 - q^2
        // phi = 1/3*arctan( sqrt(discriminant)/q )
        // root1 = m + 2*sqrt(p)*cos(phi)
        // root2 = m - sqrt(p)*(cos(phi) + sqrt(3)*sin(phi))
        // root3 = m - sqrt(p)*(cos(phi) - sqrt(3)*sin(phi))
        double m, p, q, phi, discriminant;


        SPtMatrix3x3< double > IT(  m11, m12, m13,
                                    m12, m22, m23,
                                    m13, m23, m33 );
        SPtMatrix3x3< double > I(   1,0,0, 0,1,0, 0,0,1 );
        
        
        m = ( m11 + m22 + m33 ) / 3.0;
        
        SPtMatrix3x3< double > B = IT - m*I;
    
        q = (   B.Get(0,0) * ( B.Get(1,1)*B.Get(2,2) - B.Get(2,1)*B.Get(1,2) ) 
            -   B.Get(1,0) * ( B.Get(0,1)*B.Get(2,2) - B.Get(2,1)*B.Get(0,2) )
            +   B.Get(2,0) * ( B.Get(0,1)*B.Get(1,2) - B.Get(1,1)*B.Get(0,2) )
            ) / 2.0;

        p = 0;
        for ( int r = 0; r < 3; r++ ) {
            for ( int c = 0; c < 3; c++ ) {
                p += B.Get(r,c)*B.Get(r,c);
            }
        }
        p /= 6.0;

        discriminant = p*p*p - q*q;
        if ( discriminant >= 0 ) {
            phi = atan( sqrt(discriminant)/q ) / 3.0;

            lamda1 = m + 2*sqrt(p)*cos(phi);
            lamda2 = m - sqrt((double)p)*(cos(phi) + sqrt((double)3)*sin(phi));
            lamda3 = m - sqrt((double)p)*(cos(phi) - sqrt((double)3)*sin(phi));
        }
        else {
            cout << "Discriminant (Smith's Method (1961) < 0: discriminant = " << discriminant << "\n";
        }

        // the roots of the polynomial function are the eigen values of the matrix
        if ( fabs(m12) < errThreshold && fabs(m13) < errThreshold && fabs(m23) < errThreshold ) {
            double difThreshold = fabs( m11 * 1E-06 );
            if ( fabs(m11-m22) < difThreshold && fabs(m11-m33) < difThreshold ) {   // all roots are real and at least two are equal if discriminant = 0.
                lamda1 = lamda2 = lamda3 = ( m11 + m22 + m33 ) / 3.0;
                return  0;  // return zero discriminant value
            }
            else {                              // all roots are real and unequal if discriminant < 0.
                lamda1 = m11;
                lamda2 = m22;
                lamda3 = m33;
                return   -1;    // return negative discriminant value
            }
        }

        return static_cast<double>(discriminant);
    } // END FN:    eigenValsOf3by3SymMatrixSmith1961()

void SP_Maths::FindAnEigenVector	(	double &	m11,
		double &	m12,
		double &	m13,
		double &	m22,
		double &	m23,
		double &	m33,
		double	l1,
		double	v1[3]
	)

Definition at line 39 of file sp_mathematics.cpp.

                                                                                  : 3x3 Symmetric Matrix
                            double l1,      // O/P: eigen value
                            double v1[3]    // O/P: eigen vector
                          )
    {
        double a1,          b1 = m12,       c1 = m13;
        double a2 = m12,    b2,             c2 = m23;
        double a3 = m13,    b3 = m23,       c3;

        a1 = m11 - l1;
        b2 = m22 - l1;
        c3 = m33 - l1;
        double divider = b2*a1 - b1*a2 - b3*a1 + b1*a3;
        //if ( divider != 0 ) {
        if ( fabs(divider) > 1E-40 ) {
            v1[0] = (-b2*c1 + b1*c2 + b3*c1 - b1*c3) / divider;
            v1[1] = (a2*c1 - a1*c2 - a3*c1 + a1*c3) /divider;
            v1[2] = 1;

        // make the vector a unit vector
            double vecLen = sqrt( v1[0]*v1[0] + v1[1]*v1[1] + v1[2]*v1[2] );
            v1[0] /= vecLen;
            v1[1] /= vecLen;
            v1[2] /= vecLen;
        }
        else {
            v1[0] = 0;
            v1[1] = 0;
            v1[2] = 0;
        }
    } // END FN:    FindAnEigenVector()

void SP_Maths::FindOrthogonalVectors	(	const double	iA[3],
		double	oB[3],
		double	oC[3]
	)

Definition at line 76 of file sp_mathematics.cpp.

    {
        SPtMatrix< double > A( 3, 1 ), B( 3, 1 ), C( 3, 1 );
        for ( int i = 0; i < 3; i++ ) {
            A.Set( i, 0, iA[i] );
        }

        // alpha for x, beta for y, and zeta for z
        
        double angleZeta;
        if ( iA[0] != 0 )   angleZeta = atan( iA[1] / iA[0] );
        else                angleZeta = ( iA[1] > 0 ? 1.0 : -1.0 ) * PI / 2.0;
        double angleBeta = asin( iA[2] / sqrt( iA[0]*iA[0] + iA[1]*iA[1] + iA[2]*iA[2] ) );
        //double angleAlpha = 90.0/180.0*PI;
        SPtMatrix< double > Rz( 3, 3 );
        Rz.Set( 0, 0,  cos( angleZeta ) );  Rz.Set( 0, 1,  sin( angleZeta ) );
        Rz.Set( 1, 0, -sin( angleZeta ) );  Rz.Set( 1, 1,  cos( angleZeta ) );
        Rz.Set( 2, 2, 1 );
        SPtMatrix< double > Ry( 3, 3 );
        Ry.Set( 0, 0,  cos( angleBeta ) );  Ry.Set( 0, 2, -sin( angleBeta ) );
        Ry.Set( 1, 1, 1 );
        Ry.Set( 2, 0,  sin( angleBeta ) );  Ry.Set( 2, 2,  cos( angleBeta ) );
        
        SPtMatrix< double > Ryy( 3, 3 );
        Ryy.Set( 0, 0,  cos( PI/2.0 ) );    Ryy.Set( 0, 2, -sin( PI/2.0 ) );
        Ryy.Set( 1, 1, 1 );
        Ryy.Set( 2, 0,  sin( PI/2.0 ) );    Ryy.Set( 2, 2,  cos( PI/2.0 ) );

        /*
        SPtMatrix< double > Rx( 3, 3 );
        Rx.Set( 0, 0, 1 );
        Rx.Set( 1, 1,  cos( angleAlpha ) ); Rx.Set( 1, 2,  sin( angleAlpha ) );
        Rx.Set( 2, 1, -sin( angleAlpha ) ); Rx.Set( 2, 2,  cos( angleAlpha ) );
        */
        
        SPtMatrix< double > Rym( 3, 3 );
        Rym.Set( 0, 0,  cos( -angleBeta ) );    Rym.Set( 0, 2, -sin( -angleBeta ) );
        Rym.Set( 1, 1, 1 );
        Rym.Set( 2, 0,  sin( -angleBeta ) );    Rym.Set( 2, 2,  cos( -angleBeta ) );
        SPtMatrix< double > Rzm( 3, 3 );
        Rzm.Set( 0, 0,  cos( -angleZeta ) );    Rzm.Set( 0, 1,  sin( -angleZeta ) );
        Rzm.Set( 1, 0, -sin( -angleZeta ) );    Rzm.Set( 1, 1,  cos( -angleZeta ) );
        Rzm.Set( 2, 2, 1 );

        //B = Rz * Ry * Ryy * Rym * Rzm * A;
        B = Rzm * Rym * Ryy * Ry * Rz * A;

        //B.DisplayElements();
        //(A.GetTranspose()*B).DisplayElements();

        oB[0] = B.Get( 0, 0 );
        oB[1] = B.Get( 1, 0 );
        oB[2] = B.Get( 2, 0 );

        cout << "\niA*oB = " << ( iA[0]*oB[0] + iA[1]*oB[1] + iA[2]*oB[2] ) << "\n";
        CrossProduct( iA, oB, oC );
    } // END FN:    FindOrthogonalVectors()

double SP_Maths::RootFindingByNewtonRaphsonMethod	(	void(*)(double x, double *fx, double *dfx)	funcd,
		double	x1,
		double	x2,
		double	xacc
	)

Definition at line 688 of file sp_mathematics.cpp.

    {
        const int ITERATION_LIMIT = 200;
        double x, fx, dx = 0, dfx;

        x = 0.5 * (x1 + x2);    // initial guess is the middle of the interval (x1,x2)
        for ( int i = 1; i < ITERATION_LIMIT; i++ ) {
            (*funcd)( x, &fx, &dfx );
            if ( fx != 0 && dfx != 0) {
                dx = fx/dfx;
                x -= dx;
            }
            if ( (x1-x)*(x-x2) < 0.0 ) {
                cerr << "From RootFindingByNewtonRaphsonMethod():\n"
                     << "  The root lies outside the provided interval (x1,x2)!" << endl;
                return -777;
            }
            if ( fabs(dx) < xacc )  return x;       // the convergence value is the root
        }

        cerr << "From RootFindingByNewtonRaphsonMethod():\n"
             << "  Maximum number of iterations exceeded!" << endl;
        return -777;
    } // END FN:    RootFindingByNewtonRaphsonMethod()

double SP_Maths::rootOfCubicPolynomial	(	double	a,
		double	b,
		double	c,
		double	d,
		int	rootNo
	)

Definition at line 306 of file sp_mathematics.cpp.

    {
        double a1, a2, a3;
        double Q, R, D; 
        double S, T;    // FIX IT HERE!
        double zeta;
        double rootValue;

        a1 = b/a;
        a2 = c/a;
        a3 = d/a;
        Q = (3.0*a2 - a1*a1) / 9.0;
        R = (9.0*a1*a2 - 27.0*a3 - 2.0*a1*a1*a1) / 54.0;
        D = Q*Q*Q + R*R;

        // If a1, a2, a3 are real and if D = Q^3 + R^2 is the discriminant, then
        // (i)   one root is real and two complex conjugate if D > 0
        // (ii)  all roots are real and at least two are equal if D = 0
        // (iii) all roots are real and unequal if D < 0.

        if ( D >= 0 ) {
            double Dh = sqrt(D);
            S = pow(R + Dh, 1.0/3.0);   // FIX IT HERE!
            T = R - Dh;
            if ( T >= 0 )   T =  pow(  T, 1.0/3.0 );
            else            T = -pow( -T, 1.0/3.0 );
            
            switch ( rootNo ) {
                // first root
                case 1:
                    rootValue = S + T - a1/3.0;
                    return rootValue;
                    break;

                // second root
                case 2:
                    rootValue = -(S + T)/2.0 - a1/3.0;  // + i*sqrt(3.0)*(S - T)/2.0;
                    return rootValue;
                    break;

                // third root
                case 3:
                    rootValue = -(S + T)/2.0 - a1/3.0;  // - i*sqrt(3.0)*(S - T)/2.0;
                    return rootValue;
                    break;
                
                // default
                default:
                    cerr << "ERROR:  Incorrect Root# Option!" << endl;
                    return -777;
                    break;
            }
        }

        // All roots are real and unequal if D < 0.
        //if ( D < 0 ) {
        else {
            zeta = acos( R / sqrt(-(Q*Q*Q)) );
            switch ( rootNo ) {
                // first root
                case 1:
                    rootValue = 2 * sqrt(-Q) * cos(zeta/3.0) - a1/3.0;
                    return rootValue;
                    break;

                // second root
                case 2:
                    rootValue = 2 * sqrt(-Q) * cos(zeta/3.0 + 2.0/3.0*PI) - a1/3.0;
                    return rootValue;
                    break;

                // third root
                case 3:
                    rootValue = 2 * sqrt(-Q) * cos(zeta/3.0 + 4.0/3.0*PI) - a1/3.0;
                    return rootValue;
                    break;
                
                // default
                default:
                    cerr << "ERROR:  Incorrect Root# Option!" << endl;
                    return -777;
                    break;
            }
        }
    } // END FN:    rootOfCubicPolynomial()

double SP_Maths::rootsOfCubicPolynomial	(	double	a,
		double	b,
		double	c,
		double	d,
		double &	r1,
		double &	r2,
		double &	r3
	)

Definition at line 284 of file sp_mathematics.cpp.

    {
        std::complex< double > cr1, cr2, cr3;
        double discriminant = rootsOfCubicPolynomial( a, b, c, d, cr1, cr2, cr3 );
        r1 = cr1.real();
        r2 = cr2.real();
        r3 = cr3.real();

        return discriminant;
    }

double SP_Maths::rootsOfCubicPolynomial	(	double	a,
		double	b,
		double	c,
		double	d,
		std::complex< double > &	r1,
		std::complex< double > &	r2,
		std::complex< double > &	r3
	)

Definition at line 194 of file sp_mathematics.cpp.

    {
        double a1, a2, a3;
        double Q, R, D;
        std::complex<double> S, T;
        double zeta;

        a1 = b/a;
        a2 = c/a;
        a3 = d/a;

        //Q = (3.0*a2 - a1*a1) / 9.0;
        //R = (9.0*a1*a2 - 27.0*a3 - 2.0*a1*a1*a1) / 54.0;
        // Q^3 + R^2 can factor 9*9*9 (or 27*27) out
        // So  Q = 3*a2 - a1^2
        // and R = (f1 + f2 + f3) / 2

        Q = 3.0*a2 - a1*a1;
        double Rpos = 0, Rneg = 0;
        double f[3];
        f[0] = 9.0*a1*a2;
        f[1] = -27.0*a3;
        f[2] = -2.0*a1*a1*a1;
        for ( int i = 0; i < 3; i++ ) {
            ( f[i] >= 0 ) ? ( Rpos += f[i] ) : ( Rneg += f[i] );
            //if ( f[i] >= 0 )  Rpos += f[i];
            //else              Rneg += f[i];
        }
        R = (Rpos + Rneg) / 2.0;



        D = ( Q*Q*Q + R*R ) / 729.0;
        Q /= 9.0;
        R /= 27.0;

        // If a1, a2, a3 are real and if D = Q^3 + R^2 is the discriminant, then
        // (i)   one root is real and two complex conjugate if D > 0
        // (ii)  all roots are real and at least two are equal if D = 0
        // (iii) all roots are real and unequal if D < 0.

        if ( D >= 0 ) {
            std::complex<double> Dh = sqrt(D);
            
            S = pow(R + Dh, 1.0/3.0);
            T = pow(R - Dh, 1.0/3.0);

            r1 = S + T - a1/3.0;
            std::complex<double> pureComplex( 0, sqrt(3.0) );
            //r2 = -(S + T)/2.0 - a1/3.0 + i*sqrt(3.0)*(S - T)/2.0;
            //r3 = -(S + T)/2.0 - a1/3.0 - i*sqrt(3.0)*(S - T)/2.0;
            r2 = - (S + T)/2.0 - a1/3.0 + pureComplex*(S - T)/2.0;
            r3 = - (S + T)/2.0 - a1/3.0 - pureComplex*(S - T)/2.0;
        }

        // All roots are real and unequal if D < 0.
        //if ( D < 0 ) {
        else {
            //cout << "D < 0" << endl;
            zeta = acos( R / sqrt(-(Q*Q*Q)) );
            r1 = 2 * sqrt(-Q) * cos(zeta/3.0) - a1/3.0;
            r2 = 2 * sqrt(-Q) * cos(zeta/3.0 + 2.0/3.0*PI) - a1/3.0;
            r3 = 2 * sqrt(-Q) * cos(zeta/3.0 + 4.0/3.0*PI) - a1/3.0;
        }

        // DEBUG:
        cout.precision( 10 );
        cout    << "*x^3"  << (a1>0 ? " + " : " - ") << fabs(a1) << "*x^2"<< (a1>0 ? " + " : " - ") 
                << fabs(a2) << "*x" << (a1>0 ? " + " : " - ") << fabs(a3) << "\n";
        cout << "Q = " << Q << "\n";
        cout << "R = " << R << "\n";
        cout << "Q^3 = " << Q*Q*Q << "\n";
        cout << "R^2 = " << R*R << "\n";
        cout << "D = " << D << "\n";
        cout << "Root#1: " << r1 << " Root#2: " << r2 << " Root#3: " << r3 << "\n";
        cout.precision( 6 );


        std::complex<double> cr1, cr2, cr3;
        rootsOfCubicPolynomialCodeGeneratedByMaple( a, b, c, d, cr1, cr2, cr3 );


        return D;   // return discriminant
    } // END FN:    rootsOfCubicPolynomial()

void SP_Maths::rootsOfCubicPolynomialCodeGeneratedByMaple	(	double	a,
		double	b,
		double	c,
		double	d,
		std::complex< double > &	r1,
		std::complex< double > &	r2,
		std::complex< double > &	r3
	)

Definition at line 135 of file sp_mathematics.cpp.

    {
        std::complex<double> a1;
        std::complex<double> a2;
        std::complex<double> a3;

        a1 = b/a;
        a2 = c/a;
        a3 = d/a;

        std::complex<double> div = 1.0/3.0;

        // common part for all roots
        std::complex<double> t1 = a2*a1;
        std::complex<double> t4 = a1*a1;
        std::complex<double> t5 = t4*a1;
        std::complex<double> t7 = a2*a2;
        std::complex<double> t14 = a3*a3;
        std::complex<double> t19 = sqrt(12.0*t7*a2-3.0*t7*t4-54.0*t1*a3+81.0*t14+12.0*a3*t5);
        std::complex<double> t22 = pow(36.0*t1-108.0*a3-8.0*t5+12.0*t19,div);

        // Root# 1
        r1 = t22/6.0-6.0*(a2/3.0-t4/9.0)/t22-a1/3.0;
        
        // common part for root# 2 and 3
        std::complex<double> t28 = (a2/3.0-t4/9.0)/t22;
        std::complex<double> t31 = sqrt(3.0);

        // Root# 2
        //r2 = -t22/12.0+3.0*t28-a1/3.0+(1.0/2.0*sqrt(-1.0))*t31*(t22/6.0+6.0*t28);
        std::complex<double> st1 = t31*(t22/6.0+6.0*t28);
        std::complex<double> st2 = (sqrt(-1.0));
        cout << st1 << ", " << st2 << endl; 
        r2 = -t22/12.0+3.0*t28-a1/3.0+(1.0/2.0*sqrt(-1.0))*t31*(t22/6.0+6.0*t28);

        // Root# 3
        //r3 = -t22/12.0+3.0*t28-a1/3.0+(-1.0/2.0*sqrt(-1.0))*t31*(t22/6.0+6.0*t28);
        r2 = -t22/12.0+3.0*t28-a1/3.0+(1.0/2.0*sqrt(-1.0))*t31*(t22/6.0+6.0*t28);


        // DEBUG:
        cout.precision( 10 );
//      cout    << "*x^3"  << (a1>0 ? " + " : " - ") << fabs(a1) << "*x^2"<< (a1>0 ? " + " : " - ") 
//              << fabs(a2) << "*x" << (a1>0 ? " + " : " - ") << fabs(a3) << "\n";
        cout << "Root#1: " << r1 << " Root#2: " << r2 << " Root#3: " << r3 << "\n";
        cout.precision( 6 );

    } // END FN:    rootsOfCubicPolynomial()

---

Variable Documentation

const double SP_Maths::PI = 3.141592654

Definition at line 27 of file sp_mathematics.h.