Excusatio Non Petita

The year probably was 1994, so I guess its understandable that I don't exactly recall how did I come across a copy of Numerical Recipes in C: The Art of Scientific Computing. But while skimming through the pages, I eventually got to the part on evaluating polynomials, and read something that struck me very deeply. The paragraph that described what I now know is Horner's scheme read as follows:

We assume that you know enough never to evaluate a polynomial this way:

p=c[0]+c[1]*x+c[2]*x*x+c[3]*x*x*x +c[4]*x*x*x*x;

or (even worse!),

p=c[0]+c[1]*x+c[2]*pow(x,2.0)+c[3]*pow(x,3.0)+c[4]*pow(x,4.0);

Come the (computer) revolution, all persons found guilty of such criminal behavior will be summarily executed, and their programs won’t be! It is a matter of taste, however, whether to write

p=c[0]+x*(c[1]+x*(c[2]+x*(c[3]+x*c[4])));

or

p=(((c[4]*x+c[3])*x+c[2])*x+c[1])*x+c[0];

Boy was I guilty! Who would have ever thought there was another way of evaluating a polynomial than the obvious one? And that it would be so much more elegant and efficient?

Would be nice if it was just a polynomial thing, but of course it isn't, and the more complicated the math you are trying to get your computer to do, the larger the chance you are doing it terribly bad. And although much has been written about algorithms, the gospel doesn't seem to be spreading too swiftly: as I am writing this, the first match for a search of "prime numbers python " that google returns is this. And he is not alone...

So what am I in for? The idea is to bridge the gap between straightforward brute force (a sin I have been as guilty as anyone else of) and the elegant simplicity of efficient code. Hopefully, this will allow me to escape the guillotine when the Jacobins take over...

There are some simple rules that I intend to follow, at least until I decide not to:

1. Uncomplicated complexity. The math content will be kept reasonably simple. Or maybe simple isn't the word, as things may get utterly complex. But I will try to stay out of unnecessary complications, so a high school diploma should be all you need to come along on the ride...


2. Plain python. Sure, C++ is a zillion times faster, but python takes care of most gory details, and makes it easier to concentrate on the algorithm. Plus, it's so cool and trendy!


3. No libraries. Sage, SciPy, PARI, qhull, FFTW... whatever it is you are trying to achieve computationally, there is a faster way than replicating the code you find in this blog. And interesting as that may be, it would take a lot of the charm of writing this blog, if it turned into a description of somebody else's API. So yes, we will even code our own FFTs!