// $Date: 2001/07/07 15:53:59 $ $Author: Manuel Kauers $ $Revision: 1.2 $
/**
This file contains the algorithmical parts of the parametric Gaussian Elimination package.
It provides a function elim that determines the solution set of a parametric linear
system over all the partitions of its parameter space.
The algorithm is customizable by setting particular flags. Those flags are listed below.

COLUMN_SWAP_ALLOWED
TRUE to allow column swapping during the elimination process.

USE_MARKOWITZ
TRUE to use an adapted version of the Markowitz criterion for pivot selection.

USE_COLUMN_SIMPLIFICATION
TRUE to simplify a column if no constant pivot is found. The idea is to try to
reduce the degree of the entries in a particular column in order to recieve an entry
with zero degree which then can serve as a pivot element.

SIMPLIFY_AFTER_BRANCH
TRUE to specify that affecting the constraint store will lead to a simplification
of all the polynomial entries of the remaining submatrix. Note that this option may
lead to a fillin of new variables.

USE_FURTHER_IMPROVEMENTS
TRUE to allow a number of nondocumented heuristics which slightly improve the
algorithm's efficiency.
By default, all the above flaged are set to TRUE.
It is easy to analyze the behaviour of this implementation by use of the methods defined in
Test.mu. Those methods are based on some portions of code that were inserted into
the methods of this file. In particular, the implementation here calls methods of
Test.mu to pass several statistical results to the testing architecture. Each line
of code in this file that exists for the purpose just described ends with the character
sequence "// ###". So if you want to slightly speed up the implementation,
you can feel free to delete all the lines ending with "// ###"  they are
not a part of the algorithm. (Design pattern "Template Method")
*/
/*++
elim  normalizes a matrix w.r.t. a set of constraints by means of Gaussian Elimination.
elim(A [, C])
A  matrix with entries e that are liftable to DOM_POLY via poly(e)
C  (optional) set of constraints given as a ConstraintStore object
If no C is given, C is assumed to be empty.
elim is the user interface to invoke the Gaussian Elimination on a matrix A.
The method returns a list of pairs [A', C'] where C' are specialisations of
C and A' is the Gaussian normal form of A with respect to C'.
++*/
elim
:= proc(A : Dom::Matrix(), C = null())
local L, vars, v, i, j;
save zero, width, height;
begin
// Once for all, what are the dimensions of the matrix?
// Those values are stored in global variables to avoid frequent recalculations
width := linalg::ncols(A);
height := linalg::nrows(A);
// This is a function that determines the set of all identifiers within an expression,
// e.g. vars( f(x, g(y)) ) = {x, y}.
vars := x > if nops(x) > 1 then
if type(x) = DOM_POLY then
_union(vars(op(x, 1)));
else
_union(vars(op(x, i)) $ i=1..nops(x));
end_if;
elif op(x, 1) <> x then
vars(op(x, 1));
elif type(x) = DOM_IDENT then
{x};
else
{};
end_if;
// Use the above defined function to determine all variables inside the given matrix.
v := _union(vars(A[i, j]) $ i=1..height $ j=1..width);
// add private variables for the constraint store.
v := v union {cs_var.i $ i = 1..height * width + 1};
// add the variables of the given constraint store (if it exists)
if C = null() then
// No constraints were given, initialize an empty store.
C := ConstraintStore(poly(0, [op(v)]));
v := [op(v)];
else
// The constraint store's variable list has to contain all the variables of v.
//print(v);
//print({op(varlist(C))});
if nops(v minus {op(varlist(C))}) > 0 then
error("The given constraint store does not have enough variables.");
end_if;
v := varlist(C);
end_if;
// v is a list now.
zero := poly(0, v); // this is the zero element of our polynomial ring (global definition)
// lift the matrix' elements to polynomials over R[v];
for i from 1 to height do
for j from 1 to width do
A[i, j] := poly(A[i, j], v);
if A[i, j] = FAIL then
// FIXME
error("Up to the moment, only matrices with polynomial entries are supported.\n".
"Did not succeed to convert A[".expr2text(i).", ".expr2text(j)."] to a ".
"polynomial.");
end_if;
end_for;
end_for;
// NYI: Elimination of denominators and filling them into the constraint store.
elimRec(A, C, 1, 1); // Start the recursion.
end_proc:
/*
elimRec  an auxiliary method used by elim.
elimRec(matrix, constraints, results, row, col, d)
matrix  the whole polynomial matrix that is to be normalized at the moment
constraints  ConstraintStore object representing the set of currently active constraints
row  row index of the top left element of the submatrix currently active
col  column index of the top left element of the submatrix currently active
elimRec returns a list of pairs [A', C'] where C' is a
specialization of the given constraints and A' is the normal form of the given matrix
with respect to C'.
This method is for internal use only  no checking is done.
Though the method is declared to use for recursive use, most of the work will be done
iteratively. Recursive calls are only needed if no proper pivot element is found, i.e.
the selected pivot is neither zero nor constant w.r.t. the current constraint store.
In this case, one of the two branches is done by a recursive call of elimRec, while
the calling procedure will continue to process the other branch itself.
*/
elimRec
:= proc(matrix : Dom::Matrix(), constraints : ConstraintStore, row : DOM_INT, col : DOM_INT)
local r, c, pivot, tmp, i,
old, // ###
result, b, d;
save mat;
begin
//print(constraints, matrix);
mat := matrix; // define the matrix global.
result := []; // initial list of solutions.
// Simplify the matrix's entries with respect to the current constraints
simplifySubmatrix(row, col, constraints);
// As long as possible, do the elimination job iteratively
while col <= width and row <= height do
// select a pivot element
pivot := choosePivot(constraints, row, col);
r := op(pivot, 1);
c := op(pivot, 2);
pivot := mat[r, c];
// swap rows and columns to get pivot at [row, col] reimplementing
// mat := linalg::swapRow(mat, row, r, 1..width);
// mat := linalg::swapCol(mat, col, c, 1..height);
// to avoid matrix copy
if row <> r then
for i from 1 to width do
tmp := mat[row, i]; mat[row, i] := mat[r, i]; mat[r, i] := tmp;
end_for;
end_if;
if col <> c then
for i from 1 to height do
tmp := mat[i, col]; mat[i, col] := mat[i, c]; mat[i, c] := tmp;
end_for;
end_if;
// What do we know about the selected pivot?
if isConstant(constraints, pivot) then
// 1. pivot is constant, thus we can use it for an elimination step.
beginElimination(); // ###
// Using fraction free elimination
for r from row + 1 to height do // for all rows...
tmp := mat[r, col];
if not iszero(tmp) then
for c from col + 1 to width do // ...eliminate row
old := mat[r, c]; // ###
mat[r, c] := pivot * mat[r, c]  mat[row, c] * tmp;
eliminationStep(old, mat[r, c]); // ###
end_for;
// divide row into its gccd
d := gccd(constraints, mat[r, c] $ c=col+1..width);
if degree(d) > 0 then
for c from col + 1 to width do
mat[r, c] := mat[r, c] / d;
end_for;
end_if;
mat[r, col] := zero;
end_if;
end_for;
endElimination(); // ###
row := row + 1; // next row
elif not isZero(constraints, pivot) then
// 2. pivot is neither constant nor zero, thus we have to branch on those
// two possibilities.
// We will pass the 'isConstant'case to a recursive call and do the 'isZero'
// case ourself. (Due to technical reasons, this is not easily exchangable!)
result := result . elimRec(mat, addConstant(constraints, pivot), row, col);
constraints := addZero(constraints, pivot); // 'our case'
// print(constraints);
//print(constraints, matrix);
// this is not necessarily covered by simplifySubmatrix
mat[row, col] := zero;
// Simplify the matrix's entries with respect to the current constraints
simplifySubmatrix(row, col, constraints);
// There may be other nonconstant, nonzero entries in the current column,
// thus, we will visit the current column again.
col := col  1; // incremented immediately
// else
// 3. pivot is zero, thus nothing has to be done.
end_if;
col := col + 1; // go to next column
end_while; // matches "while col <= width and row <= height do ..."
// Leaf of this search path is reached, add the solution to the list of results
// At this point, you may insert code that does further calculations with the matrix,
// e.g. determining its rank, determinant, upper echolon form, etc., and put that
// value into the result's list.
// This is the only leaf of the seach tree in the whole algorithm.
// NYI: Verify results
append(result, [mat, constraints]);
end_proc:
/*
choosePivot  manages the choose of a privot element
choosePivot(constraints, row, col)
constraints  ConstriantStore object with the set of current constriants.
row  the row of the top left element of the submatrix of interest
col  the column of the top left element of the submatrix of interest
choosePivot returns the index of a proper pivot element. The pivot index is returned as a
sequence r, c containing the row index and the column index.
The selection strategie is as follows:

Select p := row, col as temporary pivot.

If the Markowitz criterion is allowed, update p to the Markowitz' suggestion.
If p is constant w.r.t. the ConstraintStore, return p as final selection.

Otherwise, try to find a strict constant entry in p's column.
If there is one, return this item as final selection.

Otherwise, if column simplifcation is allowed, simplify p's column. After that,
try again to find a nonzero entry of degree 0 in the current column.
If there is one, return this item as final selection.

Otherwise, ask the constraint store whether one of the polynomials in the current
row represents a constant. If so, return such an entry as final selection.

Otherwise, if the markowitz suggestion was not zero w.r.t the constraint store,
return that value.

Otherwise, try to find the simplest polynomial which is not zero with respect
to the constraint store. If you find one, return this item as final selection.

Otherwise, the whole column does only contain zero entries.
In this case, return an arbitrary entry.
Note that selecting a pivot may affect the matrix, if column simplification is allowed.
That's why we don't take the matrix as a parameter value but rather as a gobal variable.
The matrix is assumed to be found in the identifier mat.
*/
choosePivot
:= proc(constraints : ConstraintStore, row : DOM_INT, col : DOM_INT)
local i, j, r, c, p, t, p1, p2, n;
begin
pivotSelection(); // ###
// 1. Trivial selection
r := row; c := col;
// 2. If allowed, ask the Markowitz criterion for help
if USE_MARKOWITZ then
beginMarkowitz(); // ###
// Due to performance reasons, there are two implementations of
// the Markowitz criterion.
p := if COLUMN_SWAP_ALLOWED then
markowitz(row, col, constraints);
else
markowitzSingleColumn(row, col, constraints);
end_if;
r := p[1]; c := p[2];
p1 := r; p2 := c;
endMarkowitz(); // ###
if isConstant(constraints, mat[r, c]) then
// Markowitz has found something constant (<> zero).
// This seems to be a good choice.
markowitzSuccess(); // ###
return ((r, c));
end_if;
end_if;
// Note that the Markowitz criterion does allways select the column
// where the further search takes place.
// 3. Return a nonzero entry of p's column that has total degree 0,
// if something like this exists.
if USE_FURTHER_IMPROVEMENTS then
for r from row to height do
if _lazy_and(degree(mat[r, c]) = 0, not iszero(mat[r, c])) then
return ((r, c)); // Success!
end_if;
end_for;
for r from row to height do
if fastIsConstant(constraints, mat[r, c]) then
return ((r, c)); // Success!
end_if;
end_for;
end_if;
// 4. If allowed, simplify the selected column, and redo the previous
// selection attempt.
if _lazy_and(USE_COLUMN_SIMPLIFICATION,
nops((if iszero(mat[i, c]) then null() else 1 end_if $ i=row..height)) > 1)
then
// We don't do anything if there are only columns with leq one
// nonzero entry, because that doesn't make much sense.
beginSimplifyColumn(); // ###
p := simplifyColumn(row, col, c, constraints);
endSimplifyColumn(); // ###
if isConstant(constraints, mat[p[1], p[2]]) then
return ((p[1], p[2])); // Success!
end_if;
end_if;
// 5. Return an entry that is constant w.r.t. the constraint store, if
// something like this exists.
for r from row to height do
if isConstant(constraints, mat[r, c]) then
return ((r, c)); // Success!
end_if;
end_for;
// 6. Return the Markowitz suggestion, if Markowitz was used and has returned
// a nonzero entry
if _lazy_and(USE_MARKOWITZ, not isZero(constraints, mat[p1, p2])) then
markowitzSuccess(); // ###
return ((p1, p2));
end_if;
// 7. If there exist entries in the current column, that do not vanish
// w.r.t. to the constraint store, select the simplest of them
p := 1;
for r from row + 1 to height do
if _lazy_or(p = 1, simpler(mat[r, c], mat[p, c])) then
if isZero(constraints, mat[r, c]) then
mat[r, c] := zero;
else
p := r;
end_if;
end_if;
end_for;
// 8. It seems that the column does only contain zeros. In this case,
// return an arbitrary pivot.
if p = 1 then
p := row;
end_if:
p, c; // return what we have decided.
end_proc:
/*
markowitz  selecting a pivot element with minimum fillin
markowitz(row, col)
row  row where the submatrix starts
col  column where the submatrix starts
markowitz determines the nonzero element of the specified submatrix that produces
the least symbolic fillin and constant fillin.
*/
markowitz
:= proc(row, col, constraints)
local i, j, k, v1, v2, v3, w1, w2, w3, r, c, R, C, cat, vcat, wcat;
begin
// tools to easily determine the total number of symbolic or nonzero entries.
cat := p > if fastIsZero(constraints, p) then
0; // impossible element
elif degree(p) <> 0 then
// symbolic fillin
if fastIsConstant(constraints, p) then
2; // symbolic, but constant
else
3; // symbolic unknown
end_if
else
1; // strict constant
end_if;
// R[i, 3]: total number of strict symbolic entries in row i
// R[i, 2]: total number of symbolic entries in row i (may be pseudo symbolic constants)
// R[i, 1]: total number of nonzero entries in row i
// R[i, 0]: total number of zeros
// R[i, 4]: m  R[i, 0]  1
R := array(row..height, 0..5, (i, j) = 0 $ j = 0..3 $ i = row..height);
C := array(col..width, 0..5, (i, j) = 0 $ j = 0..3 $ i = col..width);
for i from row to height do
for j from col to width do
r := cat(mat[i, j]);
R[i, r] := R[i, r] + 1;
C[j, r] := C[j, r] + 1;
end_for;
end_for;
// do precalculations
for i from row to height do
R[i, 4] := width  col  R[i, 0] /*  1 + 1 */;
end_for;
for j from col to width do
C[j, 4] := height  row  C[j, 0] /*  1 + 1 */;
end_for;
r := row; c := col; w1 := infinity; w2 := infinity; w3 := infinity; wcat := 0;
// where is the minimum?
for j from col to width do
for i from row to height do
vcat := cat(mat[i, j]);
case vcat
of 0 do // mat[i, j] = zero, no proper pivot
v3 := infinity; v2 := infinity; v1 := infinity;
break;
of 1 do // mat[i, j] = constant
v3 := C[j, 4] * R[i, 4]  (C[j, 1]  1) * (R[i, 1]  1); // symbolic fillin
v2 := 0; // fillin of symbolic constants
v1 := (C[j, 1]  1) * (R[i, 1]  1); // constant fillin
break;
of 2 do // mat[i, j] is a symbolic constant
v3 := C[j, 4] * R[i, 4]; // symbolic fillin
v2 := C[j, 4] * R[i, 0]; // fillin of symbolic constants
v1 := 0; // constant fillin
break;
of 3 do // mat[i, j] is symbolic
v3 := C[j, 4] * (width /*  1*/  col /* + 1*/);
v2 := 0;
v1 := 0;
break;
end_case;
if _lazy_and(wcat = 0, vcat <> 0) or // everything is better than a zero
_lazy_and(wcat = 3, vcat <> 0, vcat <> 3) then // everything (except a zero)
// is better than a symbolical
r := i; c := j; w3 := v3; w2 := v2; w1 := v1; wcat := vcat;
next;
end_if;
if _lazy_and(wcat = 3, vcat = 3) then
if _lazy_or(v3 < w3, _lazy_and(v3 = w3, simpler(mat[i, j], mat[r, c]))) then
r := i; c := j; w3 := v3; // already equal: w2=v2=w1=v1=0, wcat=vcat=3
end_if;
next;
end_if;
if _lazy_or(vcat = 3, vcat = 0, v3 > w3) then next; end_if;
// now, we are comparing only constant cat's.
// is (i, j) better than the old value?
if _lazy_or(v3 < w3, /* v3=w3 and */ v2 < w2, _lazy_and(v2 = w2, v1 < w1)) then
r := i; c := j; w3 := v3; w2 := v2; w1 := v1; wcat := vcat; // yes
next;
end_if;
// the fillins are equal, choose the simpler expression
if _lazy_and(v2 = w2, v1 = w1, simpler(mat[i, j], mat[r, c])) then
r := i; wcat := vcat; c := j; // already equal: w3 := v3; w2 := v2; w1 := v1;
end_if;
end_for;
end_for;
r, c;
end_proc:
/*
markowitzSingleColumn  Markowitz criterion where column swap is not allowed.
markowitzSingleColumn(row, col)
row  top row of submatrix
col  left col of submatrix
This is just like markowitz except that you'll always get an element of
the col column as pivot. You should use this method iff COLUMN_SWAP_ALLOWED
is FALSE.
*/
markowitzSingleColumn
:= proc(row, col, constraints)
local i, j, k, v1, v2, v3, w1, w2, w3, r, c, R, C, cat, vcat, wcat;
begin
// tools to easily determine the total number of symbolic or nonzero entries.
cat := p > if _lazy_or(iszero(p), fastIsZero(constraints, p)) then
0; // impossible element
elif degree(p) <> 0 then
// symbolic fillin
if fastIsConstant(constraints, p) then
2; // symbolic, but constant
else
3; // symbolic unknown
end_if
else
1; // strict constant
end_if;
// R[i, 3]: total number of strict symbolic entries in row i
// R[i, 2]: total number of symbolic entries in row i (may be pseudo symbolic constants)
// R[i, 1]: total number of nonzero entries in row i
// R[i, 0]: total number of zeros
// R[i, 4]: m  R[i, 0]  1
R := array(row..height, 0..5, (i, j) = 0 $ j = 0..3 $ i = row..height);
C := array(col..width, 0..5, (i, j) = 0 $ j = 0..3 $ i = col..width);
for i from row to height do
for j from col to width do
r := cat(mat[i, j]);
R[i, r] := R[i, r] + 1;
C[j, r] := C[j, r] + 1;
end_for;
end_for;
// do precalculations
for i from row to height do
R[i, 4] := width  col  R[i, 0] /*  1 + 1 */;
end_for;
for j from col to width do
C[j, 4] := height  row  C[j, 0] /*  1 + 1 */;
end_for;
r := row; c := col; w1 := infinity; w2 := infinity; w3 := infinity; wcat := 0;
// where is the minimum?
j := col;
for i from row to height do
vcat := cat(mat[i, j]);
case vcat
of 0 do // mat[i, j] = zero, no proper pivot
v3 := infinity; v2 := infinity; v1 := infinity;
break;
of 1 do // mat[i, j] = constant
v3 := C[j, 4] * R[i, 4]  (C[j, 1]  1) * (R[i, 1]  1); // symbolic fillin
v2 := 0; // fillin of symbolic constants
v1 := (C[j, 1]  1) * (R[i, 1]  1); // constant fillin
break;
of 2 do // mat[i, j] is a symbolic constant
v3 := C[j, 4] * R[i, 4]; // symbolic fillin
v2 := C[j, 4] * R[i, 0]; // fillin of symbolic constants
v1 := 0; // constant fillin
break;
of 3 do // mat[i, j] is symbolic
v3 := C[j, 4] * (width /*  1*/  col /* + 1*/);
v2 := 0;
v1 := 0;
break;
end_case;
if _lazy_and(wcat = 0, vcat <> 0) or // everything is better than a zero
_lazy_and(wcat = 3, vcat <> 0, vcat <> 3) then // everything (except a zero)
// is better than a symbolical
r := i; /* c := j;*/ w3 := v3; w2 := v2; w1 := v1; wcat := vcat;
next;
end_if;
if _lazy_and(wcat = 3, vcat = 3) then
if _lazy_or(v3 < w3, _lazy_and(v3 = w3, simpler(mat[i, j], mat[r, c]))) then
r := i; /* c := j;*/ w3 := v3; // already equal: w2=v2=w1=v1=0, wcat=vcat=3
end_if;
next;
end_if;
if _lazy_or(vcat = 3, vcat = 0, v3 > w3) then next; end_if;
// now, we are comparing only constant cat's.
// is (i, j) better than the old value?
if _lazy_or(v3 < w3, v2 < w2, _lazy_and(v2 = w2, v1 < w1)) then
r := i; /*c := j;*/ w3 := v3; w2 := v2; w1 := v1; wcat := vcat; // yes
next;
end_if;
// the fillins are equal, choose the simpler expression
if _lazy_and(v2 = w2, v1 = w1, simpler(mat[i, j], mat[r, c])) then
r := i; wcat := vcat; /*c := j;*/ // already equal: w3 := v3; w2 := v2; w1 := v1;
end_if;
end_for;
r, c;
end_proc:
/*
simplifyColumn  performs row operations in order to simplify the current column
simplifyColumn(row, startCol, col)
row  the row of the top left element of the submatrix of interest
startCol  the leftmost column that has to be affected by the row operations
col  the column of the top left element of the submatrix of interest
simplifyColumn performes as few as possible row operations to the submatrix with
topleft element [row, startCol] in order to simplify the col column.
(startCol is assumed to be an index of the whole matrix, not just for the submatrix.)
The matrix of interest is assumed to be avialable at global identifier mat.
*/
simplifyColumn
:= proc(row : DOM_INT, startCol : DOM_INT, col : DOM_INT, constraints)
local column, eros, ero, e, i, j, n, s, t, other,
old, // ###
ready, factor, id, result, tmp, r, c;
begin
// ### 1st Step: Construct the graph. ###
// For each column entry, we have a list of EROs that were performed on it.
// One ERO is represented by the triple (indexOfOtherPolynomial, factor, id),
// where id is the total number of EROs performed during the algorithm.
column := array(1..height  row + 1); // contains the polynomials
eros := array(1..height  row + 1); // contains the ERO lists
for i from row to height do
column[i  row + 1] := mat[i, col];
eros[i  row + 1] := [];
end_for;
// begin to simplify
n := height  row + 1;
id := 0;
repeat
ready := TRUE;
for i from 2 to n do
for j from 1 to i  1 do
// we will not consider constants, for they're easy enough...
if degree(column[i]) = 0 or degree(column[j]) = 0 then
i := n: // to break the outer loop
break;
end_if;
// Thus, we have two nontrivial polynomials.
// do they match?
factor := lmonomial(column[i], DegInvLexOrder)
/ lmonomial(column[j], DegInvLexOrder);
if factor = FAIL then // they don't match this way.
factor := lmonomial(column[j], DegInvLexOrder)
/ lmonomial(column[i], DegInvLexOrder);
if factor <> FAIL then
// Match! j := j  factor * i simplifies j.
result := column[j]  factor * column[i];
// Is the result really simpler? (ensures 'noetherian')
//if simpler(result, column[j]) then
ready := FALSE;
column[j] := result;
id := id + 1;
eros[j] := append(eros[j], [i, factor, id]);
//end_if;
end_if;
else
// Match! i := i  factor * j simplifies i.
result := column[i]  factor * column[j];
//if simpler(result, column[i]) then
ready := FALSE;
column[i] := result;
id := id + 1;
eros[i] := append(eros[i], [j, factor, id]);
//end_if;
end_if;
end_for;
end_for;
until ready end_repeat;
// ### 2nd Step: Find the subtree of EROs to execute. ###
// first determin the simplest entry of the column (store its index into 's')
// If there are symbolic constants, prefer them.
s := 1;
ready := FALSE; // set to TRUE if a constant entry was found
for i from 2 to n do
if _lazy_and(not iszero(column[i]), simpler(column[i], column[s])) then
if not ready then
s := i;
ready := fastIsConstant(constraints, column[i]);
elif fastIsConstant(constraints, column[i]) then
s := i;
end_if;
end_if;
end_for;
r := s; // save return value
// Delete operations that do not belong to s's subtree as follows:
// An EROlist is said to be marked if the corresponding entry in columns is
// 1 or 0 (we don't need them any more, so we can easily overwrite them).
// Start with only column[s] marked, and mark (pseudorecursively) all the
// other lists whose indices occure as 'other index' in an ERO of a marked ERO list.
// When marking a list that is not yet marked, delete all the EROs of that list which
// have an id higher than the id of the ERO that invoked the marking.
// When you're ready with the propagation of one EROlist, mark it as 0 to avoid
// further propagations of that list.
// Note, that initially all entries are polynomials, thus there is no clash with
// entries that are already 1 or 0 (since they'd appear as poly(1, [...]) or poly(0,[...])
// respectively).
column[s] := 1;
repeat
ready := TRUE;
for i from 1 to n do
if column[i] = 1 then // eros[i] is marked.
column[i] := 0; // one time is enough...
// what eros where done to that column?
e := eros[i];
for j from nops(e) downto 1 do
ero := e[j];
id := ero[3]; // this is my id
other := ero[1]; // this is who served to simplify me in step 'id'.
if _lazy_and(column[other] <> 1, column[other] <> 0) then
// the other path is not yet marked.
column[other] := 1;
// delete the other's eros with id higher than my id.
while _lazy_and(nops(eros[other]) > 0,
eros[other][nops(eros[other])][3] > id) do
delete (eros[other])[nops(eros[other])];
end_while;
if nops(eros[other]) = 0 then
// no affect that was done to the other row is interesting for us
column[other] := 0;
else
// we have still to inspect who the other row came to its form
ready := FALSE;
end_if;
end_if;
end_for;
end_if;
end_for;
until ready end_repeat;
// Delete the whole ero lists of those entries that have not been marked at all.
for i from 1 to n do
if column[i] <> 0 then // 1 impossible, otherwise we were still in the previous loop
eros[i] := [];
end_if;
end_for;
// ### 3rd Step: Execute EROs and return simplified matrix. ###
// The remaining EROs have to be performed on the matrix' rows, according to their
// id order.
ero := []; // this gets the sorted list of eros to perform.
repeat
j := 1; // index of the ERO with minimum id
for i from 1 to n do
if nops(eros[i]) > 0 then
// It is sufficient to consider the first ero of each list, since further
// entries have been appended later and thus they have higher id's.
if _lazy_or(j = 1, ((eros[j])[1])[3] > ((eros[i])[1])[3]) then
// i is our new minimum
j := i;
end_if;
end_if;
end_for;
if j <> 1 then
// append the minimum to the ero list and delete it from the graph.
s := eros[j][1]; // this is the next ero.
// instead of the id, which is no longer needed, we store there the index
// of the target row.
s[3] := j;
ero := append(ero, s);
delete (eros[j])[1];
end_if;
until j = 1 end_repeat;
// Now, finally run the eros. Note that an index transformation is needed to get
// the right row indices of the original matrix.
eros := ero;
if nops(eros) > 0 then
columnSimplifySuccess();
end_if;
for i from 1 to nops(eros) do
ero := op(eros, i); // this is our current job.
factor := ero[2]; // factor (it is <>zero by construction)
s := ero[1] + row  1; // source row (as matrix index)
t := ero[3] + row  1; // target row (as matrix index)
// If the source row is marked to skip the FF divisision
// then the target row has to be marked, too.
if _lazy_and(startCol <> 1, mat[s, startCol  1] = 0) then
mat[t, startCol  1] := 0;
end_if;
// elimination step
for j from startCol to width do
if not iszero(mat[s, j]) then
old := mat[t, j]; // ###
mat[t, j] := mat[t, j]  factor * mat[s, j];
eliminationStep(old, mat[t, j]); // ###
end_if;
end_for;
end_for;
r + row  1, col; // return the index of the simplest expression.
end_proc:
/*
simpler  compares the complexity of two polynomials
simpler(p, q)
p  one polynomial
q  one polynomial
This method compares two polynomials and returns TRUE iff p is simpler
than q.
What is the simpler polynom?
If degree(p) * monomials(p) < degree(q) * monomials(q), then p is simpler.
If those values are equal and degree(p) < degree(q), then p is simpler.
If even the degrees are equal and monomials(p) < monomials(q), then p is
simpler. If these numbers are also equal and the leading monomial of p is less than
the leading monomial of q (w.r.t. DegInvLexOrder), then p is simpler.
*/
simpler
:=proc(p : DOM_POLY, q : DOM_POLY)
local dp, dq, mp, mq, a, b;
begin
dp := degree(p);
dq := degree(q);
mp := nops(op(p, 1));
mq := nops(op(q, 1));
a := dp * mp;
b := dq * mq;
if a < b then return (TRUE)
elif a > b then return (FALSE)
end_if;
if dp < dq then return (TRUE)
elif dp > dq then return (FALSE)
end_if;
if mp < mq then return (TRUE)
elif mp > mq then return (FALSE)
end_if;
bool(lterm(p) <> lterm(p + q));
end_proc:
/*
simplifySubmatrix  simplifies the submatrix with respect to the current constraints
simplifySubmatrix(row, col, constraints)
row  top row of submatrix to simplify
col  left column of submatrix to simplify
constraints  a constraint store
simplifySubmatrix substitutes matrix entries by simpler ones that are equal w.r.t.
the constraint store. If SIMPLIFY_AFTER_BRANCH is FALSE, the method does
not do anything. If an entry is siplifyable, its row is marked to be skiped at the next
FF reduction.
*/
simplifySubmatrix
:= proc(row, col, constraints)
local i, j, r;
begin
if not SIMPLIFY_AFTER_BRANCH then
return (null());
end_if;
beginMatrixSimplification(); // ###
for i from row to height do
for j from col to width do
r := simplifyPoly(constraints, mat[i, j]);
if simpler(r, mat[i, j]) then
mat[i, j] := r;
if col <> 1 then
mat[i, col  1] := 0;
end_if;
end_if;
end_for;
end_for;
endMatrixSimplification(); // ###
null();
end_proc:
/*++
COLUMN_SWAP_ALLOWED
Boolean variable to specify whether or not column swapping is allowed during the elimination
algorithm, initial value: TRUE
++*/
COLUMN_SWAP_ALLOWED := TRUE:
/*++
USE_MARKOWITZ
Boolean variable to specify whether or not the Markowitz criterion should be used while
selecting privot elements, initial value: TRUE
++*/
USE_MARKOWITZ := TRUE:
/*++
USE_COLUMN_SIMPLIFICATION
Boolean variable to specify whether or not the simplifyColumn method should be used while
selecting privot elements, initial value: TRUE
++*/
USE_COLUMN_SIMPLIFICATION := TRUE:
/*++
SIMPLIFY_AFTER_BRANCH
Boolean variable to specify whether or not the elements of the remaining submatrix should be
simplified with respect to the current constraint store immediately after a branch has taken
place, initial value: TRUE
++*/
SIMPLIFY_AFTER_BRANCH := TRUE:
/*++
USE_FURTHER_IMPROVEMENTS
Boolean variablen to specify whether or not some nondocumented simple heuristics should be
used to further improve the elimination process. Initial value: TRUE
++*/
USE_FURTHER_IMPROVEMENTS := TRUE:
// end of file Elimination.mu