A little, open-source polynomial solver.

Purpose Of Exercise:

Historically, techniques have existed to approximate the roots of a polynomial by purely numerical means. These approaches are eventually needed, because polynomials of degree greater than 4 may arise, but have no general solution. For example, the only situation in which a polynomial of 6th degree has an exact, analytical solution could be, if it also happens to be the product of two cubics. Yet, sextics can arise in Engineering etc., which at least require a numerical solution.

In Math, the concept of approximation is relative, in that approximations will be accurate to within a certain number of decimal places or bits. Hence, "3.14" is an approximation of Pi. Yet:

3.1415926535897932384626433832795028841971693993751058209

Is also an approximation of Pi, of which the first 50 decimal places are accurate.

Just to prove it could be done, I wrote an approximation algorithm which is not in the History Books, but which is entirely home-spun. It's based on a simple premise. In general, two Bases (Say, A and B) can be interpolated between, with the equation:

C = A*(1 - t) + B*t

Where (t) is a floating point number, which needs to be between (0.0) and (1.0) in order for (C) actually to be an interpolation between (A) and (B). This can be paraphrased, that (t) states, by how much (B) should replace (A) in determining (C).

To make this more interesting, the premise could also be, that (A), (B) and (C) are already known, but that (t) needs to be derived to fit the equation. And so, the equation to do so is:

t = (C - A) / (B - A)

In this case, if (A), (B) and (C) are just any values, we also have no guarantee that (t) will be a floating-point number between (0.0) and (1.0). So, how can this be used?

General Principle:

My Algorithm assumes that the following approach will generate values of (X), such that F(X) will become successively closer to (0.0):

Two initial values of (X) will be applied, and F(X) will be computed for each of them. Then, the value of (t) will be found as described above, such that -

0 == F(X1)*(1-t) + F(X2)*t

This is an exact process, which will initially fail to achieve that -

0 ~= F((X1*(1-t)) + (X2*t))

(Because F(X) is Not a linear function.)

Instead, what will follow is that -

Error := F((X1*(1-t)) + (X2*t))

In the next iteration, we will apply:

X1' := X1*(1-t) + X2*t

X2' := Midpoint(X1', X1)

t' := (0 - F(X1')) / (F(X2') - F(X1'))

(...)

Complex Roots:

I deem that this approach will also find complex roots, with one catch. If the system has not received any complex numbers, complex interpolations will also not be derived. And so, an alteration will be made, to the idea -

X2' := Midpoint(X1', X1)

This will be replaced with -

