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

• ConstraintStore([vars])
Constructs a new, empty ConstraintStore object over the polynomial ring with vars as variable list.
Constructs a new constraint store where additionaly p = 0 is assumed.
Constructs a new constraint store where additionaly p is assumed to be constant.
• isZero(C, p)
Checks whether or not the given polynomial is zero with respect to C.
• isConstant(C, p)
Checks whether or not the given polynomial is constant with respect to C.
• isConsistent(C)
Returns TRUE iff there exists specializations for the variables of C's polynomial ring that fulfill all the constraints of C.

Apart from that, there are some "private" tools which are described below.

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

 Geddes et.al. Algorithms for Computer Clgebra, Kluwer, 1992
 Cox et.al. Ideals, Varieties, and Algorithms, Springer, 1992
 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)) 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) = 0 then return([]); end_if; S := groebner::redset_insert(S, h); end_for; end_for; for h in G do if degree(h) = 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); delete B; if iszero(h) then next; end_if; h := groebner::reduce(h, S, DegInvLexOrder, _plus); if iszero(h) then next; end_if; h := monic(h); if degree(h) = 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; 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] <> FAIL then p := p - lm / lmonomial(b[j], DegInvLexOrder) * b[j]; 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; 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, lh); lcmh[p] := 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; // = lcm(HT(h), HT(p)) hit := bool(lh * p = 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] <> 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 <> 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; // sugar of h. D := [if lh * p <> p then // if HT(h) and HT(p) are not disjoint append(p, max(p + degree(p / p), s + degree(p / lh))); else null(); end_if \$ p in D]; // Remove redundant old pairs. F := [if p / lh <> FAIL then // if HT(h) does not divide lcm(HT(p1), HT(p2)) g := p; // = HT(p1) if not contains(lcmh, g) then lcmh[g] := groebner::term_lcm(g, lh); end_if; if lcmh[g] <> p then // lcm(HT(p1), HT(h)) <> lcm(HT(p1), HT(p2)) g := p; // = HT(p2) if not contains(lcmh, g) then lcmh[g] := groebner::term_lcm(g, lh); end_if; if lcmh[g] <> p 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) / 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

• C a constraint store,
• a a list of polynomials that are known to be zero,
• b a list of polynomials that are known not to be zero,
• c a list of polynomials that are known to be constant,
• d a list of polynomials that are known not to be zero,
• e a list of pairs [p, b] where b is the partial gbasis when appending p = 0 to C.
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], p) \$ j=1..i) then cacheHit(); // ### TRUE elif _lazy_or(ConstraintStore::contains(c[i], p), // not zero ConstraintStore::containsAllFactors(c[j], 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!