path: root/src/algebra/zerodim.spad.pamphlet

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{src/algebra zerodim.spad}
\author{Marc Moreno Maza}
\maketitle

\begin{abstract}
\end{abstract}
\tableofcontents
\eject

\section{package FGLMICPK FGLMIfCanPackage}

<<package FGLMICPK FGLMIfCanPackage>>=
import GcdDomain
import Symbol
import List
)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>>=
import GcdDomain
import RegularTriangularSet
import Symbol
import Boolean
import List
)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>>=
import GcdDomain
import Symbol
import List
)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: local 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: local 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>>=
import EuclideanDomain
import CharacteristicZero
import OrderedAbelianMonoidSup
import OrderedSet
import RecursivePolynomialCategory
import SquareFreeRegularTriangularSetCategory
import List
)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 := members(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 := members(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))) =>
         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 := members(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 negative? lambda 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) =>            
           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)) =>
           ts := internalAugment(p,ts)
           wip := [lp,ts]
           toSee := cons(wip,toSee)
         lts := makeLinearAndMonic(p,xi,ts,univ?,false,false)
         for ts: local in lts repeat
           wip := [lp,ts]
           toSee := cons(wip,toSee)
       toSave

@

\section{package RURPK RationalUnivariateRepresentationPackage}

<<package RURPK RationalUnivariateRepresentationPackage>>=
import EuclideanDomain
import CharacteristicZero
import Symbol
import List
)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: local 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 members(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>>=
import OrderedRing
import EuclideanDomain
import CharacteristicZero
import RealConstant
)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 members(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) := members(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 members(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 members(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}