 |
TAPs 0.7.7.3
|
|
TAPs 0.7.7.3
|
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 | 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]) |
| void | rootsOfCubicPolynomialCodeGeneratedByMaple (double a, double b, double c, double d, std::complex< double > &r1, std::complex< double > &r2, std::complex< double > &r3) |
| double | rootsOfCubicPolynomial (double a, double b, double c, double d, std::complex< double > &r1, std::complex< double > &r2, std::complex< double > &r3) |
| double | rootsOfCubicPolynomial (double a, double b, double c, double d, double &r1, double &r2, double &r3) |
| double | rootOfCubicPolynomial (double a, double b, double c, double d, int rootNo) |
| 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 | 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 | RootFindingByNewtonRaphsonMethod (void(*funcd)(double x, double *fx, double *dfx), double x1, double x2, double xacc) |
Variables | |
| const double | PI = 3.141592654 |
| void SP_Maths::CrossProduct | ( | const double | A[3], |
| const double | B[3], | ||
| double | C[3] | ||
| ) |
Definition at line 20 of file sp_mathematics.cpp.
Referenced by FindOrthogonalVectors().
{
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()
Here is the caller graph for this function:| double SP_Maths::DotProduct | ( | const double | A[3], |
| const double | B[3] | ||
| ) |
Definition at line 30 of file sp_mathematics.cpp.
Referenced by eigenValsAndVecsOf3by3SymMatrix(), and eigenValsAndVecsOf3by3SymMatrixUsingQRFactor().
{
return A[0]*B[0] + A[1]*B[1] + A[2]*B[2];
} // END FN: DotProduct()
Here is the caller graph for this function:| 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 525 of file sp_mathematics.cpp.
References DotProduct(), eigenValsOf3by3SymMatrix(), FindAnEigenVector(), and FindOrthogonalVectors().
{
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()
Here is the call graph for this function:| 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 602 of file sp_mathematics.cpp.
References DotProduct(), MatrixOperations::EigenValsByQRFactorization(), FindAnEigenVector(), FindOrthogonalVectors(), SPtMatrix< T >::Get(), and SPtMatrix< T >::Set().
{
//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()
Here is the call graph for this function:| 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 396 of file sp_mathematics.cpp.
References rootsOfCubicPolynomial().
Referenced by eigenValsAndVecsOf3by3SymMatrix().
{
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()
Here is the call graph for this function: Here is the caller graph for this function:| 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 448 of file sp_mathematics.cpp.
References SPtMatrix3x3< T >::Get().
{
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()
Here is the call graph for this function:| void SP_Maths::FindAnEigenVector | ( | double & | m11, |
| double & | m12, | ||
| double & | m13, | ||
| double & | m22, | ||
| double & | m23, | ||
| double & | m33, | ||
| double | l1, | ||
| double | v1[3] | ||
| ) |
Definition at line 38 of file sp_mathematics.cpp.
Referenced by eigenValsAndVecsOf3by3SymMatrix(), eigenValsAndVecsOf3by3SymMatrixUsingQRFactor(), and GenMath< T >::EigenValsAndVecsOf3by3SymMatrixUsingQRFactor().
{
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()
Here is the caller graph for this function:| void SP_Maths::FindOrthogonalVectors | ( | const double | iA[3], |
| double | oB[3], | ||
| double | oC[3] | ||
| ) |
Definition at line 75 of file sp_mathematics.cpp.
References CrossProduct(), SPtMatrix< T >::Get(), PI, and SPtMatrix< T >::Set().
Referenced by eigenValsAndVecsOf3by3SymMatrix(), eigenValsAndVecsOf3by3SymMatrixUsingQRFactor(), and GenMath< T >::EigenValsAndVecsOf3by3SymMatrixUsingQRFactor().
{
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()
Here is the call graph for this function: Here is the caller graph for this function:| double SP_Maths::RootFindingByNewtonRaphsonMethod | ( | void(*)(double x, double *fx, double *dfx) | funcd, |
| double | x1, | ||
| double | x2, | ||
| double | xacc | ||
| ) |
Definition at line 687 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 305 of file sp_mathematics.cpp.
References PI.
Referenced by VolPresTriModel< T >::PreserveVolume(), VolPresPolygonalModel< T >::PreserveVolume(), VolPresPolygonalModel< T >::PreserveVolumeByVertexNormals(), and ExtendedOpenGLPNTriangleVolPresModel< T >::PreserveVolumeByVertexNormals().
{
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()
Here is the caller graph for this function:| double SP_Maths::rootsOfCubicPolynomial | ( | double | a, |
| double | b, | ||
| double | c, | ||
| double | d, | ||
| double & | r1, | ||
| double & | r2, | ||
| double & | r3 | ||
| ) |
Definition at line 283 of file sp_mathematics.cpp.
References rootsOfCubicPolynomial().
{
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;
}
Here is the call graph for this function:| 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 193 of file sp_mathematics.cpp.
References PI, and rootsOfCubicPolynomialCodeGeneratedByMaple().
Referenced by eigenValsOf3by3SymMatrix(), and rootsOfCubicPolynomial().
{
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()
Here is the call graph for this function: Here is the caller graph for this function:| 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 134 of file sp_mathematics.cpp.
Referenced by rootsOfCubicPolynomial().
{
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()
Here is the caller graph for this function:| const double SP_Maths::PI = 3.141592654 |
Definition at line 27 of file sp_mathematics.h.
Referenced by FindOrthogonalVectors(), rootOfCubicPolynomial(), and rootsOfCubicPolynomial().
1.7.4 