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

How to deal with constraint stores?

You can easily construct a new ConstraintStore object by saying C := ConstraintStore(); {} This would give you a new empty constraint store which is extensible by adding polynomials to it. All of the polynomials you want to add to a constraint store must belong to the same ring. By default, [x, y, z, cs_var.i $ i = 1..20] is used as polynomial variables, but you can pass your own variable list to the ConstraintStore constructor.

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(x^2 - 1, [x])); {x2 - 1 = 0} This statement is possible although the given polynom is strictly speaking not a member of the right polynommial ring because [x] is a sublist of [x, y, z, cs_var.i $ i =1..20] and addZero is thus able to "lift" the polynomial to a corresponding member of the right ring. If this is not possible, an error occures:

>>  C := addZero(C, poly(a^2, [a]));
Error: Illegal argument: poly(a^2, [a]) does not have proper variables.
[ConstraintStore::new]

If you want to find out if x - 1 is zero with respect to C, just type isZero(C, poly(x - 1, [x])); FALSE and you find out that this is not necessarily the case. But if we add another constraint C := addConstant(C, poly(x + 1, [x])); {x + 1 <> 0, x2 - 1 = 0} then we will get a different result: isZero(C, poly(x - 1, [x])); TRUE since (x + 1) (x + 1) x2 - 1. One may want to define x - 1 to be constant w.r.t. C, but this would lead to an inconsistent store: C := addConstant(C, poly(x - 1, [x])): isConsistent(C); FALSE

What does internally go on?

This is a secret :-)

