This is some logic for this code whose use is described here. It also largely applies to this older code.

Consider an n dimensional abstract vector space V. Consider two sets, A and B, of n vectors each. Suppose that each set spans V. Each member of A can be represented as a linear combination of the vectors in B and vice-versa. The coefficients in these two representations each constitute an n by n matrix and these matrices are inverses of each other.

In these programs a matrix is represented as a list of lists of scalars—a list of rows, each row is a list of scalars. (A scalar is an element of the field.) For an n by n matrix the inversion routine may be understood as taking n column vectors and appending n more columns, each describing a basis vector. For n = 3 the picture is:

a11  a12  a13  1  0  0
a21  a22  a23  0  1  0
a31  a32  a33  0  0  1
The 6 columns correspond to the 6 vectors of A ∪ B. The 3 scalars of a column are the coefficients of the linear combination that forms the corresponding vector in terms of the vectors of set B.

Inside the program this is a list of 3 lists of lists each of 6 scalars.

The named let called pad appends these extra columns. It omits trailing zeros which are implicitly supplied as needed.

We apply a sequence of elementary operations to modify the representation of each of these 2n vectors. We justify the operations of the program in light of a sequence of changes to the basis in which our constant 2n vectors are expressed. These changes are thus to be considered as alias transformations—a sequence of new designations for a constant set of 2n vectors in an abstract space V. Each such operation expresses the 2n old vectors in terms of a new gradually changing vector basis. This changing basis is initially the set B and after all our elementary operations the changing basis is the set A. At the end stage the left n2 scalars will form the identity matrix and the right n2 scalars will be the matrix inverse of the initial left n2 scalars. The new picture is thus:

1  0  0  b11  b12  b13
0  1  0  b21  b22  b23
0  0  1  b31  b32  b33
The b’s are the inverse matrix, for the right n columns have continued to depict the old basis vectors thru a sequence of evolving coordinate systems, while the right n rows described and still describe the original input vectors in that same sequence of coordinate systems. The right n rows thus describe the old basis in terms of the input vectors, which is literally the inverse of describing the input vectors in terms of the old basis.

Alas the code is shorter, more efficient and more obscure than these figures suggest. Not only are zeros at the ends of rows omitted but predestined zeros and ones at the beginning of rows as well. I don’t have a clear explanation of why this is more than a space saving hack but the code is both shorter and faster this way. I hesitate to make it longer, slower and consume more space, in order to make it correspond to the simple process above. The definition of matrix inverse has not changed in 170 years and I thus think the code need not evolve. (Note 2013: I now study this to extract the magic of finding x such that Mx=0 to supply that magic in a language that does not abstract as well as Scheme.)

The elementary operations are:

Named let inx takes the padded matrix and recursively: Upon completion inx returns an upper triangular matrix with the leading zeros and ones of each row suppressed. This value is named ut. Of course this value is still padded with the transformed names of the old basis vectors.

Now named let ol, to which this was all prelude, takes ut and eliminates the values above the diagonal by yet more elementary transformations.

This may be relevant.

For the logic of a report on singular matrices, see this.
We speculate now on a generalization to the Penrose inverse. To compute the inverse of A we form the doubled matrix [A, I] and derive a sequence of matrices [X, Y] preserving the invariant X = AY where [A, I] is the first in the sequence and [I', A+] is the last. I' may have some zeros on the diagonal. Perhaps X = YA is a more convenient invariant. I must understand this before I proceed but that scheme which depends on some notion of orthogonality may fail for finite fields.