Function findCanonicalForm

plugin/seq/plotPDF.cpp:1800–1872

    GaussElimination(phi, a, NQ);

    delete [] a[0];
    
    return;
}

void findCanonicalForm( double *const PHI, const double *const phi )
{
    const double &a = phi[0]; const double &b = phi[1]; const double &c = phi[2];
    const double &d = phi[3]; const double &e = phi[4]; const double &f = phi[5];
    
    // phi(x,y) = a*x*x + b*x*y + c*y*y + d*x + e*y + f
    //          = [x,y][ a, b/2; b/2, c ][x;y] + [d,e][x;y] + f
    //          = [x,y]P P^T [ a, b/2; b/2, c ]P P^T[x;y] + [d,e]P P^T[x;y] + f
    //          = [X,Y] [lambda1,0;0,lambda2] [X;Y] + [D,E][X;Y] + f
    //          = lambda1*X*X + lambda2*Y*Y + D*X + E*Y + f
    //          = lambda1 * (X + D/(2*lambda1))^2 + lambda2 * (Y + E/(2*lambda2))^2
    //              + ( -D*D/(4*lambda1) - E*E/(4*lambda2) + f )
    //          = PHI(X,Y)
    // where
    //   v1 = [v1x;v1y] (resp. v2 = [v2x;v2y]) is an eigenvector of lambda1 (resp. lambda2)
    //   P = [v1x, v2x; v1y, v2y],  satisfying P^T P = I <=> || v1 || = || v2 || = 1, v1 \perp v2 = 0
    //   [X;Y] = P^T[x;y] <=> [x;y] = P[X;Y]
    //   [D;E] = P^T[d;e] <=> [d;e] = P[D;E]

    const double det   = (a-c)*(a-c)+b*b;
    const double sqdet = sqrt( det );

    // eigenvalues of [ a, b /2; b/2, c]
    double &lambda1 = PHI[0]; double &lambda2 = PHI[1];
    lambda1 = ( (a+c) + sqdet ) / 2;
    lambda2 = ( (a+c) - sqdet ) / 2;

    // corresponding eivenvectors
    double &v1x = PHI[2]; double &v1y= PHI[3];
    double &v2x = PHI[4]; double &v2y= PHI[5];

    if( a < c ){

	const double n = sqrt( 2*det - 2*(a-c)*sqdet );
	v1x = b / n;
	v1y = (-(a-c)+sqdet) / n;
	v2x = (a-c-sqdet) / n;
	v2y = b / n;

    } else if( a > c ){

	const double n = sqrt( 2*det + 2*(a-c)*sqdet );
	v1x = (a-c+sqdet) / n;
	v1y = b / n;
	v2x = b / n;
	v2y = (-(a-c)-sqdet) / n;

    } else { // a == c

	lambda1 = (2*a+b)/2;
	lambda2 = (2*a-b)/2;
	v1x = v1y = v2x = v2y = 1/sqrt(static_cast<double>(2));  v2y *= -1;
    }