(Consult my Studienarbeit Ausarbeitung on http://www.kauers.de/ or inspect the code)

Tools

Here is the list of all the methods that the ConstraintStore data structure provides and that are defined within this file:

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 not compleet. ++*/ fastIsZero := proc(C : ConstraintStore, p : DOM_POLY) begin checkPoly(p, "fastIsZero: ".expr2text(p)); // ### // return (isZero(C, p)); if iszero(p) then TRUE elif degree(p) = 0 then FALSE elif ConstraintStore::cacheIsZero(C, p) = TRUE then TRUE else FALSE end_if; end_proc: /*++ fastIsConstant -- testing on constance by table lookup fastIsConstant(C, p) C - a constraint store p - a polynomial fastIsConstant return TRUE if it was already calculated that p is constant in C. This test is correct, but not compleet. ++*/ fastIsConstant := proc(C : ConstraintStore, p : DOM_POLY) begin // return (isConstant(C, p)); checkPoly(p, "fastIsConstant: ".expr2text(p)); // ### if iszero(p) then FALSE elif degree(p) = 0 then TRUE elif ConstraintStore::cacheIsConstant(C, p) = TRUE then TRUE else FALSE end_if; end_proc: /*++ isConsistent -- returns TRUE iff 1 <> 0 in C. isConsistent(C) C - a constraint store This method determines whether there exists a specialization of the variables of C that fulfill all the constraints in C. It holds that C is consistent if and only if 1<>0 with respect to C. (This is the condition implemented in the method.) ++ */ isConsistent :=proc(C : ConstraintStore) local one; begin one := poly(1, ConstraintStore::varlist(C)); not iszero(ConstraintStore::reduce([one, one], ConstraintStore::totalBasis(C))[1]); end_proc: /*++ simplifyPoly -- returns a simplified version of the given polynomial with respect to C. simplifyPoly(C, p) C - a constraint store p - a polynomial Given a constraint store C and a polynomial p, the method tries to find a polynomial q, that is equal to p with respect to C, but has a lower degree. There is no guarantee that the returned polynomial q is by any meaning as simple as possible. If no simpler polynomial that p is found, then p is returned itself. ++*/ simplifyPoly :=proc(C : ConstraintStore, p : DOM_POLY) local G, divisionOccured, i, j, b, r, lm, n, vars; begin // this method is similar to ConstraintStore::reduce([p, // lterm(p, DegInvLexOrder)], ConstraintStore::totalBasis(C))[1], // but it does only consider those polynomials that do not contain // slick variables. G := ConstraintStore::totalBasis(C); n := nops(op(p, 2)); // maximum number of slick variables 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 vars := op(poly(op(b[j][1], 1)), 2); // actural vars in b[j][1] if _lazy_and((contains(vars, cs_var.i) = 0) $ i=1..n, 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; end_proc: //============================================================================= // PART IV : Implementation of cascading Groebner bases //============================================================================= /**

Cascading Gröbner Bases

Our Groebner basis implementation is incremental, that means that a root Constraint Store has a Groebner basis for its polynomial (the polynomial itself) but if a Constraint Store has a parent, it does only need to have a partial Groebner basis, i.e. a set of polynomials that extend the Groebner basis of the ancestors to a groebner basis of the new constraint store.

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 //============================================================================= /**

Caching Intermediate Results

Another data type called Cache is provided to store intermediate results that may be used again in near future.

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

Initially, the cache is empty. An entry is appended to the cache when a constraint store is created. The caches behaviour is similar to a stack: when a constraint store c is pushed into the stack, the last item in the cache is deleted as long as there is one and it is not the parent of c. That is, the cache always represents a path in a constraint store tree. You can switch the use of the cache on and off by setting the USE_CACHE flag. */ ConstraintStore::cache := []: /*++ USE_CACHE TRUE to use the cache. TRUE is the default value. ++*/ USE_CACHE := TRUE: /*-- push -- pushes a new item to the stack push(C) C - a constraint store if C has no parent or the cache is empty, the cache is set to [[C,values]]. Otherwise, all the last entries are popped up to the parent of C (or the empty cache). --*/ ConstraintStore::push := proc(C : ConstraintStore) local c, entry, cache; begin if not USE_CACHE then return (null()) end_if; cache := ConstraintStore::cache; entry := [C, [], [], [], [], []]; c := ConstraintStore::indexOf(parent(C)); if _lazy_or(nops(cache) = 0, ConstraintStore::depth(C) = 0, c <= 0)then // compleetly new cache ConstraintStore::cache := [entry]; else // 1..c is reusable ConstraintStore::cache := [cache[i] $ i=1..c, entry]; end_if; null(); end_proc: /*-- cacheIsZero -- decides isZero by table lookup cacheIsZero(C, p) C - a constraint store p - a polynomial cacheIsZero returns TRUE if p is zero w.r.t. C and this fact is known to the cache, FALSE if it is known that p is not zero and UNKNOWN otherwise. This method is fast but far from compleet whereas isZero is compleet but far from fast. --*/ ConstraintStore::cacheIsZero := proc(C : ConstraintStore, p : DOM_POLY) local i, c; begin if not USE_CACHE then return (UNKNOWN) end_if; i := ConstraintStore::indexOf(C); c := ConstraintStore::cache; if i <= 0 or nops(ConstraintStore::cache) = 0 then cacheMiss(); // ### return (UNKNOWN) end_if; if _lazy_or(ConstraintStore::containsFactor(c[j][2], p) $ j=1..i) then cacheHit(); // ### TRUE elif _lazy_or(ConstraintStore::contains(c[i][3], p), // not zero ConstraintStore::containsAllFactors(c[j][4], p) $ j=1..i) then // constant => not zero cacheHit(); // ### FALSE else cacheMiss(); // ### UNKNOWN end_if; end_proc: /*-- cacheAddZero -- inserts a new zero polynomial to the cache cacheAddZero(C, p) C - a constraint store p - a polynomial If C is currently present in the cache (see push), p is add to C's list of polynomials that are known to be zero.

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 cacheAddZero, but for polynomials that are know not to be zero. --*/ ConstraintStore::cacheAddNotZero := (C, p) -> ConstraintStore::cacheAppend(C, p, 3): /*-- cacheIsConstant -- decides isConstant by table lookup cacheIsConstant(C, p) C - a constraint store p - a polynomial This is like cacheIsZero, but for "constant". --*/ ConstraintStore::cacheIsConstant := proc(C : ConstraintStore, p : DOM_POLY) local i, c; begin if not USE_CACHE then return (UNKNOWN) end_if; i := ConstraintStore::indexOf(C); c := ConstraintStore::cache; if i <= 0 or nops(ConstraintStore::cache) = 0 then cacheMiss(); // ### return (UNKNOWN) end_if; if _lazy_or(ConstraintStore::containsAllFactors(c[j][4], p) $ j=1..i) then cacheHit(); // ### TRUE elif _lazy_or(ConstraintStore::contains(c[i][5], p), // not constant ConstraintStore::containsFactor(c[j][2], p) $ j=1..i) then // 0 => not const cacheHit(); // ### FALSE else cacheMiss(); // ### UNKNOWN end_if; end_proc: /*-- cacheAddConstant -- inserts a new constant polynomial to the cache cacheAddConstant(C, p) C - a ConstraintStore p - a polynomial Just like cacheAddZero, but for polynomials that are know to be constant. --*/ ConstraintStore::cacheAddConstant := (C, p) -> ConstraintStore::cacheAppend(C, p, 4): /*-- cacheAddNotConstant -- inserts a new non-constant polynomial to the cache cacheAddNotConstant(C, p) C - a ConstraintStore p - a polynomial Just like cacheAddZero, but for polynomials that are know not to be constant. --*/ ConstraintStore::cacheAddNotConstant := (C, p) -> ConstraintStore::cacheAppend(C, p, 5): /*-- appendBasis -- inserts a new Gröbner basis partition appendBasis(C, p, b); C - a ConstraintStore p - a polynomial b - a Gröbner basis partition b should be the Gröbner basis partition for the store that is created when p = 0 is added to C. --*/ ConstraintStore::appendBasis := proc(C : ConstraintStore, p : DOM_POLY, b : DOM_LIST) local i, c, e; begin if not USE_CACHE then return (null()) end_if; i := ConstraintStore::indexOf(C); if i <= 0 then return (null()); end_if; e := ConstraintStore::cache[i][6]; p := simplifyPoly(C, p); // this should be more normalized -> more matches if _lazy_and(p <> e[j][1] $ j=1..nops(e)) then // p is not yet here ConstraintStore::cache[i] := append(e, [p, b]); end_if; null(); end_proc: /*-- contains -- determines whether or not a polynom is member of a list up to a constant factor contains(l, p) l - a list of polynomials p - a polynomial contains returns TRUE if there is some q in l and some constant r such that p = q * r and FALSE otherwise. This is an internal method used by cacheIsZero and cacheIsConstant. --*/ ConstraintStore::contains := proc(l : DOM_LIST, p : DOM_POLY) local i, q; begin for i from nops(l) downto 1 do q := l[i]; if _lazy_and(lterm(p) = lterm(q), nops(op(p, 1)) = nops(op(q, 1)), q / p <> FAIL) then return (TRUE) end_if; end_for; FALSE; end_proc: /*-- containsFactor - tests if one of the polynomials of a list divides p containsFactor(l, p) l - a list of polynomials p - a polynomial return TRUE if and only if there is a q in l such that p/q is a polynomial. --*/ ConstraintStore::containsFactor := proc(l : DOM_LIST, p : DOM_POLY) local i, q; begin for i from nops(l) downto 1 do q := l[i]; if p / q <> FAIL then return (TRUE); end_if; end_for; FALSE; end_proc: /*-- containsAllFactors - tests if p is a product of some of a lists polynomials containsAllFactors(l, p) l - a list of polynomials p - a polynomial containsAllFactors returns TRUE if and only if p=_mult(l[i]^e[i] $ i=1..nops(l)) where e is a list of nonnegative integers. No factorisation is done. --*/ ConstraintStore::containsAllFactors := proc(l : DOM_LIST, p : DOM_POLY) local facs; begin //print(p, l); if _lazy_or(degree(p) = 0, has(l, p)) then TRUE else facs := [if op(p, 2) <> op(q, 2) then null() else p / q end_if $ q in l]; facs := [if facs[i] = FAIL then null() else facs[i] end_if $ i=1..nops(facs)]; if nops(facs) = 0 then FALSE else _lazy_and(ConstraintStore::containsAllFactors(l, facs[i]) $ i=1..nops(facs)); end_if end_if; end_proc: /*-- cacheGBasis -- "calculates" a partitional Gröbner basis by table lookup cacheGBasis(C, p) C - a constraint store p - a polynomial If the cache list currently contains a partitional basis for adding p to C, this partition is returend, otherwise, NIL indicates a cache miss. --*/ ConstraintStore::cacheGBasis := proc(C : ConstraintStore, p : DOM_POLY) local i, e; begin if not USE_CACHE then return (NIL) end_if; i := ConstraintStore::indexOf(C); if i <= 0 then cacheMiss(); // ### return (NIL); end_if; p := simplifyPoly(C, p); for e in ConstraintStore::cache[i][6] do if e[1] = p then cacheHit(); // ### return (e[2]); end_if; end_for; cacheMiss(); // ### NIL; end_proc: /*-- cacheAppend -- internal method to store data into the cache cacheAppend(C, p, j) C - a ConstraintStore p - a polynomial j - an index cacheAppend is an abstraction of cacheAddZero, cacheAddNotZero, cacheAddConstant, cacheAddNotConstant, cacheAddZeroBasis, cacheAddConstantBasis. If the cache currently contains an entry for C, p is appended to the cache entry at index j. --*/ ConstraintStore::cacheAppend := proc(C : ConstraintStore, p : DOM_POLY, j : DOM_INT) local i; begin if not USE_CACHE then return (null()) end_if; i := ConstraintStore::indexOf(C); // We won't check if p is already present -- this was done in advance if i > 0 then ConstraintStore::cache[i][j] := append(ConstraintStore::cache[i][j], p); end_if; null(); end_proc: /*-- indexOf -- Returns the index of a constraint store in the current cache list indexOf(C) C - a constraint store If the cache list currently contains an entry for C, then this entries index is returned. Otherwise, -1 is returned to indicate that there is no cache entry for C. --*/ ConstraintStore::indexOf := proc(C) local i, c, ID; begin if C = NIL then return (-1); end_if; // we go from back to front since frequently used stores tend to sit at the end of // the list c := ConstraintStore::cache; ID := ConstraintStore::ID(C); for i from nops(c) downto 1 do if ConstraintStore::ID(c[i][1]) = ID then return (i); end_if; end_for; -1; end_proc: // end of file ConstraintStore.mu //============================================================================= // PART VI : Tools //============================================================================= /**

Tools

The final part of the file defines useful tools to deal with polynomials and ConstraintStores */ /*++ gccd -- determines the greatest constant common dvisor gccd(C, p1, ..., pn) C - a constraint store p1,..., pn - polynomials, at least one gccd determines the greatest common divisor of p1..pn which is constant w.r.t. C. Caution: You may suffer from a long runtime, because polynomial factorisation is used. ++*/ gccd := proc() local C, polies, factors, oldFacs, newFacs, exponents, oldExps, newExps, h, i, j, k, p, rnd, nils; begin C := args(1); polies := [if iszero(args(i)) then null() else args(i) end_if $ i = 2..args(0)]; if nops(polies) = 0 then return (poly(0, op(poly(C), 2))); end_if; // determine the factorized gcd of all the polynomials p := gcd(op(polies)); h := factor(p); oldFacs := [h[2 * i] $ i=1..(nops(h)-1) / 2]; oldExps := [h[2 * i + 1] $ i=1..(nops(h)-1) / 2]; if nops(oldFacs) = 0 then return (poly(1, op(polies[1], 2))); end_if; // throw away all the polynomials, that are not constant newFacs := []; newExps := []; for i from 1 to nops(oldFacs) do p := oldFacs[i]; if _lazy_and(p <> NIL, fastIsConstant(C, p)) then newFacs := append(newFacs, p); newExps := append(newExps, oldExps[i]); end_if; end_for; //print(C, polies); _mult(if newFacs[i] = NIL then null() else newFacs[i]^newExps[i] end_if $ i=1..nops(newFacs)); end_proc: