I think I have a Quintic equation solver

Started by
4 comments, last by Narf the Mouse 14 years, 1 month ago
The basic idea is to recursively converge on a solution, correcting for errors in convergence. Although, of course, I cannot guarantee that it will find all the real roots, nor, due to the limitations of the math available to me, will it find complex roots. However, all the real roots it finds will be accurate to a high degree of accuracy. Due to rounding errors,it may find duplicate roots that are only an extremely small degree apart (Say, +/-0.0000000009). As well, it correctly finds roots for a Linear, Cubic, Quadratic and Quartic equation. The roots it finds are for 0 = Polynomial Formula. It may also be able to solve other types of equations, but I haven't tested it. You may use it freely, provided you give credit where credit is due. :) Warning: Not very commented

        public static List<Double> EquationSolver(Func<Double[], Double, Double> f, Double[] inputs)
        {
            List<Double> list = new List<double>();
            EquationSolveWithin(f, inputs, 0, Double.MaxValue, ref list);
            EquationSolveWithin(f, inputs, 0, (-Double.MaxValue), ref list);
            return list;
        }


        public static void EquationSolveWithin(Func<Double[], Double, Double> f, Double[] inputs, Double low, Double high, ref List<Double> list)
        {
            Double current = (high / 2.0) + low;
            Double output = 0;
            Double? firstOutput = null;
            bool reverse = false;
            Double distance = high * high * high;
            Double lastDistance = high * high * high;
            Double oldLow = low;
            Double oldHigh = high;
            Double tempLow = low;
            Double tempHigh = high;
            Double? closest = null;
            Double closestDistance = 0;
            bool triggerHigh = false;
            int iterationsWithoutBetter = 0;
            bool foundCloser = false;

            do
            {
                current = ((high - low) / 2.0) + low;

                tempLow = low;
                tempHigh = high;


                lastDistance = distance;
                output = f(inputs, current);
                if (!firstOutput.HasValue)
                {
                    firstOutput = output;
                    distance = output;
                    lastDistance = output;
                    closestDistance = output;
                    closest = current;
                    foundCloser = true;
                    iterationsWithoutBetter = 0;
                }
                distance = output;
                if (high >= -100)
                {
                    int j = 0; // Break-point anchor. Does nothing else.
                    ++j;
                }

                if ((distance > 0 && lastDistance < 0) ||
                    (distance < 0 && lastDistance > 0) ||
                    (Math.Abs(distance) > Math.Abs(lastDistance)))
                    reverse = !reverse; // Corrects for overshooting.

                if ((distance <= lastDistance && firstOutput > 0) ||
                    (distance >= lastDistance && firstOutput < 0))
                {
                    if (!reverse)
                    {
                        high = ((high - current) / 2.0) + current; // Helps the equation not overshoot
                        if (high == tempHigh)
                            high = current;
                        triggerHigh = true; // Might be useful, at some point.
                    }
                    else
                    {
                        // low = current;
                        low = ((current - low) / 2) + low; // Helps the equation not overshoot
                        if (low == tempLow)
                            low = current;
                        triggerHigh = false;
                    }
                }
                else
                {
                    if (!reverse)
                    {
                        // high = current;
                        high = ((high - current) / 2) + current; // Helps the equation not overshoot
                        if (high == tempHigh)
                            high = current;
                        triggerHigh = true;
                    }
                    else
                    {
                        // low = current;
                        low = ((current - low) / 2) + low; // Helps the equation not overshoot
                        if (low == tempLow)
                            low = current;
                        triggerHigh = false;
                    }
                }

                foundCloser = false; // Occasionally, it'll bounce back and forth around a number range that won't work. This is an adaquate solution, for now.
                ++iterationsWithoutBetter;
                if (Math.Abs(distance) <= Math.Abs(closestDistance)) // Keep a record of the best found.
                {
                    closestDistance = output;
                    closest = current;
                    foundCloser = true;
                    iterationsWithoutBetter = 0;
                }

                // Occasionally, it would get stuck where High and Current were just under the optimal solution, or the relative same for Low and Current (Just over)
                if ((closest.HasValue && high < closest && oldHigh > 0) ||
                    (closest.HasValue && high > closest && oldHigh < 0))
                {
                    if (!reverse)
                        high = closest.Value * 2;
                    else
                        low = closest.Value / 2;
                }
                if ((closest.HasValue && low > closest && oldHigh > 0) ||
                    (closest.HasValue && low < closest && oldHigh < 0))
                {
                    if (!reverse)
                        low = closest.Value / 2;
                    else
                        high = closest.Value * 2;
                }

                
                // Failed solvers. Left in there for "archeological" reasons - And just in case.
                /*
                if (distance <= lastDistance)
                {
                    if ((output < 0 && high < 0) ||
                        (output > 0 && high > 0))
                        high = current;
                    else
                        low = current;
                }
                else
                {
                    if ((output < 0 && high < 0) ||
                        (output > 0 && high > 0))
                        high = current;
                    else
                        low = current;
                }
                */
                /*
                if (output > 0)
                {
                    /* if ((output >= outputLast && oldHigh < 0) ||
                        (output <= outputLast && oldHigh > 0)) */
                    /* if (oldHigh > 0)
                        high = current;
                    else
                        low = current;
                    if ((output >= outputLast && oldHigh < 0) ||
                        (output <= outputLast && oldHigh > 0))
                        high = current;
                    else
                        low = current;
                    // high = current;
                }
                else if (output < 0)
                {
                    /* if ((output >= outputLast && oldHigh < 0) ||
                        (output <= outputLast && oldHigh > 0)) */
                    /* if (oldHigh < 0)
                        high = current;
                    else
                        low = current;
                    low = current;
                }
                /* if (output > 0)
                {
                    if ((output > 0 && high > 0) ||
                        (output > 0 && high < 0))
                        high = current;
                    else
                        low = current;
                }
                else if (output < 0)
                {
                    if ((output > 0 && high < 0) ||
                        (output < 0 && high < 0))
                        low = current;
                    else
                        high = current;
                }
                */
            // The first one is the solution checker.
            // The others are to make sure it doesn't get hung up.
            } while (output != 0 && high != low && tempHigh != current && tempLow != current && (foundCloser || iterationsWithoutBetter < 50)); 

            if (!list.Contains(current) && output == 0) // If we haven't found this one and it's a solution,
            {
                list.Add(current); // Add it and check on both sides.
                EquationSolveWithin(f, inputs, oldLow, current, ref list);
                EquationSolveWithin(f, inputs, current, oldHigh, ref list);
            }
        }



