The Algol 68 routine ‘root’ finds a root of a polynomial over the complex numbers.
A polynomial is represented as value of type P which is an array of its complex coëfficients.
Here are the meanings of the identifiers in the program:
- small
- a floating point estimate of the long long precision.
- ep(p, z)
- p(z): The value of polynomial p at z.
- der(p)
- The derivative of polynomial p.
The derivative is a polynomial.
- sm(x)
- max(|real x|, |imag x|).
- root(p)
- Some root of the polynomial p.
The root routine proper initializes a complex variable z to 0 and then proceeds to nudge z to travel a downhill course towards a root.
We minimize |p(z)| by traveling in the direction −dp(z)/dz.
This was inspired by a short proof of Fundamental Theorem of Algebra in 2011 Summer version of Mathematical Intelligencer, Vol 33, #2 by Oswaldo Oliveira.
That proof notes that there are no local minima to trap us.
This is true also of any holomorphic function.
With this insight it suffices to move in that direction, but how far?
Any time that we arrive at a z value where |p(z)| is smaller is a win, but sometimes |p(z)| is larger whereupon we cut the step in half.
The logic of the proof tells us that for some small value |p(z)| will decrease.
Two things can go wrong:
- dp(z)/dz = 0 and thus gives no direction.
- In this case |p(z)| decreases in about ½ of the nearby points.
This logic is not yet in the code.
- The symmetry is too great
- and we travel back and forth along a ridge always in a down hill direction, but not noticing that the ground is even lower on either side of us.
This motivates the use of next random which steers us a little off the otherwise optimal direction given by −dp(z)/dz.
(Buridan’s Ass)
There are some interesting routines in the tests that follow the main routine:
- pfr(p)
- returns a polynomial given the set of its roots.
Duplicate roots are respected.
- rnd(s) (too)
- returns a random normal deviate.
Its argument caries the generator state which is of type REF RNFS and which should be initialized as (1<0, ~).
- div(n, d) (too)
- does synthetic division of polynomials.
It returns quotient and remainder.
- allr(p) (too)
- Finds all roots of a polynomial repetitions replicated.
Just to the right of the two occurrences of “##” is the construct “(n|” which chooses which tests to invoke.
This notation is the ‘brief form’ of the Algol 68 case statement.
The first occurrence can vary from 1 to 3; the 2nd from 1 to 6.
Synthetic Division code;
The RC4 Algorithm;
Random Normal Deviate Generator
trouble changing precision;
Note the goto’s.
Eigenvalues and steps towards eigenvectors.
Projection as tool.