// $Date: 2002/03/10 00:48:05 $ $Author: Manuel Kauers $ $Revision: 2.0 $ /** This file implements a new data type representing a set of Constraints. Thereby, a constraint is meant to be of the form polynomial = 0 or polynomial <> 0. A set of such constraints ist meant to be their logical conjungtion. Such sets are called constraint stores.
The data structure provides methods to construct constraint stores by incrementally defining polynomials to be zero or constant (i.e. non-zero). It further provides methods for reasoning about polynomials, i.e. you can determine whether a given polynomial is equal to zero or constant (or nothing at all) with respect to a particular constraint store object. This reasoning is compleet and corrent. Finally, you have the possibility to simplify a polynomial with respect to a particular constraint store object, e.g. x2 - x may be reduced to x if your constraint store contains x2 = 0. However, the simplification algorithm is not necessarily optimal, e.g. x2 = 0 implies x = 0 and thus the polynomial of the last example could have been even simplified to 0. But of course, you won't get a result that is more complicated that your input.
Note that the constraint store needs some "private" variables to handle constraints about constant polynomials or to reason about zero polynomials. Those variables are named cs_var1, cs_var2, cs_var3, ... assuming that you will hardly ever need to used such variables for your own purposes. If you want to add up to n Constraints to a store, we need n + 1 private variables.
Now, assume that we want to define the polynomial x2 - 1 to be
zero. We can achive this by stating
>> C := addZero(C, poly(a^2, [a]));
Error: Illegal argument: poly(a^2, [a]) does not have proper variables.
[ConstraintStore::new]
(Consult my Studienarbeit Ausarbeitung on http://www.kauers.de/ or inspect the code)
1. Public tools.
Apart from that, there are some "private" tools which are described below.
*/
//=============================================================================
// PART I : Definition of the data type, Constructor
//=============================================================================
ConstraintStore := newDomain("ConstraintStore"):
/*++
ConstraintStore -- Constructs a new ConstraintStore object.
ConstraintStore([varlist])
varlist - list of variables
ConstraintStore returns a new ConstraintStore object (with no constraints) over the
polynomial ring with the variables given in varlist.
You can also say ConstraintStore(polynom) which takes the variable list from polynom,
and you can give another ConstraintStore as second argument which is meant to be the
"father" of the new store, but these two possibilities are for internal use only.
++*/
ConstraintStore::new
:= proc(polynom = poly(0, [x,y,z, cs_var.i $ i=1..20]), parent = null())
local basis, dpt, result;
begin
// Convert a variable list into the zero polynomial.
if domtype(polynom) = DOM_LIST then
polynom := poly(0, polynom);
elif domtype(polynom) = DOM_EXPR then
polynom := poly(polynom);
end_if:
// Check that the first argument is (now) a polynomial.
if domtype(polynom) <> DOM_POLY then
error("Illegal argument: ".expr2text(polynom)." is neither a polynomial nor a valid " .
"variable list.");
end_if:
dpt := 0; basis := []; // depth and groebner basis partition for the new store.
if domtype(parent) = ConstraintStore then
// We are extending an existing store.
dpt := 1 + extop(parent, 4);
// Compare coefficient rings and variable lists.
//if op(poly(parent), 3) <> op(polynom, 3) then
// error("Illegal argument: ".expr2text(polynom)." has the wrong coefficient ring.");
//end_if;
if op(poly(parent), 2) <> op(polynom, 2) then
// Different variable lists. Can we lift the polynom?
if coerce(op(polynom, 2), DOM_SET) union coerce(op(poly(parent), 2), DOM_SET)
= coerce(op(poly(parent), 2), DOM_SET) then
polynom := poly(op(polynom, 1), op(poly(parent), 2));
else
error("Illegal argument: ".expr2text(polynom)." does not have proper variables.");
end_if;
end_if;
// get basis partition from cache
basis := ConstraintStore::cacheGBasis(parent, polynom);
if basis = NIL then
// Compute this stores groebner basis partition
beginBuchberger(); // ###
basis := ConstraintStore::gbasis(parent, polynom);
endBuchberger(); // ###
appendBasis(parent, polynom, basis); // put basis into cache
end_if;
else
// We are constructing an empty store.
parent := NIL;
basis := if iszero(polynom) then
[]
else
[[polynom, lterm(polynom, DegInvLexOrder),
degree(lterm(polynom, DegInvLexOrder))]]
end_if;
end_if:
result :=
new(ConstraintStore, // a ConstraintStore is a tupel, containing...
parent, // ... the father it extends (NIL if it is root)
polynom, // ... the polynomial that is defined to be zero
basis, // ... a set of polynomials that extend the parent's Groebner basis to a
// Groebner basis of the new store, and
dpt, // ... the depth of 'this' in the ConstraintStore tree.
(csid := csid + 1) // ... the individual ID for this ConstraintStore
);
// put store into cache and return
ConstraintStore::push(result);
result;
end_proc:
//=============================================================================
// PART II : (Private) Methods to access a constraint store's instance variables
//=============================================================================
// This does not look like the most elegant way to implement this...
csid := 0: // ID generator
/*--
ID -- Returns the individual ID of a constraint store
--*/
ConstraintStore::ID := C -> extop(C, 5):
/*--
parent -- Returns the ConstraintStore C extends, or NIL if C is a root
--*/
ConstraintStore::parent := C -> extop(C, 1):
parent := C -> ConstraintStore::parent(C):
/*--
poly -- Returnd the very last polynomial in C
--*/
ConstraintStore::poly := C -> extop(C, 2):
/*--
basis -- Returns the very last partition of C's incremental Groebner basis
--*/
ConstraintStore::basis := C -> extop(C, 3):
basis := C -> ConstraintStore::basis(C):
/*--
depth -- Returns the depth of C, 0 iff C is a root
--*/
ConstraintStore::depth := C -> extop(C, 4):
depth := C -> ConstraintStore::depth(C):
/*--
varlist -- Returns the variable list of the polynomials in C including the private variables
--*/
ConstraintStore::varlist := C -> op(extop(C, 2), 2):
varlist := C -> ConstraintStore::varlist(C):
/*--
polySeq -- Returns the nonzero polynomials of C and all its and all its ancestors as a sequence
--*/
ConstraintStore::polySeq := C ->
if _lazy_and(extop(C, 4) > 0, extop(C, 2) <> NIL, not iszero(extop(C, 2))) then
// print the parent's and my polynomial.
polySeq(extop(C, 1)), extop(C, 2);
elif _lazy_and(extop(C, 2) <> NIL, not iszero(extop(C, 2))) then
// print 'my' polynomial
extop(C, 2);
elif extop(C, 4) > 0 then
// print the parent's polynomials (recursively)
polySeq(extop(C, 1));
else
// this code is only reached in case of an empty store, giving usually zero.
extop(C, 2);
end_if:
polySeq := C -> ConstraintStore::polySeq(C):
/*--
totalBasis -- Return a list of all the Groebner basis partitions of C and its ancestors.
--*/
ConstraintStore::totalBasis
:= proc(C : ConstraintStore)
local tb;
begin
tb := [];
while domtype(C) = ConstraintStore do
tb := append(tb, basis(C));
C := parent(C);
end_while;
tb;
end_proc:
/*--
print -- Overwrites the default print method for printing ConstraintStores
print()
ConstraintStores are printed as sets of constraints, where constraints are of the form
polynomial = 0 or polynomial <> 0.
The polynom2constraint method is used to convert the polynomials to the above described
appearance.
--*/
ConstraintStore::print
:= proc(C : ConstraintStore)
local d, l, r;
begin
l := [p $ p in (hold(null()), polySeq(C))]: // all polynomials in a list
r := null():
for d from 1 to nops(l) do
r := (r, ConstraintStore::poly2constraint(l[d], d - 1)):
end_for;
{r};
end_proc:
/*--
poly2constraint -- "typesets" a polynomial as a Constraint expression
poly2constraint(p, d)
p - a polynomial
d - the depth of the ConstraintStore to which p belongs.
If p contains the variable cs_var.d, then p is assumed to encode an
><0-Constraint, otherwise =0-Constraint. If p is the zero polynomial, we won't
print it.
--*/
ConstraintStore::poly2constraint
:= proc(p : DOM_POLY, d : DOM_INT)
local q, v;
begin
// Assuming that p is a '\neq0-Constraint', we try to determine the original polynomial
v := op(p, 2); // This is the polynomial list
q := (p + poly(1, v)) / poly(cs_var.d, v);
if iszero(p) then
null();
elif _lazy_or(q = FAIL, degree(poly(cs_var.d, v)) <> 1) then
op(p, 1) = 0;
else
op(q, 1) <> 0;
end_if;
end_proc:
//=============================================================================
// PART III : User Inteface
//=============================================================================
/*++
addZero -- Extends an existing ConstraintStore object by an additional constraint
addZero(C, p)
C - a constraint store
p - a polynomial
addZero returns a new ConstraintStore object that extends C by the additional constraint
p = 0.
++*/
addZero
:= proc(C : ConstraintStore, p : DOM_POLY)
begin
checkPoly(p, "addZero: ".expr2text(p)); // ###
C := ConstraintStore(p, C); // this is easy...
ConstraintStore::cacheAddZero(C, p);
C;
end_proc:
/*++
addConstant -- Extends an existing ConstraintStore object by an additional constraint
addConstant(C, p)
C - a constraint store
p - a polynomial
addConstant returns a new ConstraintStore object that extends C by the additional
constraint p <> 0.
++*/
addConstant
:= proc(C : ConstraintStore, p : DOM_POLY)
local vars, q;
begin
checkPoly(p, "addConstant: ".expr2text(p)); // ###
// p <> 0 <==> p * z - 1 = 0 for a new variable z
vars := op(poly(C), 2);
// check that our private variable exists
if degree(poly(cs_var . (1 + depth(C)), vars)) <> 1 then
error("Missing private variable cs_var".expr2text(1 + depth(C)).".");
end_if;
q := p * poly(cs_var . (1 + depth(C)), vars);
if not domtype(q) = DOM_POLY then
// try to lift q to C's ring.
q := poly(op(p, 1), vars);
if q = FAIL then
error("Illegal Argument: ".expr2text(p)." is not a proper polynomial.");
end_if;
// lifting succeeded.
q := q * poly(cs_var . (1 + depth(C)), vars);
end_if;
C := ConstraintStore(q - poly(1, vars), C);
ConstraintStore::cacheAddConstant(C, p);
C;
end:
/*++
isZero -- checks whether or not the given polynomial is zero with respect to C.
isZero(C, p)
C - a constraint store
p - a polynomial
isZero returns TRUE iff C implies that p is zero.
++*/
isZero
:=proc(C : ConstraintStore, p : DOM_POLY)
local b;
begin
checkPoly(p, "isZero: ".expr2text(p)); // ###
consistencyCheck(); // ###
// trivial checks
if iszero(p) then
return (TRUE);
elif degree(p) = 0 then
return (FALSE);
end_if;
// ask the cache
b := ConstraintStore::cacheIsZero(C, p);
if b = TRUE then
return (TRUE);
elif b = FALSE then
return (FALSE);
else // b = UNKNWON
// p is zero iff p is not constant
if not isConsistent(addConstant(C, p)) then
consistencySuccess(); // ###
ConstraintStore::cacheAddZero(C, p);
TRUE;
else
ConstraintStore::cacheAddNotZero(C, p);
FALSE;
end_if;
end_if;
end_proc:
/*++
isConstant -- checks whether or not the given polynomial is constant with respect to C.
isConstant(C, p)
C - a constraint store
p - a polynomial
isConstant returns TRUE iff C implies that p is constant (i.e. never zero).
++*/
isConstant
:=proc(C : ConstraintStore, p : DOM_POLY)
local b;
begin
checkPoly(p, "isConstant: ".expr2text(p)); // ###
consistencyCheck(); // ###
// trivial checks
if iszero(p) then
return (FALSE);
elif degree(p) = 0 then
return (TRUE);
end_if;
// ask the cache
b := ConstraintStore::cacheIsConstant(C, p);
if b = TRUE then
TRUE;
elif b = FALSE then
FALSE;
else // b = UNKNOWN
// p is constant iff p is not zero
if not isConsistent(addZero(C, p)) then
consistencySuccess(); // ###
ConstraintStore::cacheAddConstant(C, p);
TRUE;
else
ConstraintStore::cacheAddNotConstant(C, p);
FALSE;
end_if;
end_if;
end_proc:
/*++
fastIsZero -- testing on zero by table lookup
fastIsZero(C, p)
C - a constraint store
p - a polynomial
fastIsZero returns TRUE if it was already calculated that p = 0
in C. This test is correct, but
Where possible, we have recycled the procedures of the &MuPAD; groebner package, for further information consult
[1] Geddes et.al. Algorithms for Computer Clgebra, Kluwer, 1992
[2] Cox et.al. Ideals, Varieties, and Algorithms, Springer, 1992
[3] Becker et.al. Groebner Bases, Springer, 1993
Throughout the implementation we use DegInvLexOrder as monomial order. */ /*-- gbasis -- Calculates a partition for an extension of a ConstraintStore gbasis(C, p) C - a ConstraintStore p - a polynomial gbasis returns a list of polynomials that extend the partitions of C and its ancestors to a Groebner basis of [C, p].
The corresponding method in &MuPAD;'s groebner package is basis. --*/ ConstraintStore::gbasis := proc(C : ConstraintStore, p : DOM_POLY) local S, G, B, P, r, h, totalBasis, monic; begin // install normalizing routine if op(p, 3) = hold(Expr) then monic := groebner::primpart; else monic := groebner::normalize; end_if; // totalBasis is a list of all the bases of C and all its ancestors. // We build a list of lists rather than one single list to avoid the creation // of too large objects. totalBasis := ConstraintStore::totalBasis(C); // lift p to an extended polynomial p := [p, lterm(p, DegInvLexOrder), degree(lterm(p, DegInvLexOrder))]; if iszero(ConstraintStore::reduce(p, totalBasis)[1]) then return ([]); end_if; // Consult Becker et.al., p. 232 or the source of groebner::basis // for information about what follows. Here, F = {p}, what makes it slightly simpler. r := ConstraintStore::update(totalBasis, [], [], p); G := op(r, 1); B := op(r, 2); // Create reducing set S := []; for P in totalBasis do for h in P do if degree(h[2]) = 0 then return([]); end_if; S := groebner::redset_insert(S, h); end_for; end_for; for h in G do if degree(h[2]) = 0 then return([h]); end_if; S := groebner::redset_insert(S, h); end_for; // do reductions while nops(B) <> 0 do h := ConstraintStore::spoly(B[1]); delete B[1]; if iszero(h[1]) then next; end_if; h := groebner::reduce(h, S, DegInvLexOrder, _plus); if iszero(h[1]) then next; end_if; h[1] := monic(h[1]); if degree(h[2]) = 0 then return([h]); end_if; r := ConstraintStore::update(totalBasis, G, B, h); G := op(r, 1); B := op(r, 2); S := groebner::redset_insert(S, h); end_while; G; end_proc: /*-- spoly -- calculates the S-polynomial of a given critical pair. spoly(c) c - a critical pair [p, q, lcm(HT(p), HT(q)), sugar] spoly returns the S-polynomial of c using groebner's s_poly --*/ ConstraintStore::spoly := c -> groebner::s_poly(c, DegInvLexOrder, _plus): /*-- reduce -- reduces a polynomial w.r.t. a given list of partial Groebner bases reduce(p, G) p - a polynomial G - a list of partitional Groebner basis reduce returns p mod <G> using the algorithm of [2, p63]. --*/ ConstraintStore::reduce := proc(p, G) local divisionOccured, i, j, b, r, lm; begin p := p[1]; r := poly(0, op(p, 2)); lm := lmonomial(p, DegInvLexOrder); while not iszero(p) do divisisionOccured := FALSE; for b in G do j := nops(b); while _lazy_and(j <> 0, not divisisionOccured) do if lm / b[j][2] <> FAIL then p := p - lm / lmonomial(b[j][1], DegInvLexOrder) * b[j][1]; lm := lmonomial(p, DegInvLexOrder); divisisionOccured := TRUE; else j := j - 1; end_if; end_while; if divisisionOccured then break; end_if; end_for; if not divisisionOccured then r := r + lm; p := p - lm; lm := lmonomial(p, DegInvLexOrder); end_if; end_while; [r, lterm(r, DegInvLexOrder), degree(lterm(r, DegInvLexOrder))]; end_proc: /*-- update -- updates a critical pair list as needed by gbasis update(totalBasis, G, B, h) totalBasis - list of compleeted Groebner basis partitions G - the part of the new partition that was calculated so far B - the current list of critical pairs h - the polynomial to append update is the reimplemenation of groebner::update for cascading Groebner bases. See [3, p230] for further information. --*/ ConstraintStore::update := proc(totalBasis, G, B, h) local lh, C, D, p, g, hit, i, j, n, F, s, lcmh; begin checkTimeout(); // ### totalBasis := append(totalBasis, G); lh := h[2]; lcmh := table(); // maps HT(p) to lcm(HT(h), HT(p)) // initialize a frequently used part of the lmch table // and create new critical pairs C := {}; for D in totalBasis do for p in D do s := groebner::term_lcm(p[2], lh); lcmh[p[2]] := s; C := _union(C, {[p, h, s]}); end_for; end_for; C := [op(C)]; // remove redundant pairs D := {}; n := nops(C); for i from 1 to n do p := C[i]; lcmp := p[3]; // = lcm(HT(h), HT(p)) hit := bool(lh * p[1][2] = lcmp); // HT(h), HT(p) are disjoint. if not hit then // if lcm(HT(h), HT(g)) does not divide lcm(HT(h), HT(p)) for all (h,g) in C bejond p // and lcm(HT(h), HT(g)) does not divide lcm(HT(h), HT(p)) for all (h,g) in D, // then it is also a hit. hit := TRUE; for j from i + 1 to n do if lcmp / C[j][3] <> FAIL then // not a hit, skip to next pair. hit := FALSE; break; end_if; end_for; if not hit then next; end_if; for g in D do if lcmp / g[3] <> FAIL then // not a hit, skip to next pair. hit := FALSE; break; end_if; end_for; if not hit then next; end_if; end_if; D := _union(D, {p}); end_for; // Remove pairs whose leading terms are disjoint. s := h[3]; // sugar of h. D := [if lh * p[1][2] <> p[3] then // if HT(h) and HT(p) are not disjoint append(p, max(p[1][3] + degree(p[3] / p[1][2]), s + degree(p[3] / lh))); else null(); end_if $ p in D]; // Remove redundant old pairs. F := [if p[3] / lh <> FAIL then // if HT(h) does not divide lcm(HT(p1), HT(p2)) g := p[1][2]; // = HT(p1) if not contains(lcmh, g) then lcmh[g] := groebner::term_lcm(g, lh); end_if; if lcmh[g] <> p[3] then // lcm(HT(p1), HT(h)) <> lcm(HT(p1), HT(p2)) g := p[2][2]; // = HT(p2) if not contains(lcmh, g) then lcmh[g] := groebner::term_lcm(g, lh); end_if; if lcmh[g] <> p[3] then // lcm(HT(p2), HT(h)) <> lcm(HT(p1), HT(p2)) // None of the conditions were met, thus p is redundant. null(); else p; end_if; // third cond. met else p; end_if; // second cond. met else p; end_if // first cond. met $ p in B]; // B := D \cup F: if nops(D) <> 0 then D := sort(D, (() -> (groebner::pair_less(args(1), args(2), DegInvLexOrder)))); // Merge D and F into B. i := 1; j := 1; B := []; while i < nops(D) and j < nops(F) do if groebner::pair_less(D[i], F[j], DegInvLexOrder) then B := append(B, D[i]); i := i + 1; else B := append(B, F[j]); j := j + 1; end_if; end_while; F := append(B, op(D, i..nops(D)), op(F, j..nops(F))); // remainder of the longer list end_if; // Update basis G := append(select(G, (() -> (args(1)[2] / lh = FAIL))), h); // Return the new basis G and the set B of critical pairs. G, F; end_proc: //============================================================================= // PART V : Caches //============================================================================= /**
There is one global cache object that is a list of cache entries. A cache entry is a list [C, a, b, c, d, e, f] with
This method is for internal use only! You may hurt compleetnes if you call this method with invalid arguments!
cacheAddZero returns null()
--*/
ConstraintStore::cacheAddZero := (C, p) -> ConstraintStore::cacheAppend(C, p, 2):
/*--
cacheAddNotZero -- inserts a new nonzero polynomial to the cache
cacheAddNotZero(C, p)
C - a ConstraintStore
p - a polynomial
Just like