diff options
| author | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
|---|---|---|
| committer | dos-reis <gdr@axiomatics.org> | 2007-08-14 05:14:52 +0000 |
| commit | ab8cc85adde879fb963c94d15675783f2cf4b183 (patch) | |
| tree | c202482327f474583b750b2c45dedfc4e4312b1d /src/algebra/zerodim.spad.pamphlet | |
| download | open-axiom-ab8cc85adde879fb963c94d15675783f2cf4b183.tar.gz | |
Initial population.
Diffstat (limited to 'src/algebra/zerodim.spad.pamphlet')
| -rw-r--r-- | src/algebra/zerodim.spad.pamphlet | 1189 |
1 files changed, 1189 insertions, 0 deletions
diff --git a/src/algebra/zerodim.spad.pamphlet b/src/algebra/zerodim.spad.pamphlet new file mode 100644 index 00000000..6b1db2a9 --- /dev/null +++ b/src/algebra/zerodim.spad.pamphlet @@ -0,0 +1,1189 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra zerodim.spad} +\author{Marc Moreno Maza} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{package FGLMICPK FGLMIfCanPackage} +<<package FGLMICPK FGLMIfCanPackage>>= +)abbrev package FGLMICPK FGLMIfCanPackage +++ Author: Marc Moreno Maza +++ Date Created: 08/02/1999 +++ Date Last Updated: 08/02/1999 +++ Description: +++ This is just an interface between several packages and domains. +++ The goal is to compute lexicographical Groebner bases +++ of sets of polynomial with type \spadtype{Polynomial R} +++ by the {\em FGLM} algorithm if this is possible (i.e. +++ if the input system generates a zero-dimensional ideal). +++ Version: 1. +FGLMIfCanPackage(R,ls): Exports == Implementation where + R: GcdDomain + ls: List Symbol + V ==> OrderedVariableList ls + N ==> NonNegativeInteger + Z ==> Integer + B ==> Boolean + Q1 ==> Polynomial R + Q2 ==> HomogeneousDistributedMultivariatePolynomial(ls,R) + Q3 ==> DistributedMultivariatePolynomial(ls,R) + E2 ==> HomogeneousDirectProduct(#ls,NonNegativeInteger) + E3 ==> DirectProduct(#ls,NonNegativeInteger) + poltopol ==> PolToPol(ls, R) + lingrobpack ==> LinGroebnerPackage(ls,R) + groebnerpack2 ==> GroebnerPackage(R,E2,V,Q2) + groebnerpack3 ==> GroebnerPackage(R,E3,V,Q3) + Exports == with + + zeroDimensional?: List(Q1) -> B + ++ \axiom{zeroDimensional?(lq1)} returns true iff + ++ \axiom{lq1} generates a zero-dimensional ideal + ++ w.r.t. the variables of \axiom{ls}. + fglmIfCan: List(Q1) -> Union(List(Q1), "failed") + ++ \axiom{fglmIfCan(lq1)} returns the lexicographical Groebner + ++ basis of \axiom{lq1} by using the {\em FGLM} strategy, + ++ if \axiom{zeroDimensional?(lq1)} holds. + groebner: List(Q1) -> List(Q1) + ++ \axiom{groebner(lq1)} returns the lexicographical Groebner + ++ basis of \axiom{lq1}. If \axiom{lq1} generates a zero-dimensional + ++ ideal then the {\em FGLM} strategy is used, otherwise + ++ the {\em Sugar} strategy is used. + + Implementation == add + + zeroDim?(lq2: List Q2): Boolean == + lq2 := groebner(lq2)$groebnerpack2 + empty? lq2 => false + #lq2 < #ls => false + lv: List(V) := [(variable(s)$V)::V for s in ls] + for q2 in lq2 while not empty?(lv) repeat + m := leadingMonomial(q2) + x := mainVariable(m)::V + if ground?(leadingCoefficient(univariate(m,x))) then + lv := remove(x, lv) + empty? lv + + zeroDimensional?(lq1: List(Q1)): Boolean == + lq2: List(Q2) := [pToHdmp(q1)$poltopol for q1 in lq1] + zeroDim?(lq2) + + fglmIfCan(lq1:List(Q1)): Union(List(Q1),"failed") == + lq2: List(Q2) := [pToHdmp(q1)$poltopol for q1 in lq1] + lq2 := groebner(lq2)$groebnerpack2 + not zeroDim?(lq2) => "failed"::Union(List(Q1),"failed") + lq3: List(Q3) := totolex(lq2)$lingrobpack + lq1 := [dmpToP(q3)$poltopol for q3 in lq3] + lq1::Union(List(Q1),"failed") + + groebner(lq1:List(Q1)): List(Q1) == + lq2: List(Q2) := [pToHdmp(q1)$poltopol for q1 in lq1] + lq2 := groebner(lq2)$groebnerpack2 + not zeroDim?(lq2) => + lq3: List(Q3) := [pToDmp(q1)$poltopol for q1 in lq1] + lq3 := groebner(lq3)$groebnerpack3 + [dmpToP(q3)$poltopol for q3 in lq3] + lq3: List(Q3) := totolex(lq2)$lingrobpack + [dmpToP(q3)$poltopol for q3 in lq3] + +@ +\section{domain RGCHAIN RegularChain} +<<domain RGCHAIN RegularChain>>= +)abbrev domain RGCHAIN RegularChain +++ Author: Marc Moreno Maza +++ Date Created: 01/1999 +++ Date Last Updated: 23/01/1999 +++ Description: +++ A domain for regular chains (i.e. regular triangular sets) over +++ a Gcd-Domain and with a fix list of variables. +++ This is just a front-end for the \spadtype{RegularTriangularSet} +++ domain constructor. +++ Version: 1. + +RegularChain(R,ls): Exports == Implementation where + R : GcdDomain + ls: List Symbol + V ==> OrderedVariableList ls + E ==> IndexedExponents V + P ==> NewSparseMultivariatePolynomial(R,V) + TS ==> RegularTriangularSet(R,E,V,P) + + Exports == RegularTriangularSetCategory(R,E,V,P) with + zeroSetSplit: (List P, Boolean, Boolean) -> List $ + ++ \spad{zeroSetSplit(lp,clos?,info?)} returns a list \spad{lts} of regular + ++ chains such that the union of the closures of their regular zero sets + ++ equals the affine variety associated with \spad{lp}. Moreover, + ++ if \spad{clos?} is \spad{false} then the union of the regular zero + ++ set of the \spad{ts} (for \spad{ts} in \spad{lts}) equals this variety. + ++ If \spad{info?} is \spad{true} then some information is + ++ displayed during the computations. See + ++ \axiomOpFrom{zeroSetSplit}{RegularTriangularSet}. + + Implementation == RegularTriangularSet(R,E,V,P) + +@ +\section{package LEXTRIPK LexTriangularPackage} +<<package LEXTRIPK LexTriangularPackage>>= +)abbrev package LEXTRIPK LexTriangularPackage +++ Author: Marc Moreno Maza +++ Date Created: 08/02/1999 +++ Date Last Updated: 08/02/1999 +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ Description: +++ A package for solving polynomial systems with finitely many solutions. +++ The decompositions are given by means of regular triangular sets. +++ The computations use lexicographical Groebner bases. +++ The main operations are \axiomOpFrom{lexTriangular}{LexTriangularPackage} +++ and \axiomOpFrom{squareFreeLexTriangular}{LexTriangularPackage}. +++ The second one provide decompositions by means of square-free regular triangular sets. +++ Both are based on the {\em lexTriangular} method described in [1]. +++ They differ from the algorithm described in [2] by the fact that +++ multiciplities of the roots are not kept. +++ With the \axiomOpFrom{squareFreeLexTriangular}{LexTriangularPackage} operation +++ all multiciplities are removed. With the other operation some multiciplities may remain. +++ Both operations admit an optional argument to produce normalized triangular sets. \newline +++ References: \newline +++ [1] D. LAZARD "Solving Zero-dimensional Algebraic Systems" +++ published in the J. of Symbol. Comput. (1992) 13, 117-131.\newline +++ [2] M. MORENO MAZA and R. RIOBOO "Computations of gcd over +++ algebraic towers of simple extensions" In proceedings of AAECC11, Paris, 1995.\newline +++ Version: 2. + +LexTriangularPackage(R,ls): Exports == Implementation where + + R: GcdDomain + ls: List Symbol + V ==> OrderedVariableList ls + E ==> IndexedExponents V + P ==> NewSparseMultivariatePolynomial(R,V) + TS ==> RegularChain(R,ls) + ST ==> SquareFreeRegularTriangularSet(R,E,V,P) + Q1 ==> Polynomial R + PS ==> GeneralPolynomialSet(R,E,V,P) + N ==> NonNegativeInteger + Z ==> Integer + B ==> Boolean + S ==> String + K ==> Fraction R + LP ==> List P + BWTS ==> Record(val : Boolean, tower : TS) + LpWTS ==> Record(val : (List P), tower : TS) + BWST ==> Record(val : Boolean, tower : ST) + LpWST ==> Record(val : (List P), tower : ST) + polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P) + quasicomppackTS ==> QuasiComponentPackage(R,E,V,P,TS) + regsetgcdpackTS ==> SquareFreeRegularTriangularSetGcdPackage(R,E,V,P,TS) + normalizpackTS ==> NormalizationPackage(R,E,V,P,TS) + quasicomppackST ==> QuasiComponentPackage(R,E,V,P,ST) + regsetgcdpackST ==> SquareFreeRegularTriangularSetGcdPackage(R,E,V,P,ST) + normalizpackST ==> NormalizationPackage(R,E,V,P,ST) + + Exports == with + + zeroDimensional?: LP -> B + ++ \axiom{zeroDimensional?(lp)} returns true iff + ++ \axiom{lp} generates a zero-dimensional ideal + ++ w.r.t. the variables involved in \axiom{lp}. + fglmIfCan: LP -> Union(LP, "failed") + ++ \axiom{fglmIfCan(lp)} returns the lexicographical Groebner + ++ basis of \axiom{lp} by using the {\em FGLM} strategy, + ++ if \axiom{zeroDimensional?(lp)} holds . + groebner: LP -> LP + ++ \axiom{groebner(lp)} returns the lexicographical Groebner + ++ basis of \axiom{lp}. If \axiom{lp} generates a zero-dimensional + ++ ideal then the {\em FGLM} strategy is used, otherwise + ++ the {\em Sugar} strategy is used. + lexTriangular: (LP, B) -> List TS + ++ \axiom{lexTriangular(base, norm?)} decomposes the variety + ++ associated with \axiom{base} into regular chains. + ++ Thus a point belongs to this variety iff it is a regular + ++ zero of a regular set in in the output. + ++ Note that \axiom{base} needs to be a lexicographical Groebner basis + ++ of a zero-dimensional ideal. If \axiom{norm?} is \axiom{true} + ++ then the regular sets are normalized. + squareFreeLexTriangular: (LP, B) -> List ST + ++ \axiom{squareFreeLexTriangular(base, norm?)} decomposes the variety + ++ associated with \axiom{base} into square-free regular chains. + ++ Thus a point belongs to this variety iff it is a regular + ++ zero of a regular set in in the output. + ++ Note that \axiom{base} needs to be a lexicographical Groebner basis + ++ of a zero-dimensional ideal. If \axiom{norm?} is \axiom{true} + ++ then the regular sets are normalized. + zeroSetSplit: (LP, B) -> List TS + ++ \axiom{zeroSetSplit(lp, norm?)} decomposes the variety + ++ associated with \axiom{lp} into regular chains. + ++ Thus a point belongs to this variety iff it is a regular + ++ zero of a regular set in in the output. + ++ Note that \axiom{lp} needs to generate a zero-dimensional ideal. + ++ If \axiom{norm?} is \axiom{true} then the regular sets are normalized. + zeroSetSplit: (LP, B) -> List ST + ++ \axiom{zeroSetSplit(lp, norm?)} decomposes the variety + ++ associated with \axiom{lp} into square-free regular chains. + ++ Thus a point belongs to this variety iff it is a regular + ++ zero of a regular set in in the output. + ++ Note that \axiom{lp} needs to generate a zero-dimensional ideal. + ++ If \axiom{norm?} is \axiom{true} then the regular sets are normalized. + + Implementation == add + + trueVariables(lp: List(P)): List Symbol == + lv: List V := variables([lp]$PS) + truels: List Symbol := [] + for s in ls repeat + if member?(variable(s)::V, lv) then truels := cons(s,truels) + reverse truels + + zeroDimensional?(lp:List(P)): Boolean == + truels: List Symbol := trueVariables(lp) + fglmpack := FGLMIfCanPackage(R,truels) + lq1: List(Q1) := [p::Q1 for p in lp] + zeroDimensional?(lq1)$fglmpack + + fglmIfCan(lp:List(P)): Union(List(P), "failed") == + truels: List Symbol := trueVariables(lp) + fglmpack := FGLMIfCanPackage(R,truels) + lq1: List(Q1) := [p::Q1 for p in lp] + foo := fglmIfCan(lq1)$fglmpack + foo case "failed" => return("failed" :: Union(List(P), "failed")) + lp := [retract(q1)$P for q1 in (foo :: List(Q1))] + lp::Union(List(P), "failed") + + groebner(lp:List(P)): List(P) == + truels: List Symbol := trueVariables(lp) + fglmpack := FGLMIfCanPackage(R,truels) + lq1: List(Q1) := [p::Q1 for p in lp] + lq1 := groebner(lq1)$fglmpack + lp := [retract(q1)$P for q1 in lq1] + + lexTriangular(base: List(P), norm?: Boolean): List(TS) == + base := sort(infRittWu?,base) + base := remove(zero?, base) + any?(ground?, base) => [] + ts: TS := empty() + toSee: List LpWTS := [[base,ts]$LpWTS] + toSave: List TS := [] + 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; v := mvar(p) + algebraic?(v,ts) => + error "lexTriangular$LEXTRIPK: should never happen !" + norm? and zero? remainder(init(p),ts).polnum => + toSee := cons([lp, ts]$LpWTS, toSee) + (not norm?) and zero? (initiallyReduce(init(p),ts)) => + toSee := cons([lp, ts]$LpWTS, toSee) + lbwt: List BWTS := invertible?(init(p),ts)$TS + while (not empty? lbwt) repeat + bwt := first lbwt; lbwt := rest lbwt + b := bwt.val; us := bwt.tower + (not b) => toSee := cons([lp, us], toSee) + lus: List TS + if norm? + then + newp := normalizedAssociate(p,us)$normalizpackTS + lus := [internalAugment(newp,us)$TS] + else + newp := p + lus := augment(newp,us)$TS + newlp := lp + while (not empty? newlp) and (mvar(first newlp) = v) repeat + newlp := rest newlp + for us in lus repeat + toSee := cons([newlp, us]$LpWTS, toSee) + algebraicSort(toSave)$quasicomppackTS + + zeroSetSplit(lp:List(P), norm?:B): List TS == + bar := fglmIfCan(lp) + bar case "failed" => + error "zeroSetSplit$LEXTRIPK: #1 not zero-dimensional" + lexTriangular(bar::(List P),norm?) + + squareFreeLexTriangular(base: List(P), norm?: Boolean): List(ST) == + base := sort(infRittWu?,base) + base := remove(zero?, base) + any?(ground?, base) => [] + ts: ST := empty() + toSee: List LpWST := [[base,ts]$LpWST] + toSave: List ST := [] + 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; v := mvar(p) + algebraic?(v,ts) => + error "lexTriangular$LEXTRIPK: should never happen !" + norm? and zero? remainder(init(p),ts).polnum => + toSee := cons([lp, ts]$LpWST, toSee) + (not norm?) and zero? (initiallyReduce(init(p),ts)) => + toSee := cons([lp, ts]$LpWST, toSee) + lbwt: List BWST := invertible?(init(p),ts)$ST + while (not empty? lbwt) repeat + bwt := first lbwt; lbwt := rest lbwt + b := bwt.val; us := bwt.tower + (not b) => toSee := cons([lp, us], toSee) + lus: List ST + if norm? + then + newp := normalizedAssociate(p,us)$normalizpackST + lus := augment(newp,us)$ST + else + lus := augment(p,us)$ST + newlp := lp + while (not empty? newlp) and (mvar(first newlp) = v) repeat + newlp := rest newlp + for us in lus repeat + toSee := cons([newlp, us]$LpWST, toSee) + algebraicSort(toSave)$quasicomppackST + + zeroSetSplit(lp:List(P), norm?:B): List ST == + bar := fglmIfCan(lp) + bar case "failed" => + error "zeroSetSplit$LEXTRIPK: #1 not zero-dimensional" + squareFreeLexTriangular(bar::(List P),norm?) + +@ +\section{package IRURPK InternalRationalUnivariateRepresentationPackage} +<<package IRURPK InternalRationalUnivariateRepresentationPackage>>= +)abbrev package IRURPK InternalRationalUnivariateRepresentationPackage +++ Author: Marc Moreno Maza +++ Date Created: 01/1999 +++ Date Last Updated: 23/01/1999 +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ Description: +++ An internal package for computing the rational univariate representation +++ of a zero-dimensional algebraic variety given by a square-free +++ triangular set. +++ The main operation is \axiomOpFrom{rur}{InternalRationalUnivariateRepresentationPackage}. +++ It is based on the {\em generic} algorithm description in [1]. \newline References: +++ [1] D. LAZARD "Solving Zero-dimensional Algebraic Systems" +++ Journal of Symbolic Computation, 1992, 13, 117-131 +++ Version: 1. + +InternalRationalUnivariateRepresentationPackage(R,E,V,P,TS): Exports == Implementation where + R : Join(EuclideanDomain,CharacteristicZero) + E : OrderedAbelianMonoidSup + V : OrderedSet + P : RecursivePolynomialCategory(R,E,V) + TS : SquareFreeRegularTriangularSetCategory(R,E,V,P) + N ==> NonNegativeInteger + Z ==> Integer + B ==> Boolean + LV ==> List V + LP ==> List P + PWT ==> Record(val: P, tower: TS) + LPWT ==> Record(val: LP, tower: TS) + WIP ==> Record(pol: P, gap: Z, tower: TS) + BWT ==> Record(val:Boolean, tower: TS) + polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,P) + normpack ==> NormalizationPackage(R,E,V,P,TS) + + Exports == with + + rur: (TS,B) -> List TS + ++ \spad{rur(ts,univ?)} returns a rational univariate representation + ++ of \spad{ts}. This assumes that the lowest polynomial in \spad{ts} + ++ is a variable \spad{v} which does not occur in the other polynomials + ++ of \spad{ts}. This variable will be used to define the simple + ++ algebraic extension over which these other polynomials will be + ++ rewritten as univariate polynomials with degree one. + ++ If \spad{univ?} is \spad{true} then these polynomials will have + ++ a constant initial. + checkRur: (TS, List TS) -> Boolean + ++ \spad{checkRur(ts,lus)} returns \spad{true} if \spad{lus} + ++ is a rational univariate representation of \spad{ts}. + + Implementation == add + + checkRur(ts: TS, lts: List TS): Boolean == + f0 := last(ts)::P + z := mvar(f0) + ts := collectUpper(ts,z) + dts: N := degree(ts) + lp := parts(ts) + dlts: N := 0 + for us in lts repeat + dlts := dlts + degree(us) + rems := [removeZero(p,us) for p in lp] + not every?(zero?,rems) => + output(us::OutputForm)$OutputPackage + return false + (dts =$N dlts)@Boolean + + convert(p:P,sqfr?:B):TS == + -- if sqfr? ASSUME p is square-free + newts: TS := empty() + sqfr? => internalAugment(p,newts) + p := squareFreePart(p) + internalAugment(p,newts) + + prepareRur(ts: TS): List LPWT == + not purelyAlgebraic?(ts)$TS => + error "rur$IRURPK: #1 is not zero-dimensional" + lp: LP := parts(ts)$TS + lp := sort(infRittWu?,lp) + empty? lp => + error "rur$IRURPK: #1 is empty" + f0 := first lp; lp := rest lp +-- not (one?(init(f0)) and one?(mdeg(f0)) and zero?(tail(f0))) => + not ((init(f0) = 1) and (mdeg(f0) = 1) and zero?(tail(f0))) => + error "rur$IRURPK: #1 has no generating root." + empty? lp => + error "rur$IRURPK: #1 has a generating root but no indeterminates" + z: V := mvar(f0) + f1 := first lp; lp := rest lp + x1: V := mvar(f1) + newf1 := x1::P - z::P + toSave: List LPWT := [] + for ff1 in irreducibleFactors([f1])$polsetpack repeat + newf0 := eval(ff1,mvar(f1),f0) + ts := internalAugment(newf1,convert(newf0,true)@TS) + toSave := cons([lp,ts],toSave) + toSave + + makeMonic(z:V,c:P,r:P,ts:TS,s:P,univ?:B): TS == + --ASSUME r is a irreducible univariate polynomial in z + --ASSUME c and s only depends on z and mvar(s) + --ASSUME c and a have main degree 1 + --ASSUME c and s have a constant initial + --ASSUME mvar(ts) < mvar(s) + lp: LP := parts(ts) + lp := sort(infRittWu?,lp) + newts: TS := convert(r,true)@TS + s := remainder(s,newts).polnum + if univ? + then + s := normalizedAssociate(s,newts)$normpack + for p in lp repeat + p := lazyPrem(eval(p,z,c),s) + p := remainder(p,newts).polnum + newts := internalAugment(p,newts) + internalAugment(s,newts) + + next(lambda:Z):Z == + if lambda < 0 then lambda := - lambda + 1 else lambda := - lambda + + makeLinearAndMonic(p: P, xi: V, ts: TS, univ?:B, check?: B, info?: B): List TS == + -- if check? THEN some VERIFICATIONS are performed + -- if info? THEN some INFORMATION is displayed + f0 := last(ts)::P + z: V := mvar(f0) + lambda: Z := 1 + ts := collectUpper(ts,z) + toSee: List WIP := [[f0,lambda,ts]$WIP] + toSave: List TS := [] + while not empty? toSee repeat + wip := first toSee; toSee := rest toSee + (f0, lambda, ts) := (wip.pol, wip.gap, wip.tower) + if check? and ((not univariate?(f0)$polsetpack) or (mvar(f0) ~= z)) + then + output("Bad f0: ")$OutputPackage + output(f0::OutputForm)$OutputPackage + c: P := lambda * xi::P + z::P + f := eval(f0,z,c); q := eval(p,z,c) + prs := subResultantChain(q,f) + r := first prs; prs := rest prs + check? and ((not zero? degree(r,xi)) or (empty? prs)) => + error "rur$IRURPK: should never happen !" + s := first prs; prs := rest prs + check? and (zero? degree(s,xi)) and (empty? prs) => + error "rur$IRURPK: should never happen !!" + if zero? degree(s,xi) then s := first prs +-- not one? degree(s,xi) => + not (degree(s,xi) = 1) => + toSee := cons([f0,next(lambda),ts]$WIP,toSee) + h := init(s) + r := squareFreePart(r) + ground?(h) or ground?(gcd(h,r)) => + for fr in irreducibleFactors([r])$polsetpack repeat + ground? fr => "leave" + toSave := cons(makeMonic(z,c,fr,ts,s,univ?),toSave) + if info? + then + output("Unlucky lambda")$OutputPackage + output(h::OutputForm)$OutputPackage + output(r::OutputForm)$OutputPackage + toSee := cons([f0,next(lambda),ts]$WIP,toSee) + toSave + + rur (ts: TS,univ?:Boolean): List TS == + toSee: List LPWT := prepareRur(ts) + toSave: List TS := [] + while not empty? toSee repeat + wip := first toSee; toSee := rest toSee + ts: TS := wip.tower + lp: LP := wip.val + empty? lp => toSave := cons(ts,toSave) + p := first lp; lp := rest lp + xi: V := mvar(p) + p := remainder(p,ts).polnum + if not univ? + then + p := primitivePart stronglyReduce(p,ts) + ground?(p) or (mvar(p) < xi) => + error "rur$IRUROK: should never happen" +-- (one? mdeg(p)) and (ground? init(p)) => + (mdeg(p) = 1) and (ground? init(p)) => + ts := internalAugment(p,ts) + wip := [lp,ts] + toSee := cons(wip,toSee) + lts := makeLinearAndMonic(p,xi,ts,univ?,false,false) + for ts in lts repeat + wip := [lp,ts] + toSee := cons(wip,toSee) + toSave + +@ +\section{package RURPK RationalUnivariateRepresentationPackage} +<<package RURPK RationalUnivariateRepresentationPackage>>= +)abbrev package RURPK RationalUnivariateRepresentationPackage +++ Author: Marc Moreno Maza +++ Date Created: 01/1999 +++ Date Last Updated: 23/01/1999 +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Description: +++ A package for computing the rational univariate representation +++ of a zero-dimensional algebraic variety given by a regular +++ triangular set. This package is essentially an interface for the +++ \spadtype{InternalRationalUnivariateRepresentationPackage} constructor. +++ It is used in the \spadtype{ZeroDimensionalSolvePackage} +++ for solving polynomial systems with finitely many solutions. +++ Version: 1. + +RationalUnivariateRepresentationPackage(R,ls): Exports == Implementation where + R : Join(EuclideanDomain,CharacteristicZero) + ls: List Symbol + N ==> NonNegativeInteger + Z ==> Integer + P ==> Polynomial R + LP ==> List P + U ==> SparseUnivariatePolynomial(R) + RUR ==> Record(complexRoots: U, coordinates: LP) + + Exports == with + + rur: (LP,Boolean) -> List RUR + ++ \spad{rur(lp,univ?)} returns a rational univariate representation + ++ of \spad{lp}. This assumes that \spad{lp} defines a regular + ++ triangular \spad{ts} whose associated variety is zero-dimensional + ++ over \spad{R}. \spad{rur(lp,univ?)} returns a list of items + ++ \spad{[u,lc]} where \spad{u} is an irreducible univariate polynomial + ++ and each \spad{c} in \spad{lc} involves two variables: one from \spad{ls}, + ++ called the coordinate of \spad{c}, and an extra variable which + ++ represents any root of \spad{u}. Every root of \spad{u} leads to + ++ a tuple of values for the coordinates of \spad{lc}. Moreover, + ++ a point \spad{x} belongs to the variety associated with \spad{lp} iff + ++ there exists an item \spad{[u,lc]} in \spad{rur(lp,univ?)} and + ++ a root \spad{r} of \spad{u} such that \spad{x} is given by the + ++ tuple of values for the coordinates of \spad{lc} evaluated at \spad{r}. + ++ If \spad{univ?} is \spad{true} then each polynomial \spad{c} + ++ will have a constant leading coefficient w.r.t. its coordinate. + ++ See the example which illustrates the \spadtype{ZeroDimensionalSolvePackage} + ++ package constructor. + rur: (LP) -> List RUR + ++ \spad{rur(lp)} returns the same as \spad{rur(lp,true)} + rur: (LP,Boolean,Boolean) -> List RUR + ++ \spad{rur(lp,univ?,check?)} returns the same as \spad{rur(lp,true)}. + ++ Moreover, if \spad{check?} is \spad{true} then the result is checked. + + Implementation == add + news: Symbol := new()$Symbol + lv: List Symbol := concat(ls,news) + V ==> OrderedVariableList(lv) + Q ==> NewSparseMultivariatePolynomial(R,V) + E ==> IndexedExponents V + TS ==> SquareFreeRegularTriangularSet(R,E,V,Q) + QWT ==> Record(val: Q, tower: TS) + LQWT ==> Record(val: List Q, tower: TS) + polsetpack ==> PolynomialSetUtilitiesPackage(R,E,V,Q) + normpack ==> NormalizationPackage(R,E,V,Q,TS) + rurpack ==> InternalRationalUnivariateRepresentationPackage(R,E,V,Q,TS) + newv: V := variable(news)::V + newq : Q := newv :: Q + + rur(lp: List P, univ?: Boolean, check?: Boolean): List RUR == + lp := remove(zero?,lp) + empty? lp => + error "rur$RURPACK: #1 is empty" + any?(ground?,lp) => + error "rur$RURPACK: #1 is not a triangular set" + ts: TS := [[newq]$(List Q)] + lq: List Q := [] + for p in lp repeat + rif: Union(Q,"failed") := retractIfCan(p)$Q + rif case "failed" => + error "rur$RURPACK: #1 is not a subset of R[ls]" + q: Q := rif::Q + lq := cons(q,lq) + lq := sort(infRittWu?,lq) + toSee: List LQWT := [[lq,ts]$LQWT] + toSave: List TS := [] + while not empty? toSee repeat + lqwt := first toSee; toSee := rest toSee + lq := lqwt.val; ts := lqwt.tower + empty? lq => + -- output(ts::OutputForm)$OutputPackage + toSave := cons(ts,toSave) + q := first lq; lq := rest lq + not (mvar(q) > mvar(ts)) => + error "rur$RURPACK: #1 is not a triangular set" + empty? (rest(ts)::TS) => + lfq := irreducibleFactors([q])$polsetpack + for fq in lfq repeat + newts := internalAugment(fq,ts) + newlq := [remainder(q,newts).polnum for q in lq] + toSee := cons([newlq,newts]$LQWT,toSee) + lsfqwt: List QWT := squareFreePart(q,ts) + for qwt in lsfqwt repeat + q := qwt.val; ts := qwt.tower + if not ground? init(q) + then + q := normalizedAssociate(q,ts)$normpack + newts := internalAugment(q,ts) + newlq := [remainder(q,newts).polnum for q in lq] + toSee := cons([newlq,newts]$LQWT,toSee) + toReturn: List RUR := [] + for ts in toSave repeat + lus := rur(ts,univ?)$rurpack + check? and (not checkRur(ts,lus)$rurpack) => + output("RUR for: ")$OutputPackage + output(ts::OutputForm)$OutputPackage + output("Is: ")$OutputPackage + for us in lus repeat output(us::OutputForm)$OutputPackage + error "rur$RURPACK: bad result with function rur$IRURPK" + for us in lus repeat + g: U := univariate(select(us,newv)::Q)$Q + lc: LP := [convert(q)@P for q in parts(collectUpper(us,newv))] + toReturn := cons([g,lc]$RUR, toReturn) + toReturn + + rur(lp: List P, univ?: Boolean): List RUR == + rur(lp,univ?,false) + + rur(lp: List P): List RUR == rur(lp,true) + +@ +\section{package ZDSOLVE ZeroDimensionalSolvePackage} +Based on triangular decompositions and the {\bf RealClosure} constructor, +the pacakge {\bf ZeroDimensionalSolvePackage} provides operations for +computing symbolically the real or complex roots of polynomial systems +with finitely many solutions. +<<package ZDSOLVE ZeroDimensionalSolvePackage>>= +)abbrev package ZDSOLVE ZeroDimensionalSolvePackage +++ Author: Marc Moreno Maza +++ Date Created: 23/01/1999 +++ Date Last Updated: 08/02/1999 +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ A package for computing symbolically the complex and real roots of +++ zero-dimensional algebraic systems over the integer or rational +++ numbers. Complex roots are given by means of univariate representations +++ of irreducible regular chains. Real roots are given by means of tuples +++ of coordinates lying in the \spadtype{RealClosure} of the coefficient ring. +++ This constructor takes three arguments. The first one \spad{R} is the +++ coefficient ring. The second one \spad{ls} is the list of variables involved +++ in the systems to solve. The third one must be \spad{concat(ls,s)} where +++ \spad{s} is an additional symbol used for the univariate representations. +++ WARNING: The third argument is not checked. +++ All operations are based on triangular decompositions. +++ The default is to compute these decompositions directly from the input +++ system by using the \spadtype{RegularChain} domain constructor. +++ The lexTriangular algorithm can also be used for computing these decompositions +++ (see the \spadtype{LexTriangularPackage} package constructor). +++ For that purpose, the operations \axiomOpFrom{univariateSolve}{ZeroDimensionalSolvePackage}, +++ \axiomOpFrom{realSolve}{ZeroDimensionalSolvePackage} and +++ \axiomOpFrom{positiveSolve}{ZeroDimensionalSolvePackage} admit an optional +++ argument. \newline Author: Marc Moreno Maza. + +++ Version: 1. + +ZeroDimensionalSolvePackage(R,ls,ls2): Exports == Implementation where + R : Join(OrderedRing,EuclideanDomain,CharacteristicZero,RealConstant) + ls: List Symbol + ls2: List Symbol + V ==> OrderedVariableList(ls) + N ==> NonNegativeInteger + Z ==> Integer + B ==> Boolean + P ==> Polynomial R + LP ==> List P + LS ==> List Symbol + Q ==> NewSparseMultivariatePolynomial(R,V) + U ==> SparseUnivariatePolynomial(R) + TS ==> RegularChain(R,ls) + RUR ==> Record(complexRoots: U, coordinates: LP) + K ==> Fraction R + RC ==> RealClosure(K) + PRC ==> Polynomial RC + REALSOL ==> List RC + URC ==> SparseUnivariatePolynomial RC + V2 ==> OrderedVariableList(ls2) + Q2 ==> NewSparseMultivariatePolynomial(R,V2) + E2 ==> IndexedExponents V2 + ST ==> SquareFreeRegularTriangularSet(R,E2,V2,Q2) + Q2WT ==> Record(val: Q2, tower: ST) + LQ2WT ==> Record(val: List(Q2), tower: ST) + WIP ==> Record(reals: List(RC), vars: List(Symbol), pols: List(Q2)) + polsetpack ==> PolynomialSetUtilitiesPackage(R,E2,V2,Q2) + normpack ==> NormalizationPackage(R,E2,V2,Q2,ST) + rurpack ==> InternalRationalUnivariateRepresentationPackage(R,E2,V2,Q2,ST) + quasicomppack ==> SquareFreeQuasiComponentPackage(R,E2,V2,Q2,ST) + lextripack ==> LexTriangularPackage(R,ls) + + Exports == with + triangSolve: (LP,B,B) -> List RegularChain(R,ls) + ++ \spad{triangSolve(lp,info?,lextri?)} decomposes the variety + ++ associated with \axiom{lp} into regular chains. + ++ Thus a point belongs to this variety iff it is a regular + ++ zero of a regular set in in the output. + ++ Note that \axiom{lp} needs to generate a zero-dimensional ideal. + ++ If \axiom{lp} is not zero-dimensional then the result is only + ++ a decomposition of its zero-set in the sense of the closure + ++ (w.r.t. Zarisky topology). + ++ Moreover, if \spad{info?} is \spad{true} then some information is + ++ displayed during the computations. + ++ See \axiomOpFrom{zeroSetSplit}{RegularTriangularSetCategory}(lp,true,info?). + ++ If \spad{lextri?} is \spad{true} then the lexTriangular algorithm is called + ++ from the \spadtype{LexTriangularPackage} constructor + ++ (see \axiomOpFrom{zeroSetSplit}{LexTriangularPackage}(lp,false)). + ++ Otherwise, the triangular decomposition is computed directly from the input + ++ system by using the \axiomOpFrom{zeroSetSplit}{RegularChain} from \spadtype{RegularChain}. + triangSolve: (LP,B) -> List RegularChain(R,ls) + ++ \spad{triangSolve(lp,info?)} returns the same as \spad{triangSolve(lp,false)} + triangSolve: LP -> List RegularChain(R,ls) + ++ \spad{triangSolve(lp)} returns the same as \spad{triangSolve(lp,false,false)} + univariateSolve: RegularChain(R,ls) -> List Record(complexRoots: U, coordinates: LP) + ++ \spad{univariateSolve(ts)} returns a univariate representation + ++ of \spad{ts}. + ++ See \axiomOpFrom{rur}{RationalUnivariateRepresentationPackage}(lp,true). + univariateSolve: (LP,B,B,B) -> List RUR + ++ \spad{univariateSolve(lp,info?,check?,lextri?)} returns a univariate + ++ representation of the variety associated with \spad{lp}. + ++ Moreover, if \spad{info?} is \spad{true} then some information is + ++ displayed during the decomposition into regular chains. + ++ If \spad{check?} is \spad{true} then the result is checked. + ++ See \axiomOpFrom{rur}{RationalUnivariateRepresentationPackage}(lp,true). + ++ If \spad{lextri?} is \spad{true} then the lexTriangular algorithm is called + ++ from the \spadtype{LexTriangularPackage} constructor + ++ (see \axiomOpFrom{zeroSetSplit}{LexTriangularPackage}(lp,false)). + ++ Otherwise, the triangular decomposition is computed directly from the input + ++ system by using the \axiomOpFrom{zeroSetSplit}{RegularChain} from \spadtype{RegularChain}. + univariateSolve: (LP,B,B) -> List RUR + ++ \spad{univariateSolve(lp,info?,check?)} returns the same as + ++ \spad{univariateSolve(lp,info?,check?,false)}. + univariateSolve: (LP,B) -> List RUR + ++ \spad{univariateSolve(lp,info?)} returns the same as + ++ \spad{univariateSolve(lp,info?,false,false)}. + univariateSolve: LP -> List RUR + ++ \spad{univariateSolve(lp)} returns the same as + ++ \spad{univariateSolve(lp,false,false,false)}. + realSolve: RegularChain(R,ls) -> List REALSOL + ++ \spad{realSolve(ts)} returns the set of the points in the regular + ++ zero set of \spad{ts} whose coordinates are all real. + ++ WARNING: For each set of coordinates given by \spad{realSolve(ts)} + ++ the ordering of the indeterminates is reversed w.r.t. \spad{ls}. + realSolve: (LP,B,B,B) -> List REALSOL + ++ \spad{realSolve(ts,info?,check?,lextri?)} returns the set of the points + ++ in the variety associated with \spad{lp} whose coordinates are all real. + ++ Moreover, if \spad{info?} is \spad{true} then some information is + ++ displayed during decomposition into regular chains. + ++ If \spad{check?} is \spad{true} then the result is checked. + ++ If \spad{lextri?} is \spad{true} then the lexTriangular algorithm is called + ++ from the \spadtype{LexTriangularPackage} constructor + ++ (see \axiomOpFrom{zeroSetSplit}{LexTriangularPackage}(lp,false)). + ++ Otherwise, the triangular decomposition is computed directly from the input + ++ system by using the \axiomOpFrom{zeroSetSplit}{RegularChain} from \spadtype{RegularChain}. + ++ WARNING: For each set of coordinates given by \spad{realSolve(ts,info?,check?,lextri?)} + ++ the ordering of the indeterminates is reversed w.r.t. \spad{ls}. + realSolve: (LP,B,B) -> List REALSOL + ++ \spad{realSolve(ts,info?,check?)} returns the same as \spad{realSolve(ts,info?,check?,false)}. + realSolve: (LP,B) -> List REALSOL + ++ \spad{realSolve(ts,info?)} returns the same as \spad{realSolve(ts,info?,false,false)}. + realSolve: LP -> List REALSOL + ++ \spad{realSolve(lp)} returns the same as \spad{realSolve(ts,false,false,false)} + positiveSolve: RegularChain(R,ls) -> List REALSOL + ++ \spad{positiveSolve(ts)} returns the points of the regular + ++ set of \spad{ts} with (real) strictly positive coordinates. + positiveSolve: (LP,B,B) -> List REALSOL + ++ \spad{positiveSolve(lp,info?,lextri?)} returns the set of the points + ++ in the variety associated with \spad{lp} whose coordinates are (real) strictly positive. + ++ Moreover, if \spad{info?} is \spad{true} then some information is + ++ displayed during decomposition into regular chains. + ++ If \spad{lextri?} is \spad{true} then the lexTriangular algorithm is called + ++ from the \spadtype{LexTriangularPackage} constructor + ++ (see \axiomOpFrom{zeroSetSplit}{LexTriangularPackage}(lp,false)). + ++ Otherwise, the triangular decomposition is computed directly from the input + ++ system by using the \axiomOpFrom{zeroSetSplit}{RegularChain} from \spadtype{RegularChain}. + ++ WARNING: For each set of coordinates given by \spad{positiveSolve(lp,info?,lextri?)} + ++ the ordering of the indeterminates is reversed w.r.t. \spad{ls}. + positiveSolve: (LP,B) -> List REALSOL + ++ \spad{positiveSolve(lp)} returns the same as \spad{positiveSolve(lp,info?,false)}. + positiveSolve: LP -> List REALSOL + ++ \spad{positiveSolve(lp)} returns the same as \spad{positiveSolve(lp,false,false)}. + squareFree: (TS) -> List ST + ++ \spad{squareFree(ts)} returns the square-free factorization of \spad{ts}. + ++ Moreover, each factor is a Lazard triangular set and the decomposition + ++ is a Kalkbrener split of \spad{ts}, which is enough here for + ++ the matter of solving zero-dimensional algebraic systems. + ++ WARNING: \spad{ts} is not checked to be zero-dimensional. + convert: Q -> Q2 + ++ \spad{convert(q)} converts \spad{q}. + convert: P -> PRC + ++ \spad{convert(p)} converts \spad{p}. + convert: Q2 -> PRC + ++ \spad{convert(q)} converts \spad{q}. + convert: U -> URC + ++ \spad{convert(u)} converts \spad{u}. + convert: ST -> List Q2 + ++ \spad{convert(st)} returns the members of \spad{st}. + + Implementation == add + news: Symbol := last(ls2)$(List Symbol) + newv: V2 := (variable(news)$V2)::V2 + newq: Q2 := newv :: Q2 + + convert(q:Q):Q2 == + ground? q => (ground(q))::Q2 + q2: Q2 := 0 + while not ground?(q) repeat + v: V := mvar(q) + d: N := mdeg(q) + v2: V2 := (variable(convert(v)@Symbol)$V2)::V2 + iq2: Q2 := convert(init(q))@Q2 + lq2: Q2 := (v2 :: Q2) + lq2 := lq2 ** d + q2 := iq2 * lq2 + q2 + q := tail(q) + q2 + (ground(q))::Q2 + + squareFree(ts:TS):List(ST) == + irred?: Boolean := false + st: ST := [[newq]$(List Q2)] + lq: List(Q2) := [convert(p)@Q2 for p in parts(ts)] + lq := sort(infRittWu?,lq) + toSee: List LQ2WT := [] + if irred? + then + lf := irreducibleFactors([first lq])$polsetpack + lq := rest lq + for f in lf repeat + toSee := cons([cons(f,lq),st]$LQ2WT, toSee) + else + toSee := [[lq,st]$LQ2WT] + toSave: List ST := [] + while not empty? toSee repeat + lqwt := first toSee; toSee := rest toSee + lq := lqwt.val; st := lqwt.tower + empty? lq => + toSave := cons(st,toSave) + q := first lq; lq := rest lq + lsfqwt: List Q2WT := squareFreePart(q,st)$ST + for sfqwt in lsfqwt repeat + q := sfqwt.val; st := sfqwt.tower + if not ground? init(q) + then + q := normalizedAssociate(q,st)$normpack + newts := internalAugment(q,st)$ST + newlq := [remainder(q,newts).polnum for q in lq] + toSee := cons([newlq,newts]$LQ2WT,toSee) + toSave + + + triangSolve(lp: LP, info?: B, lextri?: B): List TS == + lq: List(Q) := [convert(p)$Q for p in lp] + lextri? => zeroSetSplit(lq,false)$lextripack + zeroSetSplit(lq,true,info?)$TS + + triangSolve(lp: LP, info?: B): List TS == triangSolve(lp,info?,false) + + triangSolve(lp: LP): List TS == triangSolve(lp,false) + + convert(u: U): URC == + zero? u => 0 + ground? u => ((ground(u) :: K)::RC)::URC + uu: URC := 0 + while not ground? u repeat + uu := monomial((leadingCoefficient(u) :: K):: RC,degree(u)) + uu + u := reductum u + uu + ((ground(u) :: K)::RC)::URC + + coerceFromRtoRC(r:R): RC == + (r::K)::RC + + convert(p:P): PRC == + map(coerceFromRtoRC,p)$PolynomialFunctions2(R,RC) + + convert(q2:Q2): PRC == + p: P := coerce(q2)$Q2 + convert(p)@PRC + + convert(sts:ST): List Q2 == + lq2: List(Q2) := parts(sts)$ST + lq2 := sort(infRittWu?,lq2) + rest(lq2) + + realSolve(ts: TS): List REALSOL == + lsts: List ST := squareFree(ts) + lr: REALSOL := [] + lv: List Symbol := [] + toSee := [[lr,lv,convert(sts)@(List Q2)]$WIP for sts in lsts] + toSave: List REALSOL := [] + while not empty? toSee repeat + wip := first toSee; toSee := rest toSee + lr := wip.reals; lv := wip.vars; lq2 := wip.pols + (empty? lq2) and (not empty? lr) => + toSave := cons(reverse(lr),toSave) + q2 := first lq2; lq2 := rest lq2 + qrc := convert(q2)@PRC + if not empty? lr + then + for r in reverse(lr) for v in reverse(lv) repeat + qrc := eval(qrc,v,r) + lv := cons((mainVariable(qrc) :: Symbol),lv) + urc: URC := univariate(qrc)@URC + urcRoots := allRootsOf(urc)$RC + for urcRoot in urcRoots repeat + toSee := cons([cons(urcRoot,lr),lv,lq2]$WIP, toSee) + toSave + + realSolve(lp: List(P), info?:Boolean, check?:Boolean, lextri?: Boolean): List REALSOL == + lts: List TS + lq: List(Q) := [convert(p)$Q for p in lp] + if lextri? + then + lts := zeroSetSplit(lq,false)$lextripack + else + lts := zeroSetSplit(lq,true,info?)$TS + lsts: List ST := [] + for ts in lts repeat + lsts := concat(squareFree(ts), lsts) + lsts := removeSuperfluousQuasiComponents(lsts)$quasicomppack + lr: REALSOL := [] + lv: List Symbol := [] + toSee := [[lr,lv,convert(sts)@(List Q2)]$WIP for sts in lsts] + toSave: List REALSOL := [] + while not empty? toSee repeat + wip := first toSee; toSee := rest toSee + lr := wip.reals; lv := wip.vars; lq2 := wip.pols + (empty? lq2) and (not empty? lr) => + toSave := cons(reverse(lr),toSave) + q2 := first lq2; lq2 := rest lq2 + qrc := convert(q2)@PRC + if not empty? lr + then + for r in reverse(lr) for v in reverse(lv) repeat + qrc := eval(qrc,v,r) + lv := cons((mainVariable(qrc) :: Symbol),lv) + urc: URC := univariate(qrc)@URC + urcRoots := allRootsOf(urc)$RC + for urcRoot in urcRoots repeat + toSee := cons([cons(urcRoot,lr),lv,lq2]$WIP, toSee) + if check? + then + for p in lp repeat + for realsol in toSave repeat + prc: PRC := convert(p)@PRC + for rr in realsol for symb in reverse(ls) repeat + prc := eval(prc,symb,rr) + not zero? prc => + error "realSolve$ZDSOLVE: bad result" + toSave + + realSolve(lp: List(P), info?:Boolean, check?:Boolean): List REALSOL == + realSolve(lp,info?,check?,false) + + realSolve(lp: List(P), info?:Boolean): List REALSOL == + realSolve(lp,info?,false,false) + + realSolve(lp: List(P)): List REALSOL == + realSolve(lp,false,false,false) + + positiveSolve(ts: TS): List REALSOL == + lsts: List ST := squareFree(ts) + lr: REALSOL := [] + lv: List Symbol := [] + toSee := [[lr,lv,convert(sts)@(List Q2)]$WIP for sts in lsts] + toSave: List REALSOL := [] + while not empty? toSee repeat + wip := first toSee; toSee := rest toSee + lr := wip.reals; lv := wip.vars; lq2 := wip.pols + (empty? lq2) and (not empty? lr) => + toSave := cons(reverse(lr),toSave) + q2 := first lq2; lq2 := rest lq2 + qrc := convert(q2)@PRC + if not empty? lr + then + for r in reverse(lr) for v in reverse(lv) repeat + qrc := eval(qrc,v,r) + lv := cons((mainVariable(qrc) :: Symbol),lv) + urc: URC := univariate(qrc)@URC + urcRoots := allRootsOf(urc)$RC + for urcRoot in urcRoots repeat + if positive? urcRoot + then + toSee := cons([cons(urcRoot,lr),lv,lq2]$WIP, toSee) + toSave + + positiveSolve(lp: List(P), info?:Boolean, lextri?: Boolean): List REALSOL == + lts: List TS + lq: List(Q) := [convert(p)$Q for p in lp] + if lextri? + then + lts := zeroSetSplit(lq,false)$lextripack + else + lts := zeroSetSplit(lq,true,info?)$TS + lsts: List ST := [] + for ts in lts repeat + lsts := concat(squareFree(ts), lsts) + lsts := removeSuperfluousQuasiComponents(lsts)$quasicomppack + lr: REALSOL := [] + lv: List Symbol := [] + toSee := [[lr,lv,convert(sts)@(List Q2)]$WIP for sts in lsts] + toSave: List REALSOL := [] + while not empty? toSee repeat + wip := first toSee; toSee := rest toSee + lr := wip.reals; lv := wip.vars; lq2 := wip.pols + (empty? lq2) and (not empty? lr) => + toSave := cons(reverse(lr),toSave) + q2 := first lq2; lq2 := rest lq2 + qrc := convert(q2)@PRC + if not empty? lr + then + for r in reverse(lr) for v in reverse(lv) repeat + qrc := eval(qrc,v,r) + lv := cons((mainVariable(qrc) :: Symbol),lv) + urc: URC := univariate(qrc)@URC + urcRoots := allRootsOf(urc)$RC + for urcRoot in urcRoots repeat + if positive? urcRoot + then + toSee := cons([cons(urcRoot,lr),lv,lq2]$WIP, toSee) + toSave + + positiveSolve(lp: List(P), info?:Boolean): List REALSOL == + positiveSolve(lp, info?, false) + + positiveSolve(lp: List(P)): List REALSOL == + positiveSolve(lp, false, false) + + univariateSolve(ts: TS): List RUR == + toSee: List ST := squareFree(ts) + toSave: List RUR := [] + for st in toSee repeat + lus: List ST := rur(st,true)$rurpack + for us in lus repeat + g: U := univariate(select(us,newv)::Q2)$Q2 + lc: LP := [convert(q2)@P for q2 in parts(collectUpper(us,newv)$ST)$ST] + toSave := cons([g,lc]$RUR, toSave) + toSave + + univariateSolve(lp: List(P), info?:Boolean, check?:Boolean, lextri?: Boolean): List RUR == + lts: List TS + lq: List(Q) := [convert(p)$Q for p in lp] + if lextri? + then + lts := zeroSetSplit(lq,false)$lextripack + else + lts := zeroSetSplit(lq,true,info?)$TS + toSee: List ST := [] + for ts in lts repeat + toSee := concat(squareFree(ts), toSee) + toSee := removeSuperfluousQuasiComponents(toSee)$quasicomppack + toSave: List RUR := [] + if check? + then + lq2: List(Q2) := [convert(p)$Q2 for p in lp] + for st in toSee repeat + lus: List ST := rur(st,true)$rurpack + for us in lus repeat + if check? + then + rems: List(Q2) := [removeZero(q2,us)$ST for q2 in lq2] + not every?(zero?,rems) => + output(st::OutputForm)$OutputPackage + output("Has a bad RUR component:")$OutputPackage + output(us::OutputForm)$OutputPackage + error "univariateSolve$ZDSOLVE: bad RUR" + g: U := univariate(select(us,newv)::Q2)$Q2 + lc: LP := [convert(q2)@P for q2 in parts(collectUpper(us,newv)$ST)$ST] + toSave := cons([g,lc]$RUR, toSave) + toSave + + univariateSolve(lp: List(P), info?:Boolean, check?:Boolean): List RUR == + univariateSolve(lp,info?,check?,false) + + univariateSolve(lp: List(P), info?:Boolean): List RUR == + univariateSolve(lp,info?,false,false) + + univariateSolve(lp: List(P)): List RUR == + univariateSolve(lp,false,false,false) + +@ +\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>> + +<<package FGLMICPK FGLMIfCanPackage>> +<<domain RGCHAIN RegularChain>> +<<package LEXTRIPK LexTriangularPackage>> +<<package IRURPK InternalRationalUnivariateRepresentationPackage>> +<<package RURPK RationalUnivariateRepresentationPackage>> +<<package ZDSOLVE ZeroDimensionalSolvePackage>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |