diff options
Diffstat (limited to 'src/algebra/regset.spad.pamphlet')
| -rw-r--r-- | src/algebra/regset.spad.pamphlet | 1806 |
1 files changed, 1806 insertions, 0 deletions
diff --git a/src/algebra/regset.spad.pamphlet b/src/algebra/regset.spad.pamphlet new file mode 100644 index 00000000..99d0747e --- /dev/null +++ b/src/algebra/regset.spad.pamphlet @@ -0,0 +1,1806 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra regset.spad} +\author{Marc Moreno Maza} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{category RSETCAT RegularTriangularSetCategory} +<<category RSETCAT RegularTriangularSetCategory>>= +)abbrev category RSETCAT RegularTriangularSetCategory +++ Author: Marc Moreno Maza +++ Date Created: 09/03/1998 +++ Date Last Updated: 12/15/1998 +++ Basic Functions: +++ Related Constructors: +++ Also See: essai Graphisme +++ AMS Classifications: +++ Keywords: polynomial, multivariate, ordered variables set +++ Description: +++ The category of regular triangular sets, introduced under +++ the name regular chains in [1] (and other papers). +++ In [3] it is proved that regular triangular sets and towers of simple +++ extensions of a field are equivalent notions. +++ In the following definitions, all polynomials and ideals +++ are taken from the polynomial ring \spad{k[x1,...,xn]} where \spad{k} +++ is the fraction field of \spad{R}. +++ The triangular set \spad{[t1,...,tm]} is regular +++ iff for every \spad{i} the initial of \spad{ti+1} is invertible +++ in the tower of simple extensions associated with \spad{[t1,...,ti]}. +++ A family \spad{[T1,...,Ts]} of regular triangular sets +++ is a split of Kalkbrener of a given ideal \spad{I} +++ iff the radical of \spad{I} is equal to the intersection +++ of the radical ideals generated by the saturated ideals +++ of the \spad{[T1,...,Ti]}. +++ A family \spad{[T1,...,Ts]} of regular triangular sets +++ is a split of Kalkbrener of a given triangular set \spad{T} +++ iff it is a split of Kalkbrener of the saturated ideal of \spad{T}. +++ Let \spad{K} be an algebraic closure of \spad{k}. +++ Assume that \spad{V} is finite with cardinality +++ \spad{n} and let \spad{A} be the affine space \spad{K^n}. +++ For a regular triangular set \spad{T} let denote by \spad{W(T)} the +++ set of regular zeros of \spad{T}. +++ A family \spad{[T1,...,Ts]} of regular triangular sets +++ is a split of Lazard of a given subset \spad{S} of \spad{A} +++ iff the union of the \spad{W(Ti)} contains \spad{S} and +++ is contained in the closure of \spad{S} (w.r.t. Zariski topology). +++ A family \spad{[T1,...,Ts]} of regular triangular sets +++ is a split of Lazard of a given triangular set \spad{T} +++ if it is a split of Lazard of \spad{W(T)}. +++ Note that if \spad{[T1,...,Ts]} is a split of Lazard of +++ \spad{T} then it is also a split of Kalkbrener of \spad{T}. +++ The converse is false. +++ This category provides operations related to both kinds of +++ splits, the former being related to ideals decomposition whereas +++ the latter deals with varieties decomposition. +++ See the example illustrating the \spadtype{RegularTriangularSet} constructor +++ for more explanations about decompositions by means of regular triangular sets. \newline +++ References : +++ [1] M. KALKBRENER "Three contributions to elimination theory" +++ Phd Thesis, University of Linz, Austria, 1991. +++ [2] M. KALKBRENER "Algorithmic properties of polynomial rings" +++ Journal of Symbol. Comp. 1998 +++ [3] P. AUBRY, D. LAZARD and M. MORENO MAZA "On the Theories +++ of Triangular Sets" Journal of Symbol. Comp. (to appear) +++ [4] M. MORENO MAZA "A new algorithm for computing triangular +++ decomposition of algebraic varieties" NAG Tech. Rep. 4/98. +++ Version: 2 + +RegularTriangularSetCategory(R:GcdDomain, E:OrderedAbelianMonoidSup,_ + V:OrderedSet,P:RecursivePolynomialCategory(R,E,V)): + Category == + TriangularSetCategory(R,E,V,P) with + + purelyAlgebraic?: (P,$) -> Boolean + ++ \spad{purelyAlgebraic?(p,ts)} returns \spad{true} iff every + ++ variable of \spad{p} is algebraic w.r.t. \spad{ts}. + purelyTranscendental? : (P,$) -> Boolean + ++ \spad{purelyTranscendental?(p,ts)} returns \spad{true} iff every + ++ variable of \spad{p} is not algebraic w.r.t. \spad{ts} + algebraicCoefficients? : (P,$) -> Boolean + ++ \spad{algebraicCoefficients?(p,ts)} returns \spad{true} iff every + ++ variable of \spad{p} which is not the main one of \spad{p} + ++ is algebraic w.r.t. \spad{ts}. + purelyAlgebraic?: $ -> Boolean + ++ \spad{purelyAlgebraic?(ts)} returns true iff for every algebraic + ++ variable \spad{v} of \spad{ts} we have + ++ \spad{algebraicCoefficients?(t_v,ts_v_-)} where \spad{ts_v} + ++ is \axiomOpFrom{select}{TriangularSetCategory}(ts,v) and \spad{ts_v_-} is + ++ \axiomOpFrom{collectUnder}{TriangularSetCategory}(ts,v). + purelyAlgebraicLeadingMonomial?: (P, $) -> Boolean + ++ \spad{purelyAlgebraicLeadingMonomial?(p,ts)} returns true iff + ++ the main variable of any non-constant iterarted initial + ++ of \spad{p} is algebraic w.r.t. \spad{ts}. + invertibleElseSplit? : (P,$) -> Union(Boolean,List $) + ++ \spad{invertibleElseSplit?(p,ts)} returns \spad{true} (resp. + ++ \spad{false}) if \spad{p} is invertible in the tower + ++ associated with \spad{ts} or returns a split of Kalkbrener + ++ of \spad{ts}. + invertible? : (P,$) -> List Record(val : Boolean, tower : $) + ++ \spad{invertible?(p,ts)} returns \spad{lbwt} where \spad{lbwt.i} + ++ is the result of \spad{invertibleElseSplit?(p,lbwt.i.tower)} and + ++ the list of the \spad{(lqrwt.i).tower} is a split of Kalkbrener of \spad{ts}. + invertible?: (P,$) -> Boolean + ++ \spad{invertible?(p,ts)} returns true iff \spad{p} is invertible + ++ in the tower associated with \spad{ts}. + invertibleSet: (P,$) -> List $ + ++ \spad{invertibleSet(p,ts)} returns a split of Kalkbrener of the + ++ quotient ideal of the ideal \axiom{I} by \spad{p} where \spad{I} is + ++ the radical of saturated of \spad{ts}. + lastSubResultantElseSplit: (P, P, $) -> Union(P,List $) + ++ \spad{lastSubResultantElseSplit(p1,p2,ts)} returns either + ++ \spad{g} a quasi-monic gcd of \spad{p1} and \spad{p2} w.r.t. + ++ the \spad{ts} or a split of Kalkbrener of \spad{ts}. + ++ This assumes that \spad{p1} and \spad{p2} have the same maim + ++ variable and that this variable is greater that any variable + ++ occurring in \spad{ts}. + lastSubResultant: (P, P, $) -> List Record(val : P, tower : $) + ++ \spad{lastSubResultant(p1,p2,ts)} returns \spad{lpwt} such that + ++ \spad{lpwt.i.val} is a quasi-monic gcd of \spad{p1} and \spad{p2} + ++ w.r.t. \spad{lpwt.i.tower}, for every \spad{i}, and such + ++ that the list of the \spad{lpwt.i.tower} is a split of Kalkbrener of + ++ \spad{ts}. Moreover, if \spad{p1} and \spad{p2} do not + ++ have a non-trivial gcd w.r.t. \spad{lpwt.i.tower} then \spad{lpwt.i.val} + ++ is the resultant of these polynomials w.r.t. \spad{lpwt.i.tower}. + ++ This assumes that \spad{p1} and \spad{p2} have the same maim + ++ variable and that this variable is greater that any variable + ++ occurring in \spad{ts}. + squareFreePart: (P,$) -> List Record(val : P, tower : $) + ++ \spad{squareFreePart(p,ts)} returns \spad{lpwt} such that + ++ \spad{lpwt.i.val} is a square-free polynomial + ++ w.r.t. \spad{lpwt.i.tower}, this polynomial being associated with \spad{p} + ++ modulo \spad{lpwt.i.tower}, for every \spad{i}. Moreover, + ++ the list of the \spad{lpwt.i.tower} is a split + ++ of Kalkbrener of \spad{ts}. + ++ WARNING: This assumes that \spad{p} is a non-constant polynomial such that + ++ if \spad{p} is added to \spad{ts}, then the resulting set is a + ++ regular triangular set. + intersect: (P,$) -> List $ + ++ \spad{intersect(p,ts)} returns the same as + ++ \spad{intersect([p],ts)} + intersect: (List P, $) -> List $ + ++ \spad{intersect(lp,ts)} returns \spad{lts} a split of Lazard + ++ of the intersection of the affine variety associated + ++ with \spad{lp} and the regular zero set of \spad{ts}. + intersect: (List P, List $) -> List $ + ++ \spad{intersect(lp,lts)} returns the same as + ++ \spad{concat([intersect(lp,ts) for ts in lts])|} + intersect: (P, List $) -> List $ + ++ \spad{intersect(p,lts)} returns the same as + ++ \spad{intersect([p],lts)} + augment: (P,$) -> List $ + ++ \spad{augment(p,ts)} assumes that \spad{p} is a non-constant + ++ polynomial whose main variable is greater than any variable + ++ of \spad{ts}. This operation assumes also that if \spad{p} is + ++ added to \spad{ts} the resulting set, say \spad{ts+p}, is a + ++ regular triangular set. Then it returns a split of Kalkbrener + ++ of \spad{ts+p}. This may not be \spad{ts+p} itself, if for + ++ instance \spad{ts+p} is required to be square-free. + augment: (P,List $) -> List $ + ++ \spad{augment(p,lts)} returns the same as + ++ \spad{concat([augment(p,ts) for ts in lts])} + augment: (List P,$) -> List $ + ++ \spad{augment(lp,ts)} returns \spad{ts} if \spad{empty? lp}, + ++ \spad{augment(p,ts)} if \spad{lp = [p]}, otherwise + ++ \spad{augment(first lp, augment(rest lp, ts))} + augment: (List P,List $) -> List $ + ++ \spad{augment(lp,lts)} returns the same as + ++ \spad{concat([augment(lp,ts) for ts in lts])} + internalAugment: (P, $) -> $ + ++ \spad{internalAugment(p,ts)} assumes that \spad{augment(p,ts)} + ++ returns a singleton and returns it. + internalAugment: (List P, $) -> $ + ++ \spad{internalAugment(lp,ts)} returns \spad{ts} if \spad{lp} + ++ is empty otherwise returns + ++ \spad{internalAugment(rest lp, internalAugment(first lp, ts))} + extend: (P,$) -> List $ + ++ \spad{extend(p,ts)} assumes that \spad{p} is a non-constant + ++ polynomial whose main variable is greater than any variable + ++ of \spad{ts}. Then it returns a split of Kalkbrener + ++ of \spad{ts+p}. This may not be \spad{ts+p} itself, if for + ++ instance \spad{ts+p} is not a regular triangular set. + extend: (P, List $) -> List $ + ++ \spad{extend(p,lts)} returns the same as + ++ \spad{concat([extend(p,ts) for ts in lts])|} + extend: (List P,$) -> List $ + ++ \spad{extend(lp,ts)} returns \spad{ts} if \spad{empty? lp} + ++ \spad{extend(p,ts)} if \spad{lp = [p]} else + ++ \spad{extend(first lp, extend(rest lp, ts))} + extend: (List P,List $) -> List $ + ++ \spad{extend(lp,lts)} returns the same as + ++ \spad{concat([extend(lp,ts) for ts in lts])|} + zeroSetSplit: (List P, Boolean) -> List $ + ++ \spad{zeroSetSplit(lp,clos?)} returns \spad{lts} a split of Kalkbrener + ++ of the radical ideal associated with \spad{lp}. + ++ If \spad{clos?} is false, it is also a decomposition of the + ++ variety associated with \spad{lp} into the regular zero set of the \spad{ts} in \spad{lts} + ++ (or, in other words, a split of Lazard of this variety). + ++ See the example illustrating the \spadtype{RegularTriangularSet} constructor + ++ for more explanations about decompositions by means of regular triangular sets. + + add + + NNI ==> NonNegativeInteger + INT ==> Integer + LP ==> List P + PWT ==> Record(val : P, tower : $) + LpWT ==> Record(val : (List P), tower : $) + Split ==> List $ + pack ==> PolynomialSetUtilitiesPackage(R,E,V,P) + + purelyAlgebraic?(p: P, ts: $): Boolean == + ground? p => true + not algebraic?(mvar(p),ts) => false + algebraicCoefficients?(p,ts) + + purelyTranscendental?(p:P,ts:$): Boolean == + empty? ts => true + lv : List V := variables(p)$P + while (not empty? lv) and (not algebraic?(first(lv),ts)) repeat lv := rest lv + empty? lv + + purelyAlgebraicLeadingMonomial?(p: P, ts: $): Boolean == + ground? p => true + algebraic?(mvar(p),ts) and purelyAlgebraicLeadingMonomial?(init(p), ts) + + algebraicCoefficients?(p:P,ts:$): Boolean == + ground? p => true + (not ground? init(p)) and not (algebraic?(mvar(init(p)),ts)) => false + algebraicCoefficients?(init(p),ts) => + ground? tail(p) => true + mvar(tail(p)) = mvar(p) => + algebraicCoefficients?(tail(p),ts) + algebraic?(mvar(tail(p)),ts) => + algebraicCoefficients?(tail(p),ts) + false + false + + if V has Finite + then + purelyAlgebraic?(ts: $): Boolean == + empty? ts => true + size()$V = #ts => true + lp: LP := sort(infRittWu?,members(ts)) + i: NonNegativeInteger := size()$V + for p in lp repeat + v: V := mvar(p) + (i = (lookup(v)$V)::NNI) => + i := subtractIfCan(i,1)::NNI + univariate?(p)$pack => + i := subtractIfCan(i,1)::NNI + not algebraicCoefficients?(p,collectUnder(ts,v)) => + return false + i := subtractIfCan(i,1)::NNI + true + + else + + purelyAlgebraic?(ts: $): Boolean == + empty? ts => true + v: V := mvar(ts) + p: P := select(ts,v)::P + ts := collectUnder(ts,v) + empty? ts => univariate?(p)$pack + not purelyAlgebraic?(ts) => false + algebraicCoefficients?(p,ts) + + augment(p:P,lts:List $) == + toSave: Split := [] + while not empty? lts repeat + ts := first lts + lts := rest lts + toSave := concat(augment(p,ts),toSave) + toSave + + augment(lp:LP,ts:$) == + toSave: Split := [ts] + empty? lp => toSave + lp := sort(infRittWu?,lp) + while not empty? lp repeat + p := first lp + lp := rest lp + toSave := augment(p,toSave) + toSave + + augment(lp:LP,lts:List $) == + empty? lp => lts + toSave: Split := [] + while not empty? lts repeat + ts := first lts + lts := rest lts + toSave := concat(augment(lp,ts),toSave) + toSave + + extend(p:P,lts:List $) == + toSave : Split := [] + while not empty? lts repeat + ts := first lts + lts := rest lts + toSave := concat(extend(p,ts),toSave) + toSave + + extend(lp:LP,ts:$) == + toSave: Split := [ts] + empty? lp => toSave + lp := sort(infRittWu?,lp) + while not empty? lp repeat + p := first lp + lp := rest lp + toSave := extend(p,toSave) + toSave + + extend(lp:LP,lts:List $) == + empty? lp => lts + toSave: Split := [] + while not empty? lts repeat + ts := first lts + lts := rest lts + toSave := concat(extend(lp,ts),toSave) + toSave + + intersect(lp:LP,lts:List $): List $ == + -- A VERY GENERAL default algorithm + (empty? lp) or (empty? lts) => lts + lp := [primitivePart(p) for p in lp] + lp := removeDuplicates lp + lp := remove(zero?,lp) + any?(ground?,lp) => [] + toSee: List LpWT := [[lp,ts]$LpWT for ts in lts] + toSave: List $ := [] + lp: LP + p: P + ts: $ + lus: List $ + while (not empty? toSee) repeat + lpwt := first toSee; toSee := rest toSee + lp := lpwt.val; ts := lpwt.tower + empty? lp => toSave := cons(ts, toSave) + p := first lp; lp := rest lp + lus := intersect(p,ts) + toSee := concat([[lp,us]$LpWT for us in lus], toSee) + toSave + + intersect(lp: LP,ts: $): List $ == + intersect(lp,[ts]) + + intersect(p: P,lts: List $): List $ == + intersect([p],lts) + +@ +\section{package QCMPACK QuasiComponentPackage} +<<package QCMPACK QuasiComponentPackage>>= +)abbrev package QCMPACK QuasiComponentPackage +++ Author: Marc Moreno Maza +++ marc@nag.co.uk +++ Date Created: 08/30/1998 +++ Date Last Updated: 12/16/1998 +++ Basic Functions: +++ Related Constructors: +++ Also See: `tosedom.spad' +++ AMS Classifications: +++ Keywords: +++ Description: +++ A package for removing redundant quasi-components and redundant +++ branches when decomposing a variety by means of quasi-components +++ of regular triangular sets. \newline +++ References : +++ [1] D. LAZARD "A new method for solving algebraic systems of +++ positive dimension" Discr. App. Math. 33:147-160,1991 +++ [2] M. MORENO MAZA "Calculs de pgcd au-dessus des tours +++ d'extensions simples et resolution des systemes d'equations +++ algebriques" These, Universite P.etM. Curie, Paris, 1997. +++ [3] M. MORENO MAZA "A new algorithm for computing triangular +++ decomposition of algebraic varieties" NAG Tech. Rep. 4/98. +++ Version: 3. + +QuasiComponentPackage(R,E,V,P,TS): Exports == Implementation where + + R : GcdDomain + E : OrderedAbelianMonoidSup + V : OrderedSet + P : RecursivePolynomialCategory(R,E,V) + TS : RegularTriangularSetCategory(R,E,V,P) + N ==> NonNegativeInteger + Z ==> Integer + B ==> Boolean + S ==> String + LP ==> List P + PtoP ==> P -> P + PS ==> GeneralPolynomialSet(R,E,V,P) + PWT ==> Record(val : P, tower : TS) + BWT ==> Record(val : Boolean, tower : TS) + LpWT ==> Record(val : (List P), tower : TS) + Branch ==> Record(eq: List P, tower: TS, ineq: List P) + UBF ==> Union(Branch,"failed") + Split ==> List TS + Key ==> Record(left:TS, right:TS) + Entry ==> Boolean + H ==> TabulatedComputationPackage(Key, Entry) + polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P) + + Exports == with + startTable!: (S,S,S) -> Void + ++ \axiom{startTableGcd!(s1,s2,s3)} + ++ is an internal subroutine, exported only for developement. + stopTable!: () -> Void + ++ \axiom{stopTableGcd!()} + ++ is an internal subroutine, exported only for developement. + supDimElseRittWu?: (TS,TS) -> Boolean + ++ \axiom{supDimElseRittWu(ts,us)} returns true iff \axiom{ts} + ++ has less elements than \axiom{us} otherwise if \axiom{ts} + ++ has higher rank than \axiom{us} w.r.t. Riit and Wu ordering. + algebraicSort: Split -> Split + ++ \axiom{algebraicSort(lts)} sorts \axiom{lts} w.r.t + ++ \axiomOpFrom{supDimElseRittWu?}{QuasiComponentPackage}. + moreAlgebraic?: (TS,TS) -> Boolean + ++ \axiom{moreAlgebraic?(ts,us)} returns false iff \axiom{ts} + ++ and \axiom{us} are both empty, or \axiom{ts} + ++ has less elements than \axiom{us}, or some variable is + ++ algebraic w.r.t. \axiom{us} and is not w.r.t. \axiom{ts}. + subTriSet?: (TS,TS) -> Boolean + ++ \axiom{subTriSet?(ts,us)} returns true iff \axiom{ts} is + ++ a sub-set of \axiom{us}. + subPolSet?: (LP, LP) -> Boolean + ++ \axiom{subPolSet?(lp1,lp2)} returns true iff \axiom{lp1} is + ++ a sub-set of \axiom{lp2}. + internalSubPolSet?: (LP, LP) -> Boolean + ++ \axiom{internalSubPolSet?(lp1,lp2)} returns true iff \axiom{lp1} is + ++ a sub-set of \axiom{lp2} assuming that these lists are sorted + ++ increasingly w.r.t. \axiomOpFrom{infRittWu?}{RecursivePolynomialCategory}. + internalInfRittWu?: (LP, LP) -> Boolean + ++ \axiom{internalInfRittWu?(lp1,lp2)} + ++ is an internal subroutine, exported only for developement. + infRittWu?: (LP, LP) -> Boolean + ++ \axiom{infRittWu?(lp1,lp2)} + ++ is an internal subroutine, exported only for developement. + internalSubQuasiComponent?: (TS,TS) -> Union(Boolean,"failed") + ++ \axiom{internalSubQuasiComponent?(ts,us)} returns a boolean \spad{b} value + ++ if the fact that the regular zero set of \axiom{us} contains that of + ++ \axiom{ts} can be decided (and in that case \axiom{b} gives this + ++ inclusion) otherwise returns \axiom{"failed"}. + subQuasiComponent?: (TS,TS) -> Boolean + ++ \axiom{subQuasiComponent?(ts,us)} returns true iff + ++ \axiomOpFrom{internalSubQuasiComponent?}{QuasiComponentPackage} + ++ returs true. + subQuasiComponent?: (TS,Split) -> Boolean + ++ \axiom{subQuasiComponent?(ts,lus)} returns true iff + ++ \axiom{subQuasiComponent?(ts,us)} holds for one \spad{us} in \spad{lus}. + removeSuperfluousQuasiComponents: Split -> Split + ++ \axiom{removeSuperfluousQuasiComponents(lts)} removes from \axiom{lts} + ++ any \spad{ts} such that \axiom{subQuasiComponent?(ts,us)} holds for + ++ another \spad{us} in \axiom{lts}. + subCase?: (LpWT,LpWT) -> Boolean + ++ \axiom{subCase?(lpwt1,lpwt2)} + ++ is an internal subroutine, exported only for developement. + removeSuperfluousCases: List LpWT -> List LpWT + ++ \axiom{removeSuperfluousCases(llpwt)} + ++ is an internal subroutine, exported only for developement. + prepareDecompose: (LP, List(TS),B,B) -> List Branch + ++ \axiom{prepareDecompose(lp,lts,b1,b2)} + ++ is an internal subroutine, exported only for developement. + branchIfCan: (LP,TS,LP,B,B,B,B,B) -> Union(Branch,"failed") + ++ \axiom{branchIfCan(leq,ts,lineq,b1,b2,b3,b4,b5)} + ++ is an internal subroutine, exported only for developement. + + Implementation == add + + squareFreeFactors(lp: LP): LP == + lsflp: LP := [] + for p in lp repeat + lsfp := squareFreeFactors(p)$polsetpack + lsflp := concat(lsfp,lsflp) + sort(infRittWu?,removeDuplicates lsflp) + + startTable!(ok: S, ko: S, domainName: S): Void == + initTable!()$H + if (not empty? ok) and (not empty? ko) then printInfo!(ok,ko)$H + if (not empty? domainName) then startStats!(domainName)$H + void() + + stopTable!(): Void == + if makingStats?()$H then printStats!()$H + clearTable!()$H + + supDimElseRittWu? (ts:TS,us:TS): Boolean == + #ts < #us => true + #ts > #us => false + lp1 :LP := members(ts) + lp2 :LP := members(us) + while (not empty? lp1) and (not infRittWu?(first(lp2),first(lp1))) repeat + lp1 := rest lp1 + lp2 := rest lp2 + not empty? lp1 + + algebraicSort (lts:Split): Split == + lts := removeDuplicates lts + sort(supDimElseRittWu?,lts) + + moreAlgebraic?(ts:TS,us:TS): Boolean == + empty? ts => empty? us + empty? us => true + #ts < #us => false + for p in (members us) repeat + not algebraic?(mvar(p),ts) => return false + true + + subTriSet?(ts:TS,us:TS): Boolean == + empty? ts => true + empty? us => false + mvar(ts) > mvar(us) => false + mvar(ts) < mvar(us) => subTriSet?(ts,rest(us)::TS) + first(ts)::P = first(us)::P => subTriSet?(rest(ts)::TS,rest(us)::TS) + false + + internalSubPolSet?(lp1: LP, lp2: LP): Boolean == + empty? lp1 => true + empty? lp2 => false + associates?(first lp1, first lp2) => + internalSubPolSet?(rest lp1, rest lp2) + infRittWu?(first lp1, first lp2) => false + internalSubPolSet?(lp1, rest lp2) + + subPolSet?(lp1: LP, lp2: LP): Boolean == + lp1 := sort(infRittWu?, lp1) + lp2 := sort(infRittWu?, lp2) + internalSubPolSet?(lp1,lp2) + + infRittWu?(lp1: LP, lp2: LP): Boolean == + lp1 := sort(infRittWu?, lp1) + lp2 := sort(infRittWu?, lp2) + internalInfRittWu?(lp1,lp2) + + internalInfRittWu?(lp1: LP, lp2: LP): Boolean == + empty? lp1 => not empty? lp2 + empty? lp2 => false + infRittWu?(first lp1, first lp2)$P => true + infRittWu?(first lp2, first lp1)$P => false + infRittWu?(rest lp1, rest lp2)$$ + + subCase? (lpwt1:LpWT,lpwt2:LpWT): Boolean == + -- ASSUME lpwt.{1,2}.val is sorted w.r.t. infRittWu? + not internalSubPolSet?(lpwt2.val, lpwt1.val) => false + subQuasiComponent?(lpwt1.tower,lpwt2.tower) + + internalSubQuasiComponent?(ts:TS,us:TS): Union(Boolean,"failed") == + -- "failed" is false iff saturate(us) is radical + subTriSet?(us,ts) => true + not moreAlgebraic?(ts,us) => false::Union(Boolean,"failed") + for p in (members us) repeat + mdeg(p) < mdeg(select(ts,mvar(p))::P) => + return("failed"::Union(Boolean,"failed")) + for p in (members us) repeat + not zero? initiallyReduce(p,ts) => + return("failed"::Union(Boolean,"failed")) + lsfp := squareFreeFactors(initials us) + for p in lsfp repeat + not invertible?(p,ts)@B => + return(false::Union(Boolean,"failed")) + true::Union(Boolean,"failed") + + subQuasiComponent?(ts:TS,us:TS): Boolean == + k: Key := [ts, us] + e := extractIfCan(k)$H + e case Entry => e::Entry + ubf: Union(Boolean,"failed") := internalSubQuasiComponent?(ts,us) + b: Boolean := (ubf case Boolean) and (ubf::Boolean) + insert!(k,b)$H + b + + subQuasiComponent?(ts:TS,lus:Split): Boolean == + for us in lus repeat + subQuasiComponent?(ts,us)@B => return true + false + + removeSuperfluousCases (cases:List LpWT) == + #cases < 2 => cases + toSee := sort(supDimElseRittWu?(#1.tower,#2.tower),cases) + lpwt1,lpwt2 : LpWT + toSave,headmaxcases,maxcases,copymaxcases : List LpWT + while not empty? toSee repeat + lpwt1 := first toSee + toSee := rest toSee + toSave := [] + for lpwt2 in toSee repeat + if subCase?(lpwt1,lpwt2) + then + lpwt1 := lpwt2 + else + if not subCase?(lpwt2,lpwt1) + then + toSave := cons(lpwt2,toSave) + if empty? maxcases + then + headmaxcases := [lpwt1] + maxcases := headmaxcases + else + copymaxcases := maxcases + while (not empty? copymaxcases) and _ + (not subCase?(lpwt1,first(copymaxcases))) repeat + copymaxcases := rest copymaxcases + if empty? copymaxcases + then + setrest!(headmaxcases,[lpwt1]) + headmaxcases := rest headmaxcases + toSee := reverse toSave + maxcases + + removeSuperfluousQuasiComponents(lts: Split): Split == + lts := removeDuplicates lts + #lts < 2 => lts + toSee := algebraicSort lts + toSave,headmaxlts,maxlts,copymaxlts : Split + while not empty? toSee repeat + ts := first toSee + toSee := rest toSee + toSave := [] + for us in toSee repeat + if subQuasiComponent?(ts,us)@B + then + ts := us + else + if not subQuasiComponent?(us,ts)@B + then + toSave := cons(us,toSave) + if empty? maxlts + then + headmaxlts := [ts] + maxlts := headmaxlts + else + copymaxlts := maxlts + while (not empty? copymaxlts) and _ + (not subQuasiComponent?(ts,first(copymaxlts))@B) repeat + copymaxlts := rest copymaxlts + if empty? copymaxlts + then + setrest!(headmaxlts,[ts]) + headmaxlts := rest headmaxlts + toSee := reverse toSave + algebraicSort maxlts + + removeAssociates (lp:LP):LP == + removeDuplicates [primitivePart(p) for p in lp] + + branchIfCan(leq: LP,ts: TS,lineq: LP, b1:B,b2:B,b3:B,b4:B,b5:B):UBF == + -- ASSUME pols in leq are squarefree and mainly primitive + -- if b1 then CLEAN UP leq + -- if b2 then CLEAN UP lineq + -- if b3 then SEARCH for ZERO in lineq with leq + -- if b4 then SEARCH for ZERO in lineq with ts + -- if b5 then SEARCH for ONE in leq with lineq + if b1 + then + leq := removeAssociates(leq) + leq := remove(zero?,leq) + any?(ground?,leq) => + return("failed"::Union(Branch,"failed")) + if b2 + then + any?(zero?,lineq) => + return("failed"::Union(Branch,"failed")) + lineq := removeRedundantFactors(lineq)$polsetpack + if b3 + then + ps: PS := construct(leq)$PS + for q in lineq repeat + zero? remainder(q,ps).polnum => + return("failed"::Union(Branch,"failed")) + (empty? leq) or (empty? lineq) => ([leq, ts, lineq]$Branch)::UBF + if b4 + then + for q in lineq repeat + zero? initiallyReduce(q,ts) => + return("failed"::Union(Branch,"failed")) + if b5 + then + newleq: LP := [] + for p in leq repeat + for q in lineq repeat + if mvar(p) = mvar(q) + then + g := gcd(p,q) + newp := (p exquo g)::P + ground? newp => + return("failed"::Union(Branch,"failed")) + newleq := cons(newp,newleq) + else + newleq := cons(p,newleq) + leq := newleq + leq := sort(infRittWu?, removeDuplicates leq) + ([leq, ts, lineq]$Branch)::UBF + + prepareDecompose(lp: LP, lts: List(TS), b1: B, b2: B): List Branch == + -- if b1 then REMOVE REDUNDANT COMPONENTS in lts + -- if b2 then SPLIT the input system with squareFree + lp := sort(infRittWu?, remove(zero?,removeAssociates(lp))) + any?(ground?,lp) => [] + empty? lts => [] + if b1 then lts := removeSuperfluousQuasiComponents lts + not b2 => + [[lp,ts,squareFreeFactors(initials ts)]$Branch for ts in lts] + toSee: List Branch + lq: LP := [] + toSee := [[lq,ts,squareFreeFactors(initials ts)]$Branch for ts in lts] + empty? lp => toSee + for p in lp repeat + lsfp := squareFreeFactors(p)$polsetpack + branches: List Branch := [] + lq := [] + for f in lsfp repeat + for branch in toSee repeat + leq : LP := branch.eq + ts := branch.tower + lineq : LP := branch.ineq + ubf1: UBF := branchIfCan(leq,ts,lq,false,false,true,true,true)@UBF + ubf1 case "failed" => "leave" + ubf2: UBF := branchIfCan([f],ts,lineq,false,false,true,true,true)@UBF + ubf2 case "failed" => "leave" + leq := sort(infRittWu?,removeDuplicates concat(ubf1.eq,ubf2.eq)) + lineq := sort(infRittWu?,removeDuplicates concat(ubf1.ineq,ubf2.ineq)) + newBranch := branchIfCan(leq,ts,lineq,false,false,false,false,false) + branches:= cons(newBranch::Branch,branches) + lq := cons(f,lq) + toSee := branches + sort(supDimElseRittWu?(#1.tower,#2.tower),toSee) + +@ +\section{package RSETGCD RegularTriangularSetGcdPackage} +<<package RSETGCD RegularTriangularSetGcdPackage>>= +)abbrev package RSETGCD RegularTriangularSetGcdPackage +++ Author: Marc Moreno Maza (marc@nag.co.uk) +++ Date Created: 08/30/1998 +++ Date Last Updated: 12/15/1998 +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ Description: +++ An internal package for computing gcds and resultants of univariate +++ polynomials with coefficients in a tower of simple extensions of a field.\newline +++ References : +++ [1] M. MORENO MAZA and R. RIOBOO "Computations of gcd over +++ algebraic towers of simple extensions" In proceedings of AAECC11 +++ Paris, 1995. +++ [2] M. MORENO MAZA "Calculs de pgcd au-dessus des tours +++ d'extensions simples et resolution des systemes d'equations +++ algebriques" These, Universite P.etM. Curie, Paris, 1997. +++ [3] M. MORENO MAZA "A new algorithm for computing triangular +++ decomposition of algebraic varieties" NAG Tech. Rep. 4/98. +++ Version: 4. + +RegularTriangularSetGcdPackage(R,E,V,P,TS): Exports == Implementation where + + R : GcdDomain + E : OrderedAbelianMonoidSup + V : OrderedSet + P : RecursivePolynomialCategory(R,E,V) + TS : RegularTriangularSetCategory(R,E,V,P) + N ==> NonNegativeInteger + Z ==> Integer + B ==> Boolean + S ==> String + LP ==> List P + PtoP ==> P -> P + PS ==> GeneralPolynomialSet(R,E,V,P) + PWT ==> Record(val : P, tower : TS) + BWT ==> Record(val : Boolean, tower : TS) + LpWT ==> Record(val : (List P), tower : TS) + Branch ==> Record(eq: List P, tower: TS, ineq: List P) + UBF ==> Union(Branch,"failed") + Split ==> List TS + KeyGcd ==> Record(arg1: P, arg2: P, arg3: TS, arg4: B) + EntryGcd ==> List PWT + HGcd ==> TabulatedComputationPackage(KeyGcd, EntryGcd) + KeyInvSet ==> Record(arg1: P, arg3: TS) + EntryInvSet ==> List TS + HInvSet ==> TabulatedComputationPackage(KeyInvSet, EntryInvSet) + polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P) + quasicomppack ==> QuasiComponentPackage(R,E,V,P,TS) + + Exports == with + startTableGcd!: (S,S,S) -> Void + ++ \axiom{startTableGcd!(s1,s2,s3)} + ++ is an internal subroutine, exported only for developement. + stopTableGcd!: () -> Void + ++ \axiom{stopTableGcd!()} + ++ is an internal subroutine, exported only for developement. + startTableInvSet!: (S,S,S) -> Void + ++ \axiom{startTableInvSet!(s1,s2,s3)} + ++ is an internal subroutine, exported only for developement. + stopTableInvSet!: () -> Void + ++ \axiom{stopTableInvSet!()} is an internal subroutine, + ++ exported only for developement. + prepareSubResAlgo: (P,P,TS) -> List LpWT + ++ \axiom{prepareSubResAlgo(p1,p2,ts)} + ++ is an internal subroutine, exported only for developement. + internalLastSubResultant: (P,P,TS,B,B) -> List PWT + ++ \axiom{internalLastSubResultant(p1,p2,ts,inv?,break?)} + ++ is an internal subroutine, exported only for developement. + internalLastSubResultant: (List LpWT,V,B) -> List PWT + ++ \axiom{internalLastSubResultant(lpwt,v,flag)} is an internal + ++ subroutine, exported only for developement. + integralLastSubResultant: (P,P,TS) -> List PWT + ++ \axiom{integralLastSubResultant(p1,p2,ts)} + ++ is an internal subroutine, exported only for developement. + toseLastSubResultant: (P,P,TS) -> List PWT + ++ \axiom{toseLastSubResultant(p1,p2,ts)} has the same specifications as + ++ \axiomOpFrom{lastSubResultant}{RegularTriangularSetCategory}. + toseInvertible?: (P,TS) -> B + ++ \axiom{toseInvertible?(p1,p2,ts)} has the same specifications as + ++ \axiomOpFrom{invertible?}{RegularTriangularSetCategory}. + toseInvertible?: (P,TS) -> List BWT + ++ \axiom{toseInvertible?(p1,p2,ts)} has the same specifications as + ++ \axiomOpFrom{invertible?}{RegularTriangularSetCategory}. + toseInvertibleSet: (P,TS) -> Split + ++ \axiom{toseInvertibleSet(p1,p2,ts)} has the same specifications as + ++ \axiomOpFrom{invertibleSet}{RegularTriangularSetCategory}. + toseSquareFreePart: (P,TS) -> List PWT + ++ \axiom{toseSquareFreePart(p,ts)} has the same specifications as + ++ \axiomOpFrom{squareFreePart}{RegularTriangularSetCategory}. + + Implementation == add + + startTableGcd!(ok: S, ko: S, domainName: S): Void == + initTable!()$HGcd + printInfo!(ok,ko)$HGcd + startStats!(domainName)$HGcd + void() + + stopTableGcd!(): Void == + if makingStats?()$HGcd then printStats!()$HGcd + clearTable!()$HGcd + + startTableInvSet!(ok: S, ko: S, domainName: S): Void == + initTable!()$HInvSet + printInfo!(ok,ko)$HInvSet + startStats!(domainName)$HInvSet + void() + + stopTableInvSet!(): Void == + if makingStats?()$HInvSet then printStats!()$HInvSet + clearTable!()$HInvSet + + toseInvertible?(p:P,ts:TS): Boolean == + q := primitivePart initiallyReduce(p,ts) + zero? q => false + normalized?(q,ts) => true + v := mvar(q) + not algebraic?(v,ts) => + toCheck: List BWT := toseInvertible?(p,ts)@(List BWT) + for bwt in toCheck repeat + bwt.val = false => return false + return true + ts_v := select(ts,v)::P + ts_v_- := collectUnder(ts,v) + lgwt := internalLastSubResultant(ts_v,q,ts_v_-,false,true) + for gwt in lgwt repeat + g := gwt.val; + (not ground? g) and (mvar(g) = v) => + return false + true + + toseInvertible?(p:P,ts:TS): List BWT == + q := primitivePart initiallyReduce(p,ts) + zero? q => [[false,ts]$BWT] + normalized?(q,ts) => [[true,ts]$BWT] + v := mvar(q) + not algebraic?(v,ts) => + lbwt: List BWT := [] + toCheck: List BWT := toseInvertible?(init(q),ts)@(List BWT) + for bwt in toCheck repeat + bwt.val => lbwt := cons(bwt,lbwt) + newq := removeZero(q,bwt.tower) + zero? newq => lbwt := cons(bwt,lbwt) + lbwt := concat(toseInvertible?(newq,bwt.tower)@(List BWT), lbwt) + return lbwt + ts_v := select(ts,v)::P + ts_v_- := collectUnder(ts,v) + ts_v_+ := collectUpper(ts,v) + lgwt := internalLastSubResultant(ts_v,q,ts_v_-,false,false) + lbwt: List BWT := [] + for gwt in lgwt repeat + g := gwt.val; ts := gwt.tower + (ground? g) or (mvar(g) < v) => + ts := internalAugment(ts_v,ts) + ts := internalAugment(members(ts_v_+),ts) + lbwt := cons([true, ts]$BWT,lbwt) + g := mainPrimitivePart g + ts_g := internalAugment(g,ts) + ts_g := internalAugment(members(ts_v_+),ts_g) + -- USE internalAugment with parameters ?? + lbwt := cons([false, ts_g]$BWT,lbwt) + h := lazyPquo(ts_v,g) + (ground? h) or (mvar(h) < v) => "leave" + h := mainPrimitivePart h + ts_h := internalAugment(h,ts) + ts_h := internalAugment(members(ts_v_+),ts_h) + -- USE internalAugment with parameters ?? + -- CAN BE OPTIMIZED if the input tower is separable + inv := toseInvertible?(q,ts_h)@(List BWT) + lbwt := concat([bwt for bwt in inv | bwt.val],lbwt) + sort(#1.val < #2.val,lbwt) + + toseInvertibleSet(p:P,ts:TS): Split == + k: KeyInvSet := [p,ts] + e := extractIfCan(k)$HInvSet + e case EntryInvSet => e::EntryInvSet + q := primitivePart initiallyReduce(p,ts) + zero? q => [] + normalized?(q,ts) => [ts] + v := mvar(q) + toSave: Split := [] + not algebraic?(v,ts) => + toCheck: List BWT := toseInvertible?(init(q),ts)@(List BWT) + for bwt in toCheck repeat + bwt.val => toSave := cons(bwt.tower,toSave) + newq := removeZero(q,bwt.tower) + zero? newq => "leave" + toSave := concat(toseInvertibleSet(newq,bwt.tower), toSave) + toSave := removeDuplicates toSave + return algebraicSort(toSave)$quasicomppack + ts_v := select(ts,v)::P + ts_v_- := collectUnder(ts,v) + ts_v_+ := collectUpper(ts,v) + lgwt := internalLastSubResultant(ts_v,q,ts_v_-,false,false) + for gwt in lgwt repeat + g := gwt.val; ts := gwt.tower + (ground? g) or (mvar(g) < v) => + ts := internalAugment(ts_v,ts) + ts := internalAugment(members(ts_v_+),ts) + toSave := cons(ts,toSave) + g := mainPrimitivePart g + h := lazyPquo(ts_v,g) + h := mainPrimitivePart h + (ground? h) or (mvar(h) < v) => "leave" + ts_h := internalAugment(h,ts) + ts_h := internalAugment(members(ts_v_+),ts_h) + inv := toseInvertibleSet(q,ts_h) + toSave := removeDuplicates concat(inv,toSave) + toSave := algebraicSort(toSave)$quasicomppack + insert!(k,toSave)$HInvSet + toSave + + toseSquareFreePart_wip(p:P, ts: TS): List PWT == + -- ASSUME p is not constant and mvar(p) > mvar(ts) + -- ASSUME init(p) is invertible w.r.t. ts + -- ASSUME p is mainly primitive +-- one? mdeg(p) => [[p,ts]$PWT] + mdeg(p) = 1 => [[p,ts]$PWT] + v := mvar(p)$P + q: P := mainPrimitivePart D(p,v) + lgwt: List PWT := internalLastSubResultant(p,q,ts,true,false) + lpwt : List PWT := [] + sfp : P + for gwt in lgwt repeat + g := gwt.val; us := gwt.tower + (ground? g) or (mvar(g) < v) => + lpwt := cons([p,us],lpwt) + g := mainPrimitivePart g + sfp := lazyPquo(p,g) + sfp := mainPrimitivePart stronglyReduce(sfp,us) + lpwt := cons([sfp,us],lpwt) + lpwt + + toseSquareFreePart_base(p:P, ts: TS): List PWT == [[p,ts]$PWT] + + toseSquareFreePart(p:P, ts: TS): List PWT == toseSquareFreePart_wip(p,ts) + + prepareSubResAlgo(p1:P,p2:P,ts:TS): List LpWT == + -- ASSUME mvar(p1) = mvar(p2) > mvar(ts) and mdeg(p1) >= mdeg(p2) + -- ASSUME init(p1) invertible modulo ts !!! + toSee: List LpWT := [[[p1,p2],ts]$LpWT] + toSave: List LpWT := [] + v := mvar(p1) + while (not empty? toSee) repeat + lpwt := first toSee; toSee := rest toSee + p1 := lpwt.val.1; p2 := lpwt.val.2 + ts := lpwt.tower + lbwt := toseInvertible?(leadingCoefficient(p2,v),ts)@(List BWT) + for bwt in lbwt repeat + (bwt.val = true) and (degree(p2,v) > 0) => + p3 := prem(p1, -p2) + s: P := init(p2)**(mdeg(p1) - mdeg(p2))::N + toSave := cons([[p2,p3,s],bwt.tower]$LpWT,toSave) + -- p2 := initiallyReduce(p2,bwt.tower) + newp2 := primitivePart initiallyReduce(p2,bwt.tower) + (bwt.val = true) => + -- toSave := cons([[p2,0,1],bwt.tower]$LpWT,toSave) + toSave := cons([[p2,0,1],bwt.tower]$LpWT,toSave) + -- zero? p2 => + zero? newp2 => + toSave := cons([[p1,0,1],bwt.tower]$LpWT,toSave) + -- toSee := cons([[p1,p2],ts]$LpWT,toSee) + toSee := cons([[p1,newp2],bwt.tower]$LpWT,toSee) + toSave + + integralLastSubResultant(p1:P,p2:P,ts:TS): List PWT == + -- ASSUME mvar(p1) = mvar(p2) > mvar(ts) and mdeg(p1) >= mdeg(p2) + -- ASSUME p1 and p2 have no algebraic coefficients + lsr := lastSubResultant(p1, p2) + ground?(lsr) => [[lsr,ts]$PWT] + mvar(lsr) < mvar(p1) => [[lsr,ts]$PWT] + gi1i2 := gcd(init(p1),init(p2)) + ex: Union(P,"failed") := (gi1i2 * lsr) exquo$P init(lsr) + ex case "failed" => [[lsr,ts]$PWT] + [[ex::P,ts]$PWT] + + internalLastSubResultant(p1:P,p2:P,ts:TS,b1:B,b2:B): List PWT == + -- ASSUME mvar(p1) = mvar(p2) > mvar(ts) and mdeg(p1) >= mdeg(p2) + -- if b1 ASSUME init(p2) invertible w.r.t. ts + -- if b2 BREAK with the first non-trivial gcd + k: KeyGcd := [p1,p2,ts,b2] + e := extractIfCan(k)$HGcd + e case EntryGcd => e::EntryGcd + toSave: List PWT + empty? ts => + toSave := integralLastSubResultant(p1,p2,ts) + insert!(k,toSave)$HGcd + return toSave + toSee: List LpWT + if b1 + then + p3 := prem(p1, -p2) + s: P := init(p2)**(mdeg(p1) - mdeg(p2))::N + toSee := [[[p2,p3,s],ts]$LpWT] + else + toSee := prepareSubResAlgo(p1,p2,ts) + toSave := internalLastSubResultant(toSee,mvar(p1),b2) + insert!(k,toSave)$HGcd + toSave + + internalLastSubResultant(llpwt: List LpWT,v:V,b2:B): List PWT == + toReturn: List PWT := []; toSee: List LpWT; + while (not empty? llpwt) repeat + toSee := llpwt; llpwt := [] + -- CONSIDER FIRST the vanishing current last subresultant + for lpwt in toSee repeat + p1 := lpwt.val.1; p2 := lpwt.val.2; s := lpwt.val.3; ts := lpwt.tower + lbwt := toseInvertible?(leadingCoefficient(p2,v),ts)@(List BWT) + for bwt in lbwt repeat + bwt.val = false => + toReturn := cons([p1,bwt.tower]$PWT, toReturn) + b2 and positive?(degree(p1,v)) => return toReturn + llpwt := cons([[p1,p2,s],bwt.tower]$LpWT, llpwt) + empty? llpwt => "leave" + -- CONSIDER NOW the branches where the computations continue + toSee := llpwt; llpwt := [] + lpwt := first toSee; toSee := rest toSee + p1 := lpwt.val.1; p2 := lpwt.val.2; s := lpwt.val.3 + delta: N := (mdeg(p1) - degree(p2,v))::N + p3: P := LazardQuotient2(p2, leadingCoefficient(p2,v), s, delta) + zero?(degree(p3,v)) => + toReturn := cons([p3,lpwt.tower]$PWT, toReturn) + for lpwt in toSee repeat + toReturn := cons([p3,lpwt.tower]$PWT, toReturn) + (p1, p2) := (p3, next_subResultant2(p1, p2, p3, s)) + s := leadingCoefficient(p1,v) + llpwt := cons([[p1,p2,s],lpwt.tower]$LpWT, llpwt) + for lpwt in toSee repeat + llpwt := cons([[p1,p2,s],lpwt.tower]$LpWT, llpwt) + toReturn + + toseLastSubResultant(p1:P,p2:P,ts:TS): List PWT == + ground? p1 => + error"in toseLastSubResultantElseSplit$TOSEGCD : bad #1" + ground? p2 => + error"in toseLastSubResultantElseSplit$TOSEGCD : bad #2" + not (mvar(p2) = mvar(p1)) => + error"in toseLastSubResultantElseSplit$TOSEGCD : bad #2" + algebraic?(mvar(p1),ts) => + error"in toseLastSubResultantElseSplit$TOSEGCD : bad #1" + not initiallyReduced?(p1,ts) => + error"in toseLastSubResultantElseSplit$TOSEGCD : bad #1" + not initiallyReduced?(p2,ts) => + error"in toseLastSubResultantElseSplit$TOSEGCD : bad #2" + purelyTranscendental?(p1,ts) and purelyTranscendental?(p2,ts) => + integralLastSubResultant(p1,p2,ts) + if mdeg(p1) < mdeg(p2) then + (p1, p2) := (p2, p1) + if odd?(mdeg(p1)) and odd?(mdeg(p2)) then p2 := - p2 + internalLastSubResultant(p1,p2,ts,false,false) + +@ +\section{package RSDCMPK RegularSetDecompositionPackage} +<<package RSDCMPK RegularSetDecompositionPackage>>= +)abbrev package RSDCMPK RegularSetDecompositionPackage +++ Author: Marc Moreno Maza +++ Date Created: 09/16/1998 +++ Date Last Updated: 12/16/1998 +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ Description: +++ A package providing a new algorithm for solving polynomial systems +++ by means of regular chains. Two ways of solving are proposed: +++ in the sense of Zariski closure (like in Kalkbrener's algorithm) +++ or in the sense of the regular zeros (like in Wu, Wang or Lazard +++ methods). This algorithm is valid for nay type +++ of regular set. It does not care about the way a polynomial is +++ added in an regular set, or how two quasi-components are compared +++ (by an inclusion-test), or how the invertibility test is made in +++ the tower of simple extensions associated with a regular set. +++ These operations are realized respectively by the domain \spad{TS} +++ and the packages \axiomType{QCMPACK}(R,E,V,P,TS) and \axiomType{RSETGCD}(R,E,V,P,TS). +++ The same way it does not care about the way univariate polynomial +++ gcd (with coefficients in the tower of simple extensions associated +++ with a regular set) are computed. The only requirement is that these +++ gcd need to have invertible initials (normalized or not). +++ WARNING. There is no need for a user to call diectly any operation +++ of this package since they can be accessed by the domain \axiom{TS}. +++ Thus, the operations of this package are not documented.\newline +++ References : +++ [1] M. MORENO MAZA "A new algorithm for computing triangular +++ decomposition of algebraic varieties" NAG Tech. Rep. 4/98. +++ Version: 5. Same as 4 but Does NOT use any unproved criteria. + +RegularSetDecompositionPackage(R,E,V,P,TS): Exports == Implementation where + + R : GcdDomain + E : OrderedAbelianMonoidSup + V : OrderedSet + P : RecursivePolynomialCategory(R,E,V) + TS : RegularTriangularSetCategory(R,E,V,P) + N ==> NonNegativeInteger + Z ==> Integer + B ==> Boolean + LP ==> List P + PS ==> GeneralPolynomialSet(R,E,V,P) + PWT ==> Record(val : P, tower : TS) + BWT ==> Record(val : Boolean, tower : TS) + LpWT ==> Record(val : (List P), tower : TS) + Wip ==> Record(done: Split, todo: List LpWT) + Branch ==> Record(eq: List P, tower: TS, ineq: List P) + UBF ==> Union(Branch,"failed") + Split ==> List TS + iprintpack ==> InternalPrintPackage() + polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P) + quasicomppack ==> QuasiComponentPackage(R,E,V,P,TS) + regsetgcdpack ==> RegularTriangularSetGcdPackage(R,E,V,P,TS) + + Exports == with + + KrullNumber: (LP, Split) -> N + numberOfVariables: (LP, Split) -> N + algebraicDecompose: (P,TS,B) -> Record(done: Split, todo: List LpWT) + transcendentalDecompose: (P,TS,N) -> Record(done: Split, todo: List LpWT) + transcendentalDecompose: (P,TS) -> Record(done: Split, todo: List LpWT) + internalDecompose: (P,TS,N,B) -> Record(done: Split, todo: List LpWT) + internalDecompose: (P,TS,N) -> Record(done: Split, todo: List LpWT) + internalDecompose: (P,TS) -> Record(done: Split, todo: List LpWT) + decompose: (LP, Split, B, B) -> Split + decompose: (LP, Split, B, B, B, B, B) -> Split + upDateBranches: (LP,Split,List LpWT,Wip,N) -> List LpWT + convert: Record(val: List P,tower: TS) -> String + printInfo: (List Record(val: List P,tower: TS), N) -> Void + + Implementation == add + + KrullNumber(lp: LP, lts: Split): N == + ln: List N := [#(ts) for ts in lts] + n := #lp + reduce(max,ln) + + numberOfVariables(lp: LP, lts: Split): N == + lv: List V := variables([lp]$PS) + for ts in lts repeat lv := concat(variables(ts), lv) + # removeDuplicates(lv) + + algebraicDecompose(p: P, ts: TS, clos?: B): Record(done: Split, todo: List LpWT) == + ground? p => + error " in algebraicDecompose$REGSET: should never happen !" + v := mvar(p); n := #ts + ts_v_- := collectUnder(ts,v) + ts_v_+ := collectUpper(ts,v) + ts_v := select(ts,v)::P + if mdeg(p) < mdeg(ts_v) + then + lgwt := internalLastSubResultant(ts_v,p,ts_v_-,true,false)$regsetgcdpack + else + lgwt := internalLastSubResultant(p,ts_v,ts_v_-,true,false)$regsetgcdpack + lts: Split := [] + llpwt: List LpWT := [] + for gwt in lgwt repeat + g := gwt.val; us := gwt.tower + zero? g => + error " in algebraicDecompose$REGSET: should never happen !!" + ground? g => "leave" + if mvar(g) = v then lts := concat(augment(members(ts_v_+),augment(g,us)),lts) + h := leadingCoefficient(g,v) + b: Boolean := purelyAlgebraic?(us) + lsfp := squareFreeFactors(h)$polsetpack + lus := augment(members(ts_v_+),augment(ts_v,us)@Split) + for f in lsfp repeat + ground? f => "leave" + b and purelyAlgebraic?(f,us) => "leave" + for vs in lus repeat + llpwt := cons([[f,p],vs]$LpWT, llpwt) + [lts,llpwt] + + transcendentalDecompose(p: P, ts: TS,bound: N): Record(done: Split, todo: List LpWT) == + lts: Split + if #ts < bound + then + lts := augment(p,ts) + else + lts := [] + llpwt: List LpWT := [] + [lts,llpwt] + + transcendentalDecompose(p: P, ts: TS): Record(done: Split, todo: List LpWT) == + lts: Split:= augment(p,ts) + llpwt: List LpWT := [] + [lts,llpwt] + + internalDecompose(p: P, ts: TS,bound: N,clos?:B): Record(done: Split, todo: List LpWT) == + clos? => internalDecompose(p,ts,bound) + internalDecompose(p,ts) + + internalDecompose(p: P, ts: TS,bound: N): Record(done: Split, todo: List LpWT) == + -- ASSUME p not constant + llpwt: List LpWT := [] + lts: Split := [] + -- EITHER mvar(p) is null + if (not zero? tail(p)) and (not ground? (lmp := leastMonomial(p))) + then + llpwt := cons([[mvar(p)::P],ts]$LpWT,llpwt) + p := (p exquo lmp)::P + ip := squareFreePart init(p); tp := tail p + p := mainPrimitivePart p + -- OR init(p) is null or not + lbwt := invertible?(ip,ts)@(List BWT) + for bwt in lbwt repeat + bwt.val => + if algebraic?(mvar(p),bwt.tower) + then + rsl := algebraicDecompose(p,bwt.tower,true) + else + rsl := transcendentalDecompose(p,bwt.tower,bound) + lts := concat(rsl.done,lts) + llpwt := concat(rsl.todo,llpwt) + -- purelyAlgebraicLeadingMonomial?(ip,bwt.tower) => "leave" -- UNPROVED CRITERIA + purelyAlgebraic?(ip,bwt.tower) and purelyAlgebraic?(bwt.tower) => "leave" -- SAFE + (not ground? ip) => + zero? tp => llpwt := cons([[ip],bwt.tower]$LpWT, llpwt) + (not ground? tp) => llpwt := cons([[ip,tp],bwt.tower]$LpWT, llpwt) + riv := removeZero(ip,bwt.tower) + (zero? riv) => + zero? tp => lts := cons(bwt.tower,lts) + (not ground? tp) => llpwt := cons([[tp],bwt.tower]$LpWT, llpwt) + llpwt := cons([[riv * mainMonomial(p) + tp],bwt.tower]$LpWT, llpwt) + [lts,llpwt] + + internalDecompose(p: P, ts: TS): Record(done: Split, todo: List LpWT) == + -- ASSUME p not constant + llpwt: List LpWT := [] + lts: Split := [] + -- EITHER mvar(p) is null + if (not zero? tail(p)) and (not ground? (lmp := leastMonomial(p))) + then + llpwt := cons([[mvar(p)::P],ts]$LpWT,llpwt) + p := (p exquo lmp)::P + ip := squareFreePart init(p); tp := tail p + p := mainPrimitivePart p + -- OR init(p) is null or not + lbwt := invertible?(ip,ts)@(List BWT) + for bwt in lbwt repeat + bwt.val => + if algebraic?(mvar(p),bwt.tower) + then + rsl := algebraicDecompose(p,bwt.tower,false) + else + rsl := transcendentalDecompose(p,bwt.tower) + lts := concat(rsl.done,lts) + llpwt := concat(rsl.todo,llpwt) + purelyAlgebraic?(ip,bwt.tower) and purelyAlgebraic?(bwt.tower) => "leave" + (not ground? ip) => + zero? tp => llpwt := cons([[ip],bwt.tower]$LpWT, llpwt) + (not ground? tp) => llpwt := cons([[ip,tp],bwt.tower]$LpWT, llpwt) + riv := removeZero(ip,bwt.tower) + (zero? riv) => + zero? tp => lts := cons(bwt.tower,lts) + (not ground? tp) => llpwt := cons([[tp],bwt.tower]$LpWT, llpwt) + llpwt := cons([[riv * mainMonomial(p) + tp],bwt.tower]$LpWT, llpwt) + [lts,llpwt] + + decompose(lp: LP, lts: Split, clos?: B, info?: B): Split == + decompose(lp,lts,false,false,clos?,true,info?) + + convert(lpwt: LpWT): String == + ls: List String := ["<", string((#(lpwt.val))::Z), ",", string((#(lpwt.tower))::Z), ">" ] + concat ls + + printInfo(toSee: List LpWT, n: N): Void == + lpwt := first toSee + s: String := concat ["[", string((#toSee)::Z), " ", convert(lpwt)@String] + m: N := #(lpwt.val) + toSee := rest toSee + for lpwt in toSee repeat + m := m + #(lpwt.val) + s := concat [s, ",", convert(lpwt)@String] + s := concat [s, " -> |", string(m::Z), "|; {", string(n::Z),"}]"] + iprint(s)$iprintpack + void() + + decompose(lp: LP, lts: Split, cleanW?: B, sqfr?: B, clos?: B, rem?: B, info?: B): Split == + -- if cleanW? then REMOVE REDUNDANT COMPONENTS in lts + -- if sqfr? then SPLIT the system with SQUARE-FREE FACTORIZATION + -- if clos? then SOLVE in the closure sense + -- if rem? then REDUCE the current p by using remainder + -- if info? then PRINT info + empty? lp => lts + branches: List Branch := prepareDecompose(lp,lts,cleanW?,sqfr?)$quasicomppack + empty? branches => [] + toSee: List LpWT := [[br.eq,br.tower]$LpWT for br in branches] + toSave: Split := [] + if clos? then bound := KrullNumber(lp,lts) else bound := numberOfVariables(lp,lts) + while (not empty? toSee) repeat + if info? then printInfo(toSee,#toSave) + lpwt := first toSee; toSee := rest toSee + lp := lpwt.val; ts := lpwt.tower + empty? lp => + toSave := cons(ts, toSave) + p := first lp; lp := rest lp + if rem? and (not ground? p) and (not empty? ts) + then + p := remainder(p,ts).polnum + p := removeZero(p,ts) + zero? p => toSee := cons([lp,ts]$LpWT, toSee) + ground? p => "leave" + rsl := internalDecompose(p,ts,bound,clos?) + toSee := upDateBranches(lp,toSave,toSee,rsl,bound) + removeSuperfluousQuasiComponents(toSave)$quasicomppack + + upDateBranches(leq:LP,lts:Split,current:List LpWT,wip: Wip,n:N): List LpWT == + newBranches: List LpWT := wip.todo + newComponents: Split := wip.done + branches1, branches2: List LpWT + branches1 := []; branches2 := [] + for branch in newBranches repeat + us := branch.tower + #us > n => "leave" + newleq := sort(infRittWu?,concat(leq,branch.val)) + --foo := rewriteSetWithReduction(newleq,us,initiallyReduce,initiallyReduced?) + --any?(ground?,foo) => "leave" + branches1 := cons([newleq,us]$LpWT, branches1) + for us in newComponents repeat + #us > n => "leave" + subQuasiComponent?(us,lts)$quasicomppack => "leave" + --newleq := leq + --foo := rewriteSetWithReduction(newleq,us,initiallyReduce,initiallyReduced?) + --any?(ground?,foo) => "leave" + branches2 := cons([leq,us]$LpWT, branches2) + empty? branches1 => + empty? branches2 => current + concat(branches2, current) + branches := concat [branches2, branches1, current] + -- branches := concat(branches,current) + removeSuperfluousCases(branches)$quasicomppack + +@ +\section{domain REGSET RegularTriangularSet} +Several domain constructors implement regular triangular sets (or regular +chains). Among them {\bf RegularTriangularSet} and +{\bf SquareFreeRegularTriangularSet}. They also implement an algorithm +by Marc Moreno Maza for computing triangular decompositions of polynomial +systems. This method is refined in the package {\bf LazardSetSolvingPackage} +in order to produce decompositions by means of Lazard triangular sets. +<<domain REGSET RegularTriangularSet>>= +)abbrev domain REGSET RegularTriangularSet +++ Author: Marc Moreno Maza +++ Date Created: 08/25/1998 +++ Date Last Updated: 16/12/1998 +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ Description: +++ This domain provides an implementation of regular chains. +++ Moreover, the operation \axiomOpFrom{zeroSetSplit}{RegularTriangularSetCategory} +++ is an implementation of a new algorithm for solving polynomial systems by +++ means of regular chains.\newline +++ References : +++ [1] M. MORENO MAZA "A new algorithm for computing triangular +++ decomposition of algebraic varieties" NAG Tech. Rep. 4/98. +++ Version: Version 11. + +RegularTriangularSet(R,E,V,P) : Exports == Implementation where + + R : GcdDomain + E : OrderedAbelianMonoidSup + V : OrderedSet + P : RecursivePolynomialCategory(R,E,V) + N ==> NonNegativeInteger + Z ==> Integer + B ==> Boolean + LP ==> List P + PtoP ==> P -> P + PS ==> GeneralPolynomialSet(R,E,V,P) + PWT ==> Record(val : P, tower : $) + BWT ==> Record(val : Boolean, tower : $) + LpWT ==> Record(val : (List P), tower : $) + Split ==> List $ + iprintpack ==> InternalPrintPackage() + polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P) + quasicomppack ==> QuasiComponentPackage(R,E,V,P,$) + regsetgcdpack ==> RegularTriangularSetGcdPackage(R,E,V,P,$) + regsetdecomppack ==> RegularSetDecompositionPackage(R,E,V,P,$) + + Exports == RegularTriangularSetCategory(R,E,V,P) with + + internalAugment: (P,$,B,B,B,B,B) -> List $ + ++ \axiom{internalAugment(p,ts,b1,b2,b3,b4,b5)} + ++ is an internal subroutine, exported only for developement. + zeroSetSplit: (LP, B, B) -> Split + ++ \axiom{zeroSetSplit(lp,clos?,info?)} has the same specifications as + ++ \axiomOpFrom{zeroSetSplit}{RegularTriangularSetCategory}. + ++ Moreover, if \axiom{clos?} then solves in the sense of the Zariski closure + ++ else solves in the sense of the regular zeros. If \axiom{info?} then + ++ do print messages during the computations. + zeroSetSplit: (LP, B, B, B, B) -> Split + ++ \axiom{zeroSetSplit(lp,b1,b2.b3,b4)} + ++ is an internal subroutine, exported only for developement. + internalZeroSetSplit: (LP, B, B, B) -> Split + ++ \axiom{internalZeroSetSplit(lp,b1,b2,b3)} + ++ is an internal subroutine, exported only for developement. + pre_process: (LP, B, B) -> Record(val: LP, towers: Split) + ++ \axiom{pre_process(lp,b1,b2)} + ++ is an internal subroutine, exported only for developement. + + Implementation == add + + Rep ==> LP + + rep(s:$):Rep == s pretend Rep + per(l:Rep):$ == l pretend $ + + copy ts == + per(copy(rep(ts))$LP) + empty() == + per([]) + empty?(ts:$) == + empty?(rep(ts)) + parts ts == + rep(ts) + members ts == + rep(ts) + map (f : PtoP, ts : $) : $ == + construct(map(f,rep(ts))$LP)$$ + map! (f : PtoP, ts : $) : $ == + construct(map!(f,rep(ts))$LP)$$ + member? (p,ts) == + member?(p,rep(ts))$LP + unitIdealIfCan() == + "failed"::Union($,"failed") + roughUnitIdeal? ts == + false + coerce(ts:$) : OutputForm == + lp : List(P) := reverse(rep(ts)) + brace([p::OutputForm for p in lp]$List(OutputForm))$OutputForm + mvar ts == + empty? ts => error "mvar$REGSET: #1 is empty" + mvar(first(rep(ts)))$P + first ts == + empty? ts => "failed"::Union(P,"failed") + first(rep(ts))::Union(P,"failed") + last ts == + empty? ts => "failed"::Union(P,"failed") + last(rep(ts))::Union(P,"failed") + rest ts == + empty? ts => "failed"::Union($,"failed") + per(rest(rep(ts)))::Union($,"failed") + coerce(ts:$) : (List P) == + rep(ts) + + collectUpper (ts,v) == + empty? ts => ts + lp := rep(ts) + newlp : Rep := [] + while (not empty? lp) and (mvar(first(lp)) > v) repeat + newlp := cons(first(lp),newlp) + lp := rest lp + per(reverse(newlp)) + + collectUnder (ts,v) == + empty? ts => ts + lp := rep(ts) + while (not empty? lp) and (mvar(first(lp)) >= v) repeat + lp := rest lp + per(lp) + + construct(lp:List(P)) == + ts : $ := per([]) + empty? lp => ts + lp := sort(infRittWu?,lp) + while not empty? lp repeat + eif := extendIfCan(ts,first(lp)) + not (eif case $) => + error"in construct : List P -> $ from REGSET : bad #1" + ts := eif::$ + lp := rest lp + ts + + extendIfCan(ts:$,p:P) == + ground? p => "failed"::Union($,"failed") + empty? ts => + p := primitivePart p + (per([p]))::Union($,"failed") + not (mvar(ts) < mvar(p)) => "failed"::Union($,"failed") + invertible?(init(p),ts)@Boolean => + (per(cons(p,rep(ts))))::Union($,"failed") + "failed"::Union($,"failed") + + removeZero(p:P, ts:$): P == + (ground? p) or (empty? ts) => p + v := mvar(p) + ts_v_- := collectUnder(ts,v) + if algebraic?(v,ts) + then + q := lazyPrem(p,select(ts,v)::P) + zero? q => return q + zero? removeZero(q,ts_v_-) => return 0 + empty? ts_v_- => p + q: P := 0 + while positive? degree(p,v) repeat + q := removeZero(init(p),ts_v_-) * mainMonomial(p) + q + p := tail(p) + q + removeZero(p,ts_v_-) + + internalAugment(p:P,ts:$): $ == + -- ASSUME that adding p to ts DOES NOT require any split + ground? p => error "in internalAugment$REGSET: ground? #1" + first(internalAugment(p,ts,false,false,false,false,false)) + + internalAugment(lp:List(P),ts:$): $ == + -- ASSUME that adding p to ts DOES NOT require any split + empty? lp => ts + internalAugment(rest lp, internalAugment(first lp, ts)) + + internalAugment(p:P,ts:$,rem?:B,red?:B,prim?:B,sqfr?:B,extend?:B): Split == + -- ASSUME p is not a constant + -- ASSUME mvar(p) is not algebraic w.r.t. ts + -- ASSUME init(p) invertible modulo ts + -- if rem? then REDUCE p by remainder + -- if prim? then REPLACE p by its main primitive part + -- if sqfr? then FACTORIZE SQUARE FREE p over R + -- if extend? DO NOT ASSUME every pol in ts_v_+ is invertible modulo ts + v := mvar(p) + ts_v_- := collectUnder(ts,v) + ts_v_+ := collectUpper(ts,v) + if rem? then p := remainder(p,ts_v_-).polnum + -- if rem? then p := reduceByQuasiMonic(p,ts_v_-) + if red? then p := removeZero(p,ts_v_-) + if prim? then p := mainPrimitivePart p + if sqfr? + then + lsfp := squareFreeFactors(p)$polsetpack + lts: Split := [per(cons(f,rep(ts_v_-))) for f in lsfp] + else + lts: Split := [per(cons(p,rep(ts_v_-)))] + extend? => extend(members(ts_v_+),lts) + [per(concat(rep(ts_v_+),rep(us))) for us in lts] + + augment(p:P,ts:$): List $ == + ground? p => error "in augment$REGSET: ground? #1" + algebraic?(mvar(p),ts) => error "in augment$REGSET: bad #1" + -- ASSUME init(p) invertible modulo ts + -- DOES NOT ASSUME anything else. + -- THUS reduction, mainPrimitivePart and squareFree are NEEDED + internalAugment(p,ts,true,true,true,true,true) + + extend(p:P,ts:$): List $ == + ground? p => error "in extend$REGSET: ground? #1" + v := mvar(p) + not (mvar(ts) < mvar(p)) => error "in extend$REGSET: bad #1" + lts: List($) := [] + split: List($) := invertibleSet(init(p),ts) + for us in split repeat + lts := concat(augment(p,us),lts) + lts + + invertible?(p:P,ts:$): Boolean == + toseInvertible?(p,ts)$regsetgcdpack + + invertible?(p:P,ts:$): List BWT == + toseInvertible?(p,ts)$regsetgcdpack + + invertibleSet(p:P,ts:$): Split == + toseInvertibleSet(p,ts)$regsetgcdpack + + lastSubResultant(p1:P,p2:P,ts:$): List PWT == + toseLastSubResultant(p1,p2,ts)$regsetgcdpack + + squareFreePart(p:P, ts: $): List PWT == + toseSquareFreePart(p,ts)$regsetgcdpack + + intersect(p:P, ts: $): List($) == decompose([p], [ts], false, false)$regsetdecomppack + + intersect(lp: LP, lts: List($)): List($) == decompose(lp, lts, false, false)$regsetdecomppack + -- SOLVE in the regular zero sense + -- and DO NOT PRINT info + + decompose(p:P, ts: $): List($) == decompose([p], [ts], true, false)$regsetdecomppack + + decompose(lp: LP, lts: List($)): List($) == decompose(lp, lts, true, false)$regsetdecomppack + -- SOLVE in the closure sense + -- and DO NOT PRINT info + + zeroSetSplit(lp:List(P)) == zeroSetSplit(lp,true,false) + -- by default SOLVE in the closure sense + -- and DO NOT PRINT info + + zeroSetSplit(lp:List(P), clos?: B) == zeroSetSplit(lp,clos?, false) + -- DO NOT PRINT info + + zeroSetSplit(lp:List(P), clos?: B, info?: B) == + -- if clos? then SOLVE in the closure sense + -- if info? then PRINT info + -- by default USE hash-tables + -- and PREPROCESS the input system + zeroSetSplit(lp,true,clos?,info?,true) + + zeroSetSplit(lp:List(P),hash?:B,clos?:B,info?:B,prep?:B) == + -- if hash? then USE hash-tables + -- if info? then PRINT information + -- if clos? then SOLVE in the closure sense + -- if prep? then PREPROCESS the input system + if hash? + then + s1, s2, s3, dom1, dom2, dom3: String + e: String := empty()$String + if info? then (s1,s2,s3) := ("w","g","i") else (s1,s2,s3) := (e,e,e) + if info? + then + (dom1, dom2, dom3) := ("QCMPACK", "REGSETGCD: Gcd", "REGSETGCD: Inv Set") + else + (dom1, dom2, dom3) := (e,e,e) + startTable!(s1,"W",dom1)$quasicomppack + startTableGcd!(s2,"G",dom2)$regsetgcdpack + startTableInvSet!(s3,"I",dom3)$regsetgcdpack + lts := internalZeroSetSplit(lp,clos?,info?,prep?) + if hash? + then + stopTable!()$quasicomppack + stopTableGcd!()$regsetgcdpack + stopTableInvSet!()$regsetgcdpack + lts + + internalZeroSetSplit(lp:LP,clos?:B,info?:B,prep?:B) == + -- if info? then PRINT information + -- if clos? then SOLVE in the closure sense + -- if prep? then PREPROCESS the input system + if prep? + then + pp := pre_process(lp,clos?,info?) + lp := pp.val + lts := pp.towers + else + ts: $ := [[]] + lts := [ts] + lp := remove(zero?, lp) + any?(ground?, lp) => [] + empty? lp => lts + empty? lts => lts + lp := sort(infRittWu?,lp) + clos? => decompose(lp,lts, clos?, info?)$regsetdecomppack + -- IN DIM > 0 with clos? the following is false ... + for p in lp repeat + lts := decompose([p],lts, clos?, info?)$regsetdecomppack + lts + + largeSystem?(lp:LP): Boolean == + -- Gonnet and Gerdt and not Wu-Wang.2 + #lp > 16 => true + #lp < 13 => false + lts: List($) := [] + (#lp :: Z - numberOfVariables(lp,lts)$regsetdecomppack :: Z) > 3 + + smallSystem?(lp:LP): Boolean == + -- neural, Vermeer, Liu, and not f-633 and not Hairer-2 + #lp < 5 + + mediumSystem?(lp:LP): Boolean == + -- f-633 and not Hairer-2 + lts: List($) := [] + (numberOfVariables(lp,lts)$regsetdecomppack :: Z - #lp :: Z) < 2 + +-- lin?(p:P):Boolean == ground?(init(p)) and one?(mdeg(p)) + lin?(p:P):Boolean == ground?(init(p)) and (mdeg(p) = 1) + + pre_process(lp:LP,clos?:B,info?:B): Record(val: LP, towers: Split) == + -- if info? then PRINT information + -- if clos? then SOLVE in the closure sense + ts: $ := [[]]; + lts: Split := [ts] + empty? lp => [lp,lts] + lp1: List P := [] + lp2: List P := [] + for p in lp repeat + ground? (tail p) => lp1 := cons(p, lp1) + lp2 := cons(p, lp2) + lts: Split := decompose(lp1,[ts],clos?,info?)$regsetdecomppack + probablyZeroDim?(lp)$polsetpack => + largeSystem?(lp) => return [lp2,lts] + if #lp > 7 + then + -- Butcher (8,8) + Wu-Wang.2 (13,16) + lp2 := crushedSet(lp2)$polsetpack + lp2 := remove(zero?,lp2) + any?(ground?,lp2) => return [lp2, lts] + lp3 := [p for p in lp2 | lin?(p)] + lp4 := [p for p in lp2 | not lin?(p)] + if clos? + then + lts := decompose(lp4,lts, clos?, info?)$regsetdecomppack + else + lp4 := sort(infRittWu?,lp4) + for p in lp4 repeat + lts := decompose([p],lts, clos?, info?)$regsetdecomppack + lp2 := lp3 + else + lp2 := crushedSet(lp2)$polsetpack + lp2 := remove(zero?,lp2) + any?(ground?,lp2) => return [lp2, lts] + if clos? + then + lts := decompose(lp2,lts, clos?, info?)$regsetdecomppack + else + lp2 := sort(infRittWu?,lp2) + for p in lp2 repeat + lts := decompose([p],lts, clos?, info?)$regsetdecomppack + lp2 := [] + return [lp2,lts] + smallSystem?(lp) => [lp2,lts] + mediumSystem?(lp) => [crushedSet(lp2)$polsetpack,lts] + lp3 := [p for p in lp2 | lin?(p)] + lp4 := [p for p in lp2 | not lin?(p)] + if clos? + then + lts := decompose(lp4,lts, clos?, info?)$regsetdecomppack + else + lp4 := sort(infRittWu?,lp4) + for p in lp4 repeat + lts := decompose([p],lts, clos?, info?)$regsetdecomppack + if clos? + then + lts := decompose(lp3,lts, clos?, info?)$regsetdecomppack + else + lp3 := sort(infRittWu?,lp3) + for p in lp3 repeat + lts := decompose([p],lts, clos?, info?)$regsetdecomppack + lp2 := [] + return [lp2,lts] + +@ +\section{License} +<<license>>= +--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd. +--All rights reserved. +-- +--Redistribution and use in source and binary forms, with or without +--modification, are permitted provided that the following conditions are +--met: +-- +-- - Redistributions of source code must retain the above copyright +-- notice, this list of conditions and the following disclaimer. +-- +-- - Redistributions in binary form must reproduce the above copyright +-- notice, this list of conditions and the following disclaimer in +-- the documentation and/or other materials provided with the +-- distribution. +-- +-- - Neither the name of The Numerical ALgorithms Group Ltd. nor the +-- names of its contributors may be used to endorse or promote products +-- derived from this software without specific prior written permission. +-- +--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS +--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED +--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A +--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER +--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, +--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, +--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR +--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF +--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING +--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +@ +<<*>>= +<<license>> + +<<category RSETCAT RegularTriangularSetCategory>> +<<package QCMPACK QuasiComponentPackage>> +<<package RSETGCD RegularTriangularSetGcdPackage>> +<<package RSDCMPK RegularSetDecompositionPackage>> +<<domain REGSET RegularTriangularSet>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |