// $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.

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 non-constant, non-zero 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:

  1. Select p := row, col as temporary pivot.
  2. 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.
  3. Otherwise, try to find a strict constant entry in p's column. If there is one, return this item as final selection.
  4. Otherwise, if column simplifcation is allowed, simplify p's column. After that, try again to find a non-zero entry of degree 0 in the current column. If there is one, return this item as final selection.
  5. 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.
  6. Otherwise, if the markowitz suggestion was not zero w.r.t the constraint store, return that value.
  7. 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.
  8. 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; // fill-in 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; // fill-in 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 top-left 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 ERO-list 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 (pseudo-recursively) 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 ERO-list, 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 non-documented simple heuristics should be used to further improve the elimination process. Initial value: TRUE ++*/ USE_FURTHER_IMPROVEMENTS := TRUE: // end of file Elimination.mu