Advertisement
I didn't look at your code, but if you want to numerically find all roots( including complex) of an arbitrary real coefficient polynomial then you probably want Bairstow's Method. If I remember correctly it's basically an iterative 2D Taylor Approximation that finds quadratic factors(which can get you complex roots via quadratic formula). Then you perform synthetic division and repeat until done. The Wikipedia article should have enough info, but I've got a java implementation if you're interested in seeing the code. The whole "impossible to solve quintic polynomials" thing only refers to a non-iterative formula for exact answers, and I'm pretty sure that's been proven impossible.
According to Wikipedia, a small sub-set are solvable, although the math is way over my head and I didn't examine it very much.

And yeah, I made this because I knew they were solvable, just not strictly mathematically. And I was pretty sure someone had at least done this before - However, it was still pretty awesome, as I have taken no college math or programming courses - Or college courses at all.

I'll look up Taylor Approximation. I would be interested in the code, as well.

Edit: Found the Wikipedia article. I get the basic idea - Each iteration gets you closer to the actual result, much like mine - Only using a mathematical function to converge, instead of iterative "guesses".
Here's the code: www.alrecenk.com/stuff/rootfinder.zip
You're in luck. I wrote it for a class so it's actually fairly readable, well aside from all the math anyways.
Thanks; I'll look at it after I get some sleep.

Current status: Now should calculate for N = Polynomial Equation.
Current status: Seems able to solve for anything with one variable to find.
Two variables would be rather beyond it.

This topic is closed to new replies.

Advertisement