diff options
Diffstat (limited to 'src/algebra/sregset.spad.pamphlet')
| -rw-r--r-- | src/algebra/sregset.spad.pamphlet | 1606 |
1 files changed, 1606 insertions, 0 deletions
diff --git a/src/algebra/sregset.spad.pamphlet b/src/algebra/sregset.spad.pamphlet new file mode 100644 index 00000000..19886881 --- /dev/null +++ b/src/algebra/sregset.spad.pamphlet @@ -0,0 +1,1606 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra sregset.spad} +\author{Marc Moreno Maza} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{category SFRTCAT SquareFreeRegularTriangularSetCategory} +<<category SFRTCAT SquareFreeRegularTriangularSetCategory>>= +)abbrev category SFRTCAT SquareFreeRegularTriangularSetCategory +++ Author: Marc Moreno Maza +++ Date Created: 09/03/1996 +++ Date Last Updated: 09/10/1998 +++ Basic Functions: +++ Related Constructors: +++ Also See: essai Graphisme +++ AMS Classifications: +++ Keywords: polynomial, multivariate, ordered variables set +++ Description: +++ The category of square-free regular triangular sets. +++ A regular triangular set \spad{ts} is square-free if +++ the gcd of any polynomial \spad{p} in \spad{ts} and +++ \spad{differentiate(p,mvar(p))} w.r.t. +++ \axiomOpFrom{collectUnder}{TriangularSetCategory}(ts,\axiomOpFrom{mvar}{RecursivePolynomialCategory}(p)) +++ has degree zero w.r.t. \spad{mvar(p)}. Thus any square-free regular +++ set defines a tower of square-free simple extensions.\newline +++ References : +++ [1] D. LAZARD "A new method for solving algebraic systems of +++ positive dimension" Discr. App. Math. 33:147-160,1991 +++ [2] M. KALKBRENER "Algorithmic properties of polynomial rings" +++ Habilitation Thesis, ETZH, Zurich, 1995. +++ [3] M. MORENO MAZA "A new algorithm for computing triangular +++ decomposition of algebraic varieties" NAG Tech. Rep. 4/98. + + + +SquareFreeRegularTriangularSetCategory(R:GcdDomain,E:OrderedAbelianMonoidSup,_ + V:OrderedSet,P:RecursivePolynomialCategory(R,E,V)): + Category == + RegularTriangularSetCategory(R,E,V,P) + +@ +\section{package SFQCMPK SquareFreeQuasiComponentPackage} +<<package SFQCMPK SquareFreeQuasiComponentPackage>>= +)abbrev package SFQCMPK SquareFreeQuasiComponentPackage +++ Author: Marc Moreno Maza +++ Date Created: 09/23/1998 +++ Date Last Updated: 12/16/1998 +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ Description: +++ A internal 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 +++ Tech. Report (PoSSo project) +++ [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: 1. + +SquareFreeQuasiComponentPackage(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) + SQUAREFREE ==> SquareFreeRegularTriangularSetCategory(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 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?(ts,us)}{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) + + if TS has SquareFreeRegularTriangularSetCategory(R,E,V,P) + then + + internalSubQuasiComponent?(ts:TS,us:TS): Union(Boolean,"failed") == + 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 + b: B := invertible?(p,ts)$TS + not b => + return(false::Union(Boolean,"failed")) + true::Union(Boolean,"failed") + + else + + internalSubQuasiComponent?(ts:TS,us:TS): Union(Boolean,"failed") == + 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? reduceByQuasiMonic(p,ts) => + return("failed"::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 SFRGCD SquareFreeRegularTriangularSetGcdPackage} +<<package SFRGCD SquareFreeRegularTriangularSetGcdPackage>>= +)abbrev package SFRGCD SquareFreeRegularTriangularSetGcdPackage +++ Author: Marc Moreno Maza +++ Date Created: 09/23/1998 +++ Date Last Updated: 10/01/1998 +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ Description: +++ A internal package for computing gcds and resultants of univariate polynomials +++ with coefficients in a tower of simple extensions of a field. +++ There is no need to use directly this package since its main operations are +++ available from \spad{TS}. \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: 1. + +SquareFreeRegularTriangularSetGcdPackage(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) + iprintpack ==> InternalPrintPackage() + polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P) + quasicomppack ==> SquareFreeQuasiComponentPackage(R,E,V,P,TS) + + SQUAREFREE ==> SquareFreeRegularTriangularSetCategory(R,E,V,P) + + Exports == with + startTableGcd!: (S,S,S) -> Void + stopTableGcd!: () -> Void + startTableInvSet!: (S,S,S) -> Void + stopTableInvSet!: () -> Void + stosePrepareSubResAlgo: (P,P,TS) -> List LpWT + stoseInternalLastSubResultant: (P,P,TS,B,B) -> List PWT + stoseInternalLastSubResultant: (List LpWT,V,B) -> List PWT + stoseIntegralLastSubResultant: (P,P,TS) -> List PWT + stoseLastSubResultant: (P,P,TS) -> List PWT + stoseInvertible?: (P,TS) -> B + stoseInvertible?_sqfreg: (P,TS) -> List BWT + stoseInvertibleSet_sqfreg: (P,TS) -> Split + stoseInvertible?_reg: (P,TS) -> List BWT + stoseInvertibleSet_reg: (P,TS) -> Split + stoseInvertible?: (P,TS) -> List BWT + stoseInvertibleSet: (P,TS) -> Split + stoseSquareFreePart: (P,TS) -> List PWT + + 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 + + stoseInvertible?(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 := stoseInvertible?(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 := stoseInternalLastSubResultant(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 + + stosePrepareSubResAlgo(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 := stoseInvertible?(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],bwt.tower]$LpWT,toSee) + toSee := cons([[p1,newp2],bwt.tower]$LpWT,toSee) + toSave + + stoseIntegralLastSubResultant(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] + + stoseInternalLastSubResultant(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 := stoseIntegralLastSubResultant(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 := stosePrepareSubResAlgo(p1,p2,ts) + toSave := stoseInternalLastSubResultant(toSee,mvar(p1),b2) + insert!(k,toSave)$HGcd + toSave + + stoseInternalLastSubResultant(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 := stoseInvertible?(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 + + stoseLastSubResultant(p1:P,p2:P,ts:TS): List PWT == + ground? p1 => + error"in stoseLastSubResultantElseSplit$SFRGCD : bad #1" + ground? p2 => + error"in stoseLastSubResultantElseSplit$SFRGCD : bad #2" + not (mvar(p2) = mvar(p1)) => + error"in stoseLastSubResultantElseSplit$SFRGCD : bad #2" + algebraic?(mvar(p1),ts) => + error"in stoseLastSubResultantElseSplit$SFRGCD : bad #1" + not initiallyReduced?(p1,ts) => + error"in stoseLastSubResultantElseSplit$SFRGCD : bad #1" + not initiallyReduced?(p2,ts) => + error"in stoseLastSubResultantElseSplit$SFRGCD : bad #2" + purelyTranscendental?(p1,ts) and purelyTranscendental?(p2,ts) => + stoseIntegralLastSubResultant(p1,p2,ts) + if mdeg(p1) < mdeg(p2) then + (p1, p2) := (p2, p1) + if odd?(mdeg(p1)) and odd?(mdeg(p2)) then p2 := - p2 + stoseInternalLastSubResultant(p1,p2,ts,false,false) + + stoseSquareFreePart_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 := stoseInternalLastSubResultant(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 + + stoseSquareFreePart_base(p:P, ts: TS): List PWT == [[p,ts]$PWT] + + stoseSquareFreePart(p:P, ts: TS): List PWT == stoseSquareFreePart_wip(p,ts) + + stoseInvertible?_sqfreg(p:P,ts:TS): List BWT == + --iprint("+")$iprintpack + 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 := stoseInvertible?_sqfreg(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(stoseInvertible?_sqfreg(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 := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false) + lbwt: List BWT := [] + lts, lts_g, lts_h: Split + for gwt in lgwt repeat + g := gwt.val; ts := gwt.tower + (ground? g) or (mvar(g) < v) => + lts := augment(ts_v,ts)$TS + lts := augment(members(ts_v_+),lts)$TS + for ts in lts repeat + lbwt := cons([true, ts]$BWT,lbwt) + g := mainPrimitivePart g + lts_g := augment(g,ts)$TS + lts_g := augment(members(ts_v_+),lts_g)$TS + -- USE stoseInternalAugment with parameters ?? + for ts_g in lts_g repeat + lbwt := cons([false, ts_g]$BWT,lbwt) + h := lazyPquo(ts_v,g) + (ground? h) or (mvar(h) < v) => "leave" + h := mainPrimitivePart h + lts_h := augment(h,ts)$TS + lts_h := augment(members(ts_v_+),lts_h)$TS + -- USE stoseInternalAugment with parameters ?? + for ts_h in lts_h repeat + lbwt := cons([true, ts_h]$BWT,lbwt) + sort(#1.val < #2.val,lbwt) + + stoseInvertibleSet_sqfreg(p:P,ts:TS): Split == + --iprint("*")$iprintpack + 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 := stoseInvertible?_sqfreg(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(stoseInvertibleSet_sqfreg(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 := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false) + lts, lts_h: Split + for gwt in lgwt repeat + g := gwt.val; ts := gwt.tower + (ground? g) or (mvar(g) < v) => + lts := augment(ts_v,ts)$TS + lts := augment(members(ts_v_+),lts)$TS + toSave := concat(lts,toSave) + g := mainPrimitivePart g + h := lazyPquo(ts_v,g) + h := mainPrimitivePart h + (ground? h) or (mvar(h) < v) => "leave" + lts_h := augment(h,ts)$TS + lts_h := augment(members(ts_v_+),lts_h)$TS + toSave := concat(lts_h,toSave) + toSave := algebraicSort(toSave)$quasicomppack + insert!(k,toSave)$HInvSet + toSave + + stoseInvertible?_reg(p:P,ts:TS): List BWT == + --iprint("-")$iprintpack + 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 := stoseInvertible?_reg(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(stoseInvertible?_reg(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 := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false) + lbwt: List BWT := [] + lts, lts_g, lts_h: Split + for gwt in lgwt repeat + g := gwt.val; ts := gwt.tower + (ground? g) or (mvar(g) < v) => + lts := augment(ts_v,ts)$TS + lts := augment(members(ts_v_+),lts)$TS + for ts in lts repeat + lbwt := cons([true, ts]$BWT,lbwt) + g := mainPrimitivePart g + lts_g := augment(g,ts)$TS + lts_g := augment(members(ts_v_+),lts_g)$TS + -- USE internalAugment with parameters ?? + for ts_g in lts_g repeat + lbwt := cons([false, ts_g]$BWT,lbwt) + h := lazyPquo(ts_v,g) + (ground? h) or (mvar(h) < v) => "leave" + h := mainPrimitivePart h + lts_h := augment(h,ts)$TS + lts_h := augment(members(ts_v_+),lts_h)$TS + -- USE internalAugment with parameters ?? + for ts_h in lts_h repeat + inv := stoseInvertible?_reg(q,ts_h)@(List BWT) + lbwt := concat([bwt for bwt in inv | bwt.val],lbwt) + sort(#1.val < #2.val,lbwt) + + stoseInvertibleSet_reg(p:P,ts:TS): Split == + --iprint("/")$iprintpack + 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 := stoseInvertible?_reg(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(stoseInvertibleSet_reg(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 := stoseInternalLastSubResultant(ts_v,q,ts_v_-,false,false) + lts, lts_h: Split + for gwt in lgwt repeat + g := gwt.val; ts := gwt.tower + (ground? g) or (mvar(g) < v) => + lts := augment(ts_v,ts)$TS + lts := augment(members(ts_v_+),lts)$TS + toSave := concat(lts,toSave) + g := mainPrimitivePart g + h := lazyPquo(ts_v,g) + h := mainPrimitivePart h + (ground? h) or (mvar(h) < v) => "leave" + lts_h := augment(h,ts)$TS + lts_h := augment(members(ts_v_+),lts_h)$TS + for ts_h in lts_h repeat + inv := stoseInvertibleSet_reg(q,ts_h) + toSave := removeDuplicates concat(inv,toSave) + toSave := algebraicSort(toSave)$quasicomppack + insert!(k,toSave)$HInvSet + toSave + + if TS has SquareFreeRegularTriangularSetCategory(R,E,V,P) + then + + stoseInvertible?(p:P,ts:TS): List BWT == stoseInvertible?_sqfreg(p,ts) + + stoseInvertibleSet(p:P,ts:TS): Split == stoseInvertibleSet_sqfreg(p,ts) + + else + + stoseInvertible?(p:P,ts:TS): List BWT == stoseInvertible?_reg(p,ts) + + stoseInvertibleSet(p:P,ts:TS): Split == stoseInvertibleSet_reg(p,ts) + +@ +\section{package SRDCMPK SquareFreeRegularSetDecompositionPackage} +<<package SRDCMPK SquareFreeRegularSetDecompositionPackage>>= +)abbrev package SRDCMPK SquareFreeRegularSetDecompositionPackage +++ Author: Marc Moreno Maza +++ Date Created: 09/23/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 provided: +++ 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- +++ Moreno 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 \spad{QCMPPK(R,E,V,P,TS)} and \spad{RSETGCD(R,E,V,P,TS)}. +++ The same way it does not care about the way univariate polynomial +++ gcds (with coefficients in the tower of simple extensions associated +++ with a regular set) are computed. The only requirement is that these +++ gcds 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 \axiomType{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: 2. Does not use any unproved criteria. + +SquareFreeRegularSetDecompositionPackage(R,E,V,P,TS): Exports == Implementation where + + R : GcdDomain + E : OrderedAbelianMonoidSup + V : OrderedSet + P : RecursivePolynomialCategory(R,E,V) + TS : SquareFreeRegularTriangularSetCategory(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 ==> SquareFreeQuasiComponentPackage(R,E,V,P,TS) + regsetgcdpack ==> SquareFreeRegularTriangularSetGcdPackage(R,E,V,P,TS) + + Exports == with + + KrullNumber: (LP, Split) -> N + numberOfVariables: (LP, Split) -> N + algebraicDecompose: (P,TS) -> 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): 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 + lgwt: List PWT + if mdeg(p) < mdeg(ts_v) + then + lgwt := stoseInternalLastSubResultant(ts_v,p,ts_v_-,true,false)$regsetgcdpack + else + lgwt := stoseInternalLastSubResultant(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" + h := leadingCoefficient(g,v) + lus := augment(members(ts_v_+),augment(ts_v,us)$TS)$TS + lsfp := squareFreeFactors(h)$polsetpack + for f in lsfp repeat + ground? f => "leave" + for vs in lus repeat + llpwt := cons([[f,p],vs]$LpWT, llpwt) + n < #us => + error " in algebraicDecompose$REGSET: should never happen !!!" + mvar(g) = v => + lts := concat(augment(members(ts_v_+),augment(g,us)$TS)$TS,lts) + [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)$TS + else + lts := [] + llpwt: List LpWT := [] + [lts,llpwt] + + transcendentalDecompose(p: P, ts: TS): Record(done: Split, todo: List LpWT) == + lts: Split:= augment(p,ts)$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: List BWT := stoseInvertible?_sqfreg(ip,ts)$regsetgcdpack + for bwt in lbwt repeat + bwt.val => + if algebraic?(mvar(p),bwt.tower) + then + rsl := algebraicDecompose(p,bwt.tower) + else + rsl := transcendentalDecompose(p,bwt.tower,bound) + lts := concat(rsl.done,lts) + llpwt := concat(rsl.todo,llpwt) + (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: List BWT := stoseInvertible?_sqfreg(ip,ts)$regsetgcdpack + for bwt in lbwt repeat + bwt.val => + if algebraic?(mvar(p),bwt.tower) + then + rsl := algebraicDecompose(p,bwt.tower) + else + rsl := transcendentalDecompose(p,bwt.tower) + lts := concat(rsl.done,lts) + llpwt := concat(rsl.todo,llpwt) + (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 SREGSET SquareFreeRegularTriangularSet} +<<domain SREGSET SquareFreeRegularTriangularSet>>= +)abbrev domain SREGSET SquareFreeRegularTriangularSet +++ 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 square-free regular chains. +++ Moreover, the operation \axiomOpFrom{zeroSetSplit}{SquareFreeRegularTriangularSetCategory} +++ 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: 2 + +SquareFreeRegularTriangularSet(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 ==> SquareFreeQuasiComponentPackage(R,E,V,P,$) + regsetgcdpack ==> SquareFreeRegularTriangularSetGcdPackage(R,E,V,P,$) + regsetdecomppack ==> SquareFreeRegularSetDecompositionPackage(R,E,V,P,$) + + Exports == SquareFreeRegularTriangularSetCategory(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} + ++ from \spadtype{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$SREGSET: #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 SREGSET : bad #1" + ts := eif::$ + lp := rest lp + ts + + extendIfCan(ts:$,p:P) == + ground? p => "failed"::Union($,"failed") + empty? ts => + p := squareFreePart primitivePart p + (per([p]))::Union($,"failed") + not (mvar(ts) < mvar(p)) => "failed"::Union($,"failed") + invertible?(init(p),ts)@Boolean => + lts: Split := augment(p,ts) + #lts ~= 1 => "failed"::Union($,"failed") + (first lts)::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$SREGSET: 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 + lts: Split + if sqfr? + then + lts: Split := [] + lsfp := squareFreeFactors(p)$polsetpack + for f in lsfp repeat + (ground? f) or (mvar(f) < v) => "leave" + lpwt := squareFreePart(f,ts_v_-) + for pwt in lpwt repeat + sfp := pwt.val; us := pwt.tower + lts := cons( per(cons(pwt.val, rep(pwt.tower))), lts) + 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$SREGSET: ground? #1" + algebraic?(mvar(p),ts) => error "in augment$SREGSET: 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$SREGSET: ground? #1" + v := mvar(p) + not (mvar(ts) < mvar(p)) => error "in extend$SREGSET: bad #1" + split: List($) := invertibleSet(init(p),ts) + lts: List($) := [] + for us in split repeat + lts := concat(augment(p,us),lts) + lts + + invertible?(p:P,ts:$): Boolean == + stoseInvertible?(p,ts)$regsetgcdpack + + invertible?(p:P,ts:$): List BWT == + stoseInvertible?_sqfreg(p,ts)$regsetgcdpack + + invertibleSet(p:P,ts:$): Split == + stoseInvertibleSet_sqfreg(p,ts)$regsetgcdpack + + lastSubResultant(p1:P,p2:P,ts:$): List PWT == + stoseLastSubResultant(p1,p2,ts)$regsetgcdpack + + squareFreePart(p:P, ts: $): List PWT == + stoseSquareFreePart(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) + + 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 not 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 SFRTCAT SquareFreeRegularTriangularSetCategory>> +<<package SFQCMPK SquareFreeQuasiComponentPackage>> +<<package SFRGCD SquareFreeRegularTriangularSetGcdPackage>> +<<package SRDCMPK SquareFreeRegularSetDecompositionPackage>> +<<domain SREGSET SquareFreeRegularTriangularSet>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |