Can't figure out how to calculate 4th degree polynomials

I'm attempting to create a function that calculates the 4 roots of a 4th degree polynomial with included complex numbers. In my search for a formula, I came across a rather simple one contained in this discussion, described by Tito Piezas III towards the bottom of the page.

Now, I believe that the real error isn't really fault in my code (as I'm sure will be annoying to proofread) but in my understanding of the method involved. My problem is, that the quadratic roots are complex, and I don't know how to use the complex numbers when calculating the quartic roots, programatically.

He suggests deriving the quartic roots using the roots of two quadratic equations. I attempted to mimic the formula the best I can with my code below. The idea is that I calculate the two quadratic roots (under the premise they are positive only --- I don't know how otherwise), then, using those results, I can calculate the qurtic roots which I then save the real and complex values into x1,x2,x3,x4 into r1,r2,r3,r4,c1,c2,c3,c4 respectively. But, when calculating the quadratic roots, u, a value used later to calculate quartic roots: is complex!

Here is an image of his formula and steps. Blow is my code with captions on most steps.

double a, b, c, d;
double c1, c2, c3, c4;      //complex values
double r1, r2, r3, r4;      //real values

//      x^4+ax^3+bx^2+cx+d=0
a = 3;     
b = 4; 
c = 5;     //<--- example coefficients
d = 6;
if (a != 0) {

    double u,v1,v2;       
    double x,y,z;        //essentially a,b,c that he uses

    x=1;
    y= -2*b*b*b+9*a*b*c-27*c*c-27*a*a*d+72*b*d;
    z= Math.pow((b*b-3*a*c+12*d),3);

      //calculation of the v roots
    v1 = -y+(Math.sqrt(y*y-4*x*z))/(2*x);  // < negative root
    v2 = -y-(Math.sqrt(y*y-4*x*z))/(2*x);  // < negative root

//---calculations after this are invalid since v1 and v2 are NaN---       

    u = (a*a)/4 + ((-2*b+Math.pow(v1,1/3)+Math.pow(v2,1/3))/3);

    double x12sub,x34sub;

    x12sub= 3*a*a-8*b-4*u+((-a*a*a+4*a*b-8*c)/(Math.sqrt(u)));
    x34sub= 3*a*a-8*b-4*u-((-a*a*a+4*a*b-8*c)/(Math.sqrt(u)));

    r1 = -(1/4)*a +(1/2)*(Math.sqrt(u));
    r2 = -(1/4)*a +(1/2)*(Math.sqrt(u));
    r3 = -(1/4)*a -(1/2)*(Math.sqrt(u));
    r4 = -(1/4)*a -(1/2)*(Math.sqrt(u));

//--casting results into their orderly variables--

    if(x12sub<0){
        x12sub= x12sub*-1;
        x12sub = Math.sqrt(x12sub);
        x12sub = x12sub*(1/4);
        c1=x12sub;
        c2=x12sub;
    }
    else{
        r1=r1+x12sub;
        r2=r2-x12sub;
    }
    if(x34sub<0){
        x34sub= x34sub*-1;
        x34sub = Math.sqrt(x34sub);
        x34sub = x34sub*(1/4);
        c3=x34sub;
        c4=x34sub;  
    }
    else{
        r3=r3+x34sub;
        r4=r4+x34sub;
    }
}

I'm open for ANY solution. Even ones which involve the use of libraries that could help me. Thanks for the help.

Solution

Try using the Efficient Java Matrix Library. You can download the jars here: https://sourceforge.net/projects/ejml/files/v0.28/

You need to have this method in your class:

public static Complex64F[] findRoots(double... coefficients) {
    int N = coefficients.length-1;

    // Construct the companion matrix
    DenseMatrix64F c = new DenseMatrix64F(N,N);

    double a = coefficients[N];
    for( int i = 0; i < N; i++ ) {
        c.set(i,N-1,-coefficients[i]/a);
    }
    for( int i = 1; i < N; i++ ) {
        c.set(i,i-1,1);
    }

    // use generalized eigenvalue decomposition to find the roots
    EigenDecomposition<DenseMatrix64F> evd =  DecompositionFactory.eig(N,false);

    evd.decompose(c);

    Complex64F[] roots = new Complex64F[N];

    for( int i = 0; i < N; i++ ) {
        roots[i] = evd.getEigenvalue(i);
    }

    return roots;
}

Then you can use this to find, for example, the roots of x^2 + 4x + 4:

Complex64F[] c = findRoots(4, 4, 1);
    for(Complex64F f : c)
        System.out.println(f.toString());

This will print out:

-2
-2

Which is the desired result.