Conceptual and perhaps fatal bugs remain in the stuff below.

This Algol 68 code is showing enough signs of life to deserve some comments and explanation. I use Algol 68 here for two important reasons:

The compiler that I have supports arbitrary floating precision which this algorithm notriously needs.
The syntax and semantics of the subscripts, such as omitting them, makes the program expressions near enough to how I think so as to have made debugging possible.

Much of the program consists of routines copied from other contexts and munged here mainly to do complex arithmetic. The code aspires to find eigenvalues and vectors of arbitrary square complex matrices. It has jut now succeeded for a real diagonal matrix and so I pause to document and celebrate as the method makes that case non trivial.

As one learns about eigenstuff this is one of the constructive ways to show the existence of these values and as such recommends itself as a computer algorithm.

The determinant |M−λI| is an nth degree polynomial in λ when M is a complex n by n square matrix. Finding the coefficients of that polynomial could be done by formal evaluation but I use another trick here which may well be less efficient. Lagrange interpolation produces the polynomial coefficients from knowing the value of the polynomial for n+1 distinct arguments. I choose these complex values for λ from a hexagonal grid in the complex plane. We view the n+1 polynomial values as a vector in n+1 dimensions. I evaluate |M−λI| for each of these λ’s and multiply by the matrix L that transforms that vector (“v” in routine “cp”) to the vector of polynomial coefficients (“pc” in routine “cp”).

Below we speak of the vector space of polynomials over X and two distinct bases therein.

powers of X: The jth basis vector is X^j for 0≤j≤n.
points: The jth basis vector is the polynomial that is 1 for point x_j and 0 at the other points. The points are chosen in the complex plane from a triangular grid where minimum distance between points is 1.

Routine functionality in order:

determ: Calculate determinant as it consumes its square complex matrix argument.
det: Calculate determinant without destroying its argument.
inv: Inverts a complex square matrix.
mm: Matrix multiply, on sabbatical.
ev: Computes value of polynomial given coefficients and value of formal parameter.
der: Returns formal derivative of polynomial.
sm: A crude real approximation of magnitude of complex number.
rr: A generator of random reals between −½ and ½.
root: finds one root of a complex polynomial given its coefficients.
div: Synthetic division of two polynomials, returns quotient and remainder.
allr: returns all the roots.
lagr(n): Routine returns an n by n matrix transformation from ‘points’ basis to the ‘powers’ basis. Notice that this matrix does not depend on matrix whose eigenstuff is sought. Nor does its inverse.
cp: computes the characteristic polynomial.
aeval: computes all the eigenvalues.
tot: computes all eigenvalues and eigenvectors.

Notes on a misconception:

Mx = λx
0  1
2  0
evals = ±√2
call first evc = (p, q)
q = (√2)p
2p = (√2)q
wlog p=1. ∴ q = √2 ∴ p=1
We have verified our 1st evc: (1, √2).
----
q = -(√2)p
2p = -(√2)q
wlog p=1. ∴ q = -√2 ∴ p=1
We have verified our 2nd evc (1, -√2).


Mx = λx
0  1 | 1      √2
2  0 | √2     2

0  1 | 1      -√2
2  0 | -√2     2
Now I believe it!

M=
0  1
2  0
evals=
-√2 √2
mm= (for evc=-√2)
√2 1
 2 √2
nv*10^-9=
-.5 .35355
.70710 -.5