Gaussian_elim.at

When doing a pivot step, we either find a pivot, or not. In the former case
we record the pivot row number, as well as the row operation to be applied
after the cyclic row permutation bringing the pivot to the echelon position is
done (row number equal to the number of previous pivots). The latter row
operation is left-multiplication by a matrix differing by only one column from
the identity, and that column is recorded. Since we are doing (full)
Gauss-Jordan reduction and the row permutation is implied by the pivot number,
this column has the inverse of the pivot value on its main diagonal position,
and all other coefficients can be arbitrary. In case no pivot is found (the
part of the column where a pivot might be found is zero) we instead record the
nonzero initial part of the column as the |no_pivot| variant of |Pivot_info|.

Once accumulated, the different |Pivot_info| values are combined, taking the
|pivot| fields together (extracting them for the union type), while
transforming the information of non-pivot cases into a matrix whose columns
are independent and span the kernel. A third component of the result is a
function that transforms a row-reduced RHS vector, with as many components as
there are pivots, into a full size vector solution, by mapping each entry to
the corresponfing pivot position, and interspersing zeros at the non-pivot
positions. This function |spread| is constructed by repeated modification,
each newer version first calling the older one and then taking one more entry
to its proper destination. We avoid inadvertently building a recursive call,
by first binding the old function value to a new name |f|, and calling a
function of that name from the newly extended function.

Now put the above together for solving a rational linear system |A*x=b|