Simple, numerical computing.

Sometimes, I get these ideas, and I just have to sit down, write it down and get it out of my head. If it's a good idea, I prefer to do it in such a way as to have some use of it, and sometimes I do it just to see if it will work. I'm talking about ideas for computer programs, of course, and when I say write it down, I mean writing it down in a programming language and then possibly run the program. I've had one of these days today.

Yesterday, a flatmate of mine was sitting with a ninth degree polynomial he needed a solution for. A ninth degree polynomial is a hairy thing, and you don't really want to try solving it by hand, for one thing, Niels Henrik Abel proved (if I recall correctly), that there can be no explicit formula or correct algebraic formula for finding a solution of any polynomial of degree 5 or more. That is, no nice and simple x = (-b ± sqrt(b^2 - 4 ac)) / 2a, like we do for degree 2 polynomials, or even something akin to the more complex formula for degree 3 polynomials.

On the other hand, in most circumstances, we can find a solution for a polynomial, provided there is one, by numerical methods, like Newtons method, for instance. However, Newtons method is impractical to do by hand, considering that it is iterative, requires you to derive the entire polynomial and calculate the value of it and its derivative for every iteration. It might not even converge, in which case... You've wasted a lot of time. In other words, doing it by hand is repeatetive, error-prone, boring and possibly useless. Whenever I have something that fits that description, I know that my computer can do it more accurately, faster and better than I can. So I thought about it when I went to bed last night, and I had a head filled with ideas when I woke up in the morning.

The basic problem of how to implement this, is how to represent a polynomial in an effective and tidy way. A polynomial is of the form an * x^n + ... a0 * x^0 = 0. This seems to hint towards implementing it as a list of lists, or at the very least, a container of containers. However, this might not be what one wants, and it is a more complex datastructure than what you need for this problem. Consider the row vector [an ... a0] and the coloumn vector (x^n, ... , x^0). The dot product of these two vectors is an * x^n ... a0 * x0. If we reverse them both, we see that, in an environment where containers are zero-indexed, we get a row vector [a0 ... an] where the index of a given coefficient corresponds directly to the power we want to raise x, so that the value of the polynomial in a given point can be calculated with a simple for loop like so:

function evaluate at point(polynomial, point):
    result = 0
    for index, coefficient in polynomial:
        add coefficient * point^index to result
    return result

This is remarkably simple and logical, but there is another added benefit. Newtons method does, as earlier mentioned require derivation. Deriving a polynomial is a simple but tedious process, f(x) = x^n, f'(x) = n * x ^ (n - 1). With our representation, we can do this easily with a loop as well, we multiply every coefficient by its index, then chop off the first element. To see why this is so, consider: C = [a0 a1 ... an]; after multiplying every element by its index, C' = [0 a1 2a2 3a3 ... nan]. Now we chop off the first element, which means C' = [a1 2a2 3a3 ... nan] where a1 is the same coefficient as a1 in C, its index is 0, so every element in the list is now correctly derived, from which it follows that the list we now have corresponds to the derivative of the list we started out with. Newtons method should now be trivial to implement, however there are a few more considerations to be done here. Firstly, Newtons method does not always converge, or it might converge very slowly. And by the definition of the method, we have to choose a guess such that for our polynomial P(x), P'(guess) differs from 0, because P(guess) / P'(guess) will obviously not be defined otherwise.

We can however, make better guesses by thinking a little. Assuming that we don't have any a * x ^ n where n < 0, our polynomial and its derivative should both be continous. That is, if we have y and z such that P(y) < 0 < P(z), x such that P(x) = 0 is in between somewhere and y < z or z < y is true (from the Intermediate value theroem). And if we guess closer to the actual answer, we're better off because we will require fewer iterations, in addition to not getting screwed over not converging quite so often. However, since we can still be unlucky and not converge, we will put a roof on how many iterations we can actually do. Next, we need some sort of tolerance, that is, how precise do we want our answer to be? When we have clarified these issues, we can quite easily implement Newtons method such that it can work on any polynomial.

function derive(Polynomial):
    for index, coefficient in Polynomial:
        Polynomial[index] = coefficient * index
    remove first from Polynomial
    return Polynomial

function newtons method(Poly)
    guess = make initial guess
    Poly' = derive(Poly) 
    while iterations < roof and abs(evaluate in point(Poly, guess)) > tolerance:
        guess = guess - (evaluate in point(Poly, guess) / evaluate in point(Poly', guess))
    if abs(evaluate in point(Poly, guess)) <= tolerance:
        return guess

I wrote this program in C, representing a polynomial with pointers, and using (more or less) the exact methods I've described here, and I've found it to work very well. Unfortunately, it turns out my flatmate is done with the huge polynomials for now, so the program will gather dust, until I decide to clean up my hard-disk, read the code again and wonder 'what the hell is this guy doing?'

I did, however, have a whole lot of fun when coming up with it, and later when writing it, and I do feel enlightened.