MATH_CubicPolynom.hpp

Go to the documentation of this file.

#ifndef MATH_CubicPolynom_HPP
#define MATH_CubicPolynom_HPP

#include "MATH_QuadraticPolynom.h"

template<class T>
void MATH_CubicPolynom::solve(const tReal& a,const tReal& b,const tReal& c,const tReal& d,
                              int& n, T solutions[3])  {
   
    //ax^3+bx^2+cx+d=0
    if (a==0) {
        MATH_QuadraticPolynom::solve<T>(b,c,d,n,solutions);
        return ;
    }

    //X^3+pX+q=0
    //X=x-dx
    tReal dx=-b/(3*a);
    tReal p=c/a-b*b/(3*a*a);
    tReal q=(b/(27*a))*(2*b*b/(a*a)-9*c/a)+d/a;
    
    tReal delta=-(4*p*p*p+27*q*q);
   
    if (delta<0) {
        delta=sqrt(-delta/27);
        tReal u=(-q+delta)/2;
        if (u>0) u=pow(u,1./3.);
        else u=-pow(-u,1./3.);
        
        tReal v=(-q-delta)/2;
        if (v>0) v=pow(v,1./3.);
        else v=-pow(-v,1./3.);
        n=1;
        solutions[0]=(T)(u+v+dx);
        
    } else if (delta>0) {
        n=3;
        tReal theta=acos((-q/2)*sqrt(27/(-p*p*p)))/3.;
        tReal r=2*sqrt(-p/3);
        for (int k=0;k<3;k++) {
            solutions[k]=(T) (r*cos(theta+2*k*M_PI/3)+dx);
        }
        if (solutions[0]>solutions[1]) std::swap(solutions[0],solutions[1]);
        if (solutions[1]>solutions[2]) {
            std::swap(solutions[1],solutions[2]);
            if (solutions[0]>solutions[1]) std::swap(solutions[0],solutions[1]);
        }
        
    } else {
        if ((p==0) && (q==0)) {
            n=1;
            solutions[0]=(T) dx;
        } else {
            n=2;
            solutions[0]=(T) dx+3*q/p;
            solutions[1]=(T) dx-1.5*q/p;
            if (solutions[0]>solutions[1]) std::swap(solutions[0],solutions[1]);
        }
    }
    
}

#endif