X2' := X1' + ((X1 - X1') * i * 0.5)

In other words, the geometry with which (X2') trails behind Δ(X1), will be rotated 90⁰ in the complex plane. This will guarantee that complex solutions are offered, the imaginary part of which would simply become progressively smaller, as the Error does, if the true solution had an imaginary component of (0). In that case, ΔError will also become progressively smaller.

Multiplicities:

One scenario which acts as an enemy of my solver is any situation where the derivative of F(X) becomes close to zero, close to a root. The reason for that is the fact that near such values of (X), considerable changes in (X) will result in almost no change in F(X), instead forming pairs of F(X) that are both close to zero. Here, F(X) seems to form a tangent with the X-axis. Hence, there will be no ΔF(X) to base further interpolation on. As it happens, the derivative of a polynomial of degree (n) is a polynomial of degree (n-1). Thus, if we can find roots of the polynomial of the (n-1) degree, we've also found values for (X) at which F(X) has a derivative of zero.

Problem: Just as F(X) has a derivative of zero for some value of (X), so does F'(X). Thus, to find the roots of F(X) reliably, and the degree of F(X) is (n), we'd need to find the roots of F'(X) reliably, and recursively all the roots down to those of F^'(n-1)(X). But, since I did not want to add recursion to my algorithm, I found a more robust solution. I could first define an operation which 'normalizes' the polynomial(s), which means, to divide each polynomial by the coefficient of its highest exponent, so that that coefficient becomes equal to (1.0), but a normalized polynomial will have the same roots as the original polynomial. And then, down to the work, a subroutine can find the derivative of each polynomial and normalize it successively, right down to a resulting linear equation.

The sub-task of normalizing each polynomial was essential, since to find all the derivatives of a polynomial of the n^th degree, would otherwise start to generate coefficients of the highest exponent, that would reach (n!), the factorial of (n), which could become huge. What will happen with my algorithm instead is, that while the highest coefficient stays equal to (1.0), all the lower coefficients become correspondingly weak.

For every value of (X) where the derived polynomial is close to zero, a check takes place - that has some intricate logic - of whether its parent polynomial is also close to zero. And if so, that root is factored out of the parent polynomial the required number of times, since the situation represents a multiplicity of the root.

In some cases, roots in the derivatives actually correspond to roots in F(X), in which case the original polynomial has a root with multiplicity (which has been factored out as many times as required, before the solver is ever run on any parent polynomial).

Caveat when Compiling under MSYS2, MINGW, Etc.:

It is a fact which I only discovered long after I had published this source code, that the header files between my Debian 9 computer, and the MSYS2 Linux-compiler designed for Windows, have a subtle difference, that can cause compiling to come to a halt. Under MSYS2, the error message can get printed, that 'numeric_limits' is not a member of 'std'. If this did happen to the reader, then the only way he or she could have continued with compiling the code would have been, to use a text editor, and add the following line to the includes:

        #include <limits>

Many readers might not even have bothered, or known which source or header file to add this to. As of May 14, 2022, this problem is fixed in the version of the source code found on my site.

On many systems, including "Visual Studio 2019", the other header files I included, also pulled in that one. But on some compilers, they do not.

Refinements in the source code:

Over the years, I've made modifications that improve the performance, that include a polishing algorithm...

Additionally, I've created a GUI-based variation of the program, which is Not more accurate, but which might be easier for people to use, who do not like the Command-Line Interface of the first version. The following are two screen-shots of the GUI-based version:

Theoretically, one might suspect that I'd have no idea, what the solutions are supposed to be, for a polynomial of the 12th degree. But in this case, I've chosen that the coefficient of (x¹²) will be (1.0), while the constant component will be (-1.0), and all the other coefficients are (0.0). This means that the roots will also be the 12, 12th roots of (+1.0), 2 of which are Real, 2 of which are purely Imaginary, and another 8 of which are Complex.

What's left to be determined is, whether I'd know the decimal expansion of cos(30⁰) off the top of my head. No. I know that this is the square root of (3) divided by (2), but find it difficult to compute in my head, in decimal notation. But I can see that this constant repeats itself symmetrically, just as sin(30⁰) repeats itself. BTW, sin(30⁰) == ½.

Location of the code:

https://dirkmittler.homeip.net/binaries/

I don't want to give people links which will just download any software directly when clicked. The link above will lead to a (lengthy) table of contents. In that table, the files which might interest the reader would be:

Roots_Dirk_Mittler.tar.gz

Roots_Dirk_Mittler.zip

RRoots_Dirk_Mittler.tar.gz

RRoots_Dirk_Mittler.zip

LRoots_Dirk_Mittler.tar.gz

LRoots_Dirk_Mittler.zip

Dirk_Roots_GUI_1.tar.gz

Dirk_Roots_GUI_1.zip

Update April 7, 2024:

I've created a smartphone-optimized version of the GUI project, the source code for which can be found here...

Dirk_Roots_GUI_2.tar.gz

Dirk_Roots_GUI_2.zip

There's also a compiled Android package for arm64 -CPU, mobile devices...

Dirk_Roots_GUI_2-arm64-v8a.apk

I do not really recommend the armeabi-v7a version, because it's missing certain updates, and I've lost my ability to compile to this 32-bit platform.

- Enjoy,

Dirk Mittler

The polishing strategy:

After my main solver finds a sufficiently close approximation of one root, of course the thing to do is, to factor it out of the original polynomial, in a process known as "Deflation". What one may find after one has used polynomial division [by (x - r)] is, that some non-zero remainder is left, in spite of a first attempt to make (r) close to a root. Well, at some point I realized that the magnitude of this remainder, with respect to the term being factored out, is proportional to the already-approximated quotient, of the deflation. Thus, my polishing strategy is, to divide the remainder by the value of the quotient, and then to add the result to (r), to arrive at a new value of (r), and a new quotient.

Unfortunately, this strategy shares a weakness with that of the main solver, which is, that polishing will not work, if the quotient is equal to zero. This is the case revisited, where the original polynomial has a multiplicity, as the quotient shares a root with it.

And so, as much as my main solver does, my polishing strategy relies on multiplicities already having been factored out of the polynomials.

Complex coefficients:

Because one of the things which my program does is to find complex roots, and because those tend to occur in conjugate pairs, one of the initial questions to my mind were, whether the process of deflation should try to factor an entire quadratic out of the polynomial in one step. And what I found was, that the code would just be simpler, to allow each coefficient to be a complex number internally, to factor out one complex root, resulting in a quotient with complex coefficients, and then to factor out the conjugate of the previous root, so that the next coefficients will go back to being real numbers, if that is what they were in the original polynomial.

One might also say, that the complex roots will cease to occur in conjugate pairs, if the coefficients were complex.

In order to try to minimize the errors as much as possible, my algorithm will have polishing turned on for the first root, but, if the value of the deflated polynomial is close to zero at the conjugate, will turn polishing off when the conjugate is next factored out.

The original version of my program, which can only be used from the command-line and with some difficulty, exposes this ability to specify complex coefficients to the user. However, both the simplified CLI version and the GUI version no longer do so, so that their ability to work with complex coefficients is only internal.

The versions of the program which still accept complex coefficients can be found in the files 'RRoots_Dirk_Mittler.tar.gz', 'RRoots_Dirk_Mittler.zip', 'LRoots_Dirk_Mittler.tar.gz' and 'LRoots_Dirk_Mittler.tar.gz',. Here is an example of what that looks like, in my computer's Terminal Window, named "Konsole":

Large Coefficients:

One of the tasks which some people might want to fulfill with such a utility is, to be able to specify coefficients which have a 'large' absolute, where large could mean (-1e+30). However, the original program needed to be modified in order to be able to accommodate that. As of August 18, 2022 I have published the version of the program that does. It's a Command-Line -driven version, that works as 'RRoots_Dirk_Mittler.tar.gz' did, but is named 'LRoots_Dirk_Mittler.tar.gz' and 'LRoots_Dirk_Mittler.zip'.

A trade-off this version offers is, that its internal constant (epsilon) has been set more coarse than it was in ''RRoots_Dirk_Mittler.tar.gz'. Essentially, what this parameter determines is, how far apart two roots need to be, so that the initial search algorithm can recognize them as distinct. But then, polishing, if it has been left enabled, will still maximize the precision of the found roots.

Memory Leaks:

On August 29, 2022, I ran 'valgrind' on my code, and found that it had memory leaks. 'valgrind' is a program of the type which is known as a "profiler". In fact, it's the prescribed profiler to be using under Linux. Profilers are programs, which give programmers metrics on running programs. They are useful for debugging code. One question which the reader may ask himself or herself is, 'Is it crucial, that no code ever have memory leaks?' And the short answer is, 'It depends.' If the memory leak happens once, directly before the code exits, the O/S kernel tends to clean up whatever memory the programmer's source-code failed to deallocate, so that basically, the answer is No. However, if the programmer's allocation happens frequently as the program is running, then the result will be the kind of program, which just continues to consume more and more RAM, until the user closes and restarts it. So here, there should really be no memory leaks. If it's a program that just runs once from the command-line, then the answer is more or less, 'It matters, because what the program does could consume a lot of memory, before it quits.'

The first reason for which the code I published before today leaked memory was the fact, that the function 'destroy_deriv_list()' uses the '.deg' property of the first polynomial, in the derivatives list, to determine what the original degree of the polynomial was. The problem next becomes, the fact that the function 'divide()' reduces this property's numeric value, every time it factors out a root. Hence, what the degree of the original polynomial was, was forgotten as soon as it was solved, and when the time had come to destroy the derivatives list

I solved this problem in my source code, by giving my 'Polynomial' struct an additional property, that just saves what its original degree was. I then modified 'create_deriv_list()' and 'destroy_deriv_list()', first to set this property, and then to use it, in order to determine how many derivatives to delete. These updates to the source code were not compiled into the numerous binaries that should result each time. If the reader wants the benefit, then he or she need only recompile the provided source code...

When I next modified the GUI version of the program, what 'valgrind' told me was, that in addition to some start-time and stop-time oddities, there was really only one more memory leak for me to fix. This had to do with changing Line 75 of the source file 'plot.cpp', from:

To:

And so, what I have achieved in the published source code is, that the GUI application can be operated by the user as often as desired, resulting in no cumulative memory leaks. I'll just leave it to the reader again, to recompile that, in order to reap the benefits... What should be familiar is the fact that, even if the programmer has a C++ compiler, he or she will always need the Qt SDK, to be able to compile any Qt, GUI application. Also, on my Linux computer, the GUI version of my program produces some errors according to 'valgrind', that stem from how the GUI library communicates with my X-server. These little bugs are not due to my code at all, and may not be reproducible by the reader.

~

Mastodon