5

I have tried implementing a method that can solve a quartic polynomial given a, b, c, d, e, using this method: https://math.stackexchange.com/a/786/127747

It works for some solutions were there are 1 or 2 real roots, but the problem is that, sometimes the square or cubic roots involved might cause NaN values to appear in the intermediate variables if they take negatives as input, for example Math.sqrt(-9), which then messes with the final answer, making all of the roots NaN in the end of the method.

Is there any fast analytical way to only get all real roots of a quartic polynomial in Java, given variables/coefficients a, b, c, d and e, which does not involve some Complex library etc?

Edit: (Any understandable language works, but preferably Java, and if that's not the case, I will anyways make a port, and edit the answer to append it)

Edit 2: Here is my current code, where s is p from the equation, and q are just variables to optimize it a little bit, so that the same calculations aren't done twice:

public static double[] solveRealQuarticRoots(double a, double b, double c, double d, double e) {
    double s1 = 2 * c * c * c - 9 * b * c * d + 27 * (a * d * d + b * b * e) - 72 * a * c * e,
        q1 = c * c - 3 * b * d + 12 * a * e;
        s2 = s1 + Math.sqrt(-4 * q1 * q1 * q1 + s1 * s1),
        q2 = Math.cbrt(s2 / 2),
        s3 = q1 / (3 * a * q2) + q2 / (3 * a),
        s4 = Math.sqrt((b * b) / (4 * a * a) - (2 * c) / (3 * a) + s3),
        s5 = (b * b) / (2 * a * a) - (4 * c) / (3 * a) - s3,
        s6 = (-(b * b * b) / (a * a * a) + (4 * b * c) / (a * a) - (8 * d) / a) / (4 * s4);

    double[] roots = new double[4];
    for (int i = 0; i < 3; i++)
        roots[i] = -b / (4 * a) + (i > 1 ? -1 : 1) * (s4 / 2) - (i % 2 == 0 ? -1 : 1) * (Math.sqrt(s5 + (i > 1 ? -1 : 1) * s6) / 2);

    return roots;
}
Community
  • 1
  • 1
  • Have you tried BigInteger? http://stackoverflow.com/questions/849813/large-numbers-in-java – Tschallacka Jun 21 '16 at 12:49
  • @MichaelDibbets When I said "no Complex library etc" I also meant no BigInteger, since, it's not fast as I want it to be, it involves alot of long methods like add/subtract, and it's not really necessary. –  Jun 21 '16 at 12:53
  • @Bathsheba I edited the q –  Jun 21 '16 at 13:08
  • @SJuan76 Because there are more places in the formula where a radical is involved than the main discriminants in the end of it. –  Jun 21 '16 at 13:10
  • Try adding your current code. See http://stackoverflow.com/help/mcve With your unwillingless to use efficient code from established libraries, where smart people really took out the kinks you will undoubtly fall against, we need to see your code to see if we can help. – Tschallacka Jun 21 '16 at 13:30
  • @MichaelDibbets I added some code to clearify it –  Jun 21 '16 at 13:41
  • I'm no math expert, but the first place where it goes afoul where I experimented is at `Math.sqrt(-4 * q1 * q1 * q1 + s1 * s1),`. For that I found this answer of John Skeet http://stackoverflow.com/questions/7461803/why-quadratic-equations-root-result-is-nan-java – Tschallacka Jun 21 '16 at 14:02

2 Answers2

1

You can use Descartes' rule of signs and Sturm's theorem to give upper bounds on the number of roots. Conceivable this could tell you how many roots to expect. Descartes' rule is fast as it only involves comparing sign changes.

However you can add some branching in your current code to avoid the branches where you will get a complex result. 

public static double[] solveRealQuarticRoots(double a, double b, double c, double d, double e) {
    double s1 = 2 * c * c * c - 9 * b * c * d + 27 * (a * d * d + b * b * e) - 72 * a * c * e,
        q1 = c * c - 3 * b * d + 12 * a * e;
    double discrim1 = -4 * q1 * q1 * q1 + s1 * s1;
    if(discrim1 >0) {
        double s2 = s1 + Math.sqrt(discrim1);
        q2 = Math.cbrt(s2 / 2),
        s3 = q1 / (3 * a * q2) + q2 / (3 * a),
        discrim2 = (b * b) / (4 * a * a) - (2 * c) / (3 * a) + s3;
        if(discrim2>0) {
            double s4 = Math.sqrt(discrim2);
            double s5 = (b * b) / (2 * a * a) - (4 * c) / (3 * a) - s3;
            double s6 = (-(b * b * b) / (a * a * a) + (4 * b * c) / (a * a) - (8 * d) / a) / (4 * s4);
            double discrim3 = (s5 - s6), 
                   discrim4 = (s5 + s6);
            // actual root values, may not be set
            double r1, r2, r3, r4; 

            if(discrim3 > 0) {
                 double sqrt1 = Math.sqrt(s5-s6);
                 r1 = -b / (4 * a) - s4/2 + sqrt1 / 2;
                 r2 = -b / (4 * a) - s4/2 - sqrt1 / 2;
            } else if(discrib3 == 0) {
                 // repeated root case
                 r1 = -b / (4 * a) - s4/2;
            }
            if(discrim4 > 0) {
                 double sqrt2 = Math.sqrt(s5+s6);
                 r3 = -b / (4 * a) + s4/2 + sqrt2 / 2;
                 r4 = -b / (4 * a) + s4/2 - sqrt2 / 2;
            } else if(discrim4 ==0) {
                 r3 = -b / (4 * a) + s4/2;
            }
            if(discrim3 > 0 && discrim4 > 0) 
                 return {r1,r2,r3,r4};
            else if( discrim3 > 0 && discrim4 == 0 )
                 return {r1,r2,r3};
            else if( discrim3 > 0 && discrim4 < 0 )
                 return {r1,r2};
            else if( discrim3 == 0 && discrim4 > 0 )
                 return {r1,r3,r4};
            else if( discrim3 == 0 && discrim4 == 0 )
                 return {r1,r3};
            else if( discrim3 == 0 && discrim4 < 0 )
                 return {r1};
            else if( discrim3 < 0 && discrim4 > 0 )
                 return {r3,r4};
            else if( discrim3 < 0 && discrim4 == 0 )
                 return {r3};
            else if( discrim3 < 0 && discrim4 < 0 )
                 return new double[0];
       } 
   }
   return new double[0];

}

Looking further down the answers of the math exchange post there is a reference to "Finding Real Roots of Quartics" by Don Herbison-Evans, Michel Daoud Yacoub & Gustavo Fraidenraich. Which you can find here. In that paper he takes account of numerical issues.

Don't discount numerical methods. Newtons is pretty fast and can converge quickly.

Salix alba
  • 7,536
  • 2
  • 32
  • 38
0

I am not aware of a way to calculate only the real roots of a polynomial, but you could give the Apache Commons Math library a try. The LaguerreSolver can calculate all complex roots for a given polynomial.

T. Neidhart
  • 6,060
  • 2
  • 15
  • 38