Cubic-solver.cpp

//Phaw! This certainly needs cleaning up, but it does work.

#include <complex>
#include <iostream>
#include <cmath>

using namespace std;

typedef long double R;
typedef complex<double> cplex;


void Gaussian_Elimination(cplex a00,cplex a01,cplex a02,
                          cplex a10,cplex a11,cplex a12,
                          cplex a20,cplex a21,cplex a22, cplex e1, cplex e2, cplex e3) {

}

struct real_cubic {
    R a,b,c,d;
    real_cubic(R _a,R _b,R _c,R _d) 
        : a(_a),b(_b),c(_c),d(_d){
    }
    R operator()(R x) const {
        return a*x*x*x + b*x*x + c*x + d;
    }
    R solve() const {
        R delta2=(b*b - 3*a*c)/(9*a*a);
        R h = -2*a*sqrt(delta2)*sqrt(delta2)*sqrt(delta2);
        R Xn = -b/(3*a);
        R Yn = (*this)(Xn);

        cout << "delta2=" << delta2 << endl;
        cout << "h=" << h << endl;
        cout << "Xn=" << Xn <<" Yn=" << Yn << endl;

        if(Yn*Yn < h*h) {
            cout << "Yn**2 < h**2 => Three distinct real roots" << endl;
            R theta=acos(Yn/h)/3;
            cout << "theta = " << theta/M_PI/2*360 << " degrees" << endl;

            R alpha=Xn + 2 * sqrt(delta2) * cos(theta);
            R beta =Xn + 2 * sqrt(delta2) * cos(theta + 2*M_PI/3);
            R gamma=Xn + 2 * sqrt(delta2) * cos(theta + 4*M_PI/3);

            cout << "(alpha,beta,gamma) = (" << alpha << "," << beta << "," << gamma << ")" << endl;

            cout << "f(alpha)=" << (*this)(alpha) << endl;
            cout << "f(beta)=" << (*this)(beta) << endl;
            cout << "f(gamma)=" << (*this)(gamma) << endl;
        } else {
            cout << "n**2 > h**2 => Only one real root";
        }

        return 0;
    }
};



template<typename F>
complex<F> acos(complex<F> z) {
    return log(z + complex<F>(0,1) * sqrt(complex<F>(1,0)-z*z))/complex<F>(0,1);
}

template<typename C>
struct cubic {
    C a,b,c,d;
    cubic(C _a,C _b,C _c,C _d) 
        : a(_a),b(_b),c(_c),d(_d){
    }
    C operator()(C x) const {
        return a*x*x*x + b*x*x + c*x + d;
    }
    C solve(C* lambda1,C* lambda2,C* lambda3) const {
        C delta2=(b*b - 3.0*a*c)/(9.0*a*a);
        C h = -2.0*a*sqrt(delta2)*sqrt(delta2)*sqrt(delta2);
        C Xn = -b/(3.0*a);
        C Yn = (*this)(Xn);

        cout << "delta2=" << delta2 << endl;
        cout << "h=" << h << endl;
        cout << "Xn=" << Xn <<" Yn=" << Yn << endl;

        cout << "Yn/h=" << Yn/h << endl;

        //Cind the complex arccosine

        C theta=acos(Yn/h)/C(3.0,0.0);
        cout << "theta = " << theta  << endl;

        C alpha=Xn + C(2.0,0.0) * sqrt(delta2) * cos(theta);
        C beta =Xn + C(2.0,0.0) * sqrt(delta2) * cos(theta + C(2*M_PI/3,0));
        C gamma=Xn + C(2.0,0.0) * sqrt(delta2) * cos(theta + C(4*M_PI/3,0));

        cout << "f(alpha)=" << (*this)(alpha) << endl;
        cout << "f(beta)=" << (*this)(beta) << endl;
        cout << "f(gamma)=" << (*this)(gamma) << endl;

        *lambda1 = alpha;
        *lambda2 = beta;
        *lambda3 = gamma;

        return 0;
    }
};


int main() {
    real_cubic f(1.0,7,14,8);

    cout << "f(5)=" << f(5) << endl;
    f.solve();

    {
        cubic<cplex> g(cplex( 3.0,1.0),
                       cplex( 7.0,1.0),
                       cplex(14.0,1.0),
                       cplex( 8.0,4.0));
        cplex lambda1,lambda2,lambda3;
        g.solve(&lambda1,&lambda2,&lambda3);
        cout << endl;
        cout << "Solution1 = " << lambda1 << endl;
        cout << "Solution2 = " << lambda2 << endl;
        cout << "Solution3 = " << lambda3 << endl;
    }
    {
        //3x3 matrix terms
        cplex
            a = cplex(1.0,0.0) ,b = cplex(0.3,0.0) ,c = cplex(0.9,0.0) ,
            d = cplex(0.3,0.0) ,e = cplex(2.0,0.0) ,f = cplex(0.4,0.0) ,
            g = cplex(0.8,0.0) ,h = cplex(0.2,0.0) ,i = cplex(3.0,0.0) ;
    
        cplex poly_a = cplex(-1,0);
        cplex poly_b = a+e+i;
        cplex poly_c = d*b+g*c+f*h-a*e-a*i-e*i;
        cplex poly_d = a*e*i - a*f*h - d*b*i + d*c*h + g*b*f - g*c*e;
    
        cubic<cplex> characteristic_polynomial(poly_a,poly_b,poly_c,poly_d);
        cplex lambda1,lambda2,lambda3;
        characteristic_polynomial.solve(&lambda1,&lambda2,&lambda3);
    
        cout << endl;
        cout << "Eigenvalue1 = " << lambda1 << endl;
        cout << "Eigenvalue2 = " << lambda2 << endl;
        cout << "Eigenvalue3 = " << lambda3 << endl;    


        //Assume the first element of the eigenvector is 1;
        cplex x1 = (f*(lambda1-a) + a*c)/(c*(lambda1 - e) + b*f);
        cplex x2 = (-a-(e-lambda1)*x1)/f;
        
        cout << endl;
        cout << "First Eigenvector: (" << 1 << "," << x1 << "," << x2 << ")" << endl;

        cplex 
            t1 = a+b*x1+c*x2,
            t2 = d+e*x1+f*x2,
            t3 = g+h*x1+i*x2; 

        cout << "First Eigenvector tranformed: (" << t1/lambda1 << "," << t2/lambda1 << "," << t3/lambda1 << ")" << endl;

    }


    return 0;
}