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

\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra pleqn.spad}
\author{William Sit}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\begin{verbatim}
++ This package completely solves a parametric linear system of equations
++ by decomposing the set of all parametric values for which the linear
++ system is consistent into a union of quasi-algebraic  sets (which need
++ not be irredundant, but most of the time is). Each quasi-algebraic
++ set is described by a list of polynomials that vanish on the set, and
++ a list of polynomials that vanish at no point of the set.
++ For each quasi-algebraic set, the solution of the linear system
++ is given, as a particular solution and  a basis of the homogeneous
++ system.
++ The parametric linear system should be given in matrix form, with
++ a coefficient matrix and a right hand side vector. The entries
++ of the coefficient matrix and right hand side vector should be
++ polynomials in the parametric variables, over a Euclidean domain
++ of characteristic zero.
++
++ If the system is homogeneous, the right hand side need not be given.
++ The right hand side can also be  replaced by an indeterminate vector,
++ in which case, the conditions required for consistency will also be
++ given.
++
++ The package has other facilities for saving results to external
++ files, as well as solving the system for a specified minimum rank.
++ Altogether there are 12 mode maps for psolve, as explained below.

-- modified to conform with new runtime system 06/04/90
-- updated with comments for MB, 02/16/94
-- also cleaned up some unnecessary arguments in regime routine
--
-- MB: In order to allow the rhs to be indeterminate, while working
-- mainly with the parametric variables on the lhs (less number of
-- variables), certain conversions of internal representation from
-- GR to Polynomial R and Fraction Polynomial R are done.  At the time
-- of implementation, I thought that this may be more efficient.  I
-- have not done any comparison since that requires rewriting the
-- package.  My own application needs to call this package quite often,
-- and most computations involves only polynomials in the parametric
-- variables.

-- The 12 modes of psolve are probably not all necessary.  Again, I
-- was thinking that if there are many regimes and many ranks, then
-- the output is quite big, and it may be nice to be able to save it
-- and read the results in later to continue computing rather than
-- recomputing.  Because of the combinatorial nature of the algorithm
-- (computing all subdeterminants!), it does not take a very big matrix
-- to get into many regimes.  But I now have second thoughts of this
-- design, since most of the time, the results are just intermediate,
-- passed to be further processed.  On the other hand, there is probably
-- no penalty in leaving the options as is.
\end{verbatim}
\section{package PLEQN ParametricLinearEquations}
<<package PLEQN ParametricLinearEquations>>=
)abbrev package PLEQN ParametricLinearEquations
++  Author: William Sit, spring 89
ParametricLinearEquations(R,Var,Expon,GR):
            Declaration == Definition where

  R         :  Join(EuclideanDomain, CharacteristicZero)
  -- Warning: does not work if R is a field! because of Fraction R
  Var       :  Join(OrderedSet,ConvertibleTo (Symbol))
  Expon     :  OrderedAbelianMonoidSup
  GR        :  PolynomialCategory(R,Expon,Var)
  F        == Fraction R
  FILE ==> FileCategory
  FNAME ==> FileName
  GB  ==> EuclideanGroebnerBasisPackage
--  GBINTERN ==> GroebnerInternalPackage
  I   ==> Integer
  L   ==> List
  M   ==> Matrix
  NNI ==> NonNegativeInteger
  OUT ==> OutputForm
  P   ==> Polynomial
  PI  ==> PositiveInteger
  SEG ==> Segment
  SM  ==> SquareMatrix
  S   ==> String
  V   ==> Vector
  mf    ==> MultivariateFactorize(Var,Expon,R,GR)
  rp    ==> GB(R,Expon,Var,GR)
  gb    ==> GB(R,Expon,Var,GR)
  PR    ==> P R
  GF    ==> Fraction PR
  plift ==> PolynomialCategoryLifting(Expon,Var,R,GR,GF)
  Inputmode ==> Integer
  groebner ==> euclideanGroebner
  redPol  ==> euclideanNormalForm

-- MB: The following macros are data structures to store mostly
-- intermediate results
-- Rec stores a subdeterminant with corresponding row and column indices
-- Fgb is a Groebner basis for the ideal generated by the subdeterminants
--     of a given rank.
-- Linsys specifies a linearly independent system of a given system
--     assuming a given rank, using given row and column indices
-- Linsoln stores the solution to the parametric linear system as a basis
--     and a particular solution (for a given regime)
-- Rec2 stores the rank, and a list of subdeterminants of that rank,
--     and a Groebner basis for the ideal they generate.
-- Rec3 stores a regime and the corresponding solution; the regime is
--     given by a list of equations (eqzro) and one inequation (neqzro)
--     describing the quasi-algebraic set which is the regime; the
--     additional consistency conditions due to the rhs is given by wcond.
-- Ranksolns stores a list of regimes and their solutions, and the number
--     of regimes all together.
-- Rec8 (temporary) stores a quasi-algebraic set with an indication
-- whether it is empty (sysok = false) or not (sysok = true).

-- I think psolve should be renamed parametricSolve, or even
-- parametricLinearSolve.  On the other hand, may be just solve will do.
-- Please feel free to change it to conform with system conventions.
-- Most psolve routines return a list of regimes and solutions,
-- except those that output to file when the number of regimes is
-- returned instead.
-- This version has been tested on the pc version 1.608 March 13, 1992

  Rec   ==> Record(det:GR,rows:L I,cols:L I)
  Eqns  ==> L Rec
  Fgb   ==> L GR  -- groebner basis
  Linsoln   ==> Record(partsol:V GF,basis:L V GF)
  Linsys    ==> Record(mat:M GF,vec:L GF,rank:NNI,rows:L I,cols:L I)
  Rec2  ==> Record(rank:NNI,eqns:Eqns,fgb:Fgb)
  RankConds ==> L Rec2
  Rec3  ==> Record(eqzro:L GR, neqzro:L GR,wcond:L PR, bsoln:Linsoln)
  Ranksolns ==> Record(rgl:L Rec3,rgsz:I)
  Rec8 ==> Record(sysok:Boolean, z0:L GR, n0:L GR)


  Declaration == with
      psolve: (M GR, L GR) -> L Rec3
        ++ psolve(c,w) solves c z = w for all possible ranks
        ++ of the matrix c and given right hand side vector w
        -- this is mode 1
      psolve: (M GR, L Symbol) -> L Rec3
        ++ psolve(c,w) solves c z = w for all possible ranks
        ++ of the matrix c and indeterminate right hand side w
        -- this is mode 2
      psolve:  M GR        -> L Rec3
        ++ psolve(c) solves the homogeneous linear system
        ++ c z = 0 for all possible ranks of the matrix c
        -- this is mode 3
      psolve: (M GR, L GR, PI) -> L Rec3
        ++ psolve(c,w,k) solves c z = w for all possible ranks >= k
        ++ of the matrix c and given right hand side vector w
        -- this is mode 4
      psolve: (M GR, L Symbol, PI) -> L Rec3
        ++ psolve(c,w,k) solves c z = w for all possible ranks >= k
        ++ of the matrix c and indeterminate right hand side w
        -- this is mode 5
      psolve: (M GR,       PI) -> L Rec3
        ++ psolve(c) solves the homogeneous linear system
        ++ c z = 0 for all possible ranks >= k of the matrix c
        -- this is mode 6
      psolve: (M GR, L GR, S) -> I
        ++ psolve(c,w,s) solves c z = w for all possible ranks
        ++ of the matrix c and given right hand side vector w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 7
      psolve: (M GR, L Symbol, S) -> I
        ++ psolve(c,w,s) solves c z = w for all possible ranks
        ++ of the matrix c and indeterminate right hand side w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 8
      psolve: (M GR,       S) -> I
        ++ psolve(c,s) solves c z = 0 for all possible ranks
        ++ of the matrix c and given right hand side vector w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 9
      psolve: (M GR, L GR, PI, S) -> I
        ++ psolve(c,w,k,s) solves c z = w for all possible ranks >= k
        ++ of the matrix c and given right hand side w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 10
      psolve: (M GR, L Symbol, PI, S) -> I
        ++ psolve(c,w,k,s) solves c z = w for all possible ranks >= k
        ++ of the matrix c and indeterminate right hand side w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 11
      psolve: (M GR,       PI, S) -> I
        ++ psolve(c,k,s) solves c z = 0 for all possible ranks >= k
        ++ of the matrix c,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 12

      wrregime   : (L Rec3, S) -> I
        ++ wrregime(l,s) writes a list of regimes to a file named s
        ++ and returns the number of regimes written
      rdregime   : S -> L Rec3
        ++ rdregime(s) reads in a list from a file with name s

        -- for internal use only --
      -- these are exported so my other packages can use them

      bsolve: (M GR, L GF, NNI, S, Inputmode) -> Ranksolns
        ++ bsolve(c, w, r, s, m) returns a list of regimes and
        ++ solutions of the system c z = w for ranks at least r;
        ++ depending on the mode m chosen, it writes the output to
        ++ a file given by the string s.
      dmp2rfi: GR -> GF
        ++ dmp2rfi(p) converts p to target domain
      dmp2rfi: M GR -> M GF
        ++ dmp2rfi(m) converts m to target domain
      dmp2rfi: L GR -> L GF
        ++ dmp2rfi(l) converts l to target domain
      se2rfi:  L Symbol -> L GF
        ++ se2rfi(l) converts l to target domain
      pr2dmp: PR -> GR
        ++ pr2dmp(p) converts p to target domain
      hasoln: (Fgb, L GR) -> Rec8
        ++ hasoln(g, l) tests whether the quasi-algebraic set
        ++ defined by p = 0 for p in g and q ~= 0 for q in l
        ++ is empty or not and returns a simplified definition
        ++ of the quasi-algebraic set
        -- this is now done in QALGSET package
      ParCondList: (M GR,NNI) -> RankConds
        ++ ParCondList(c,r) computes a list of subdeterminants of each
        ++ rank >= r of the matrix c and returns
        ++ a groebner basis for the
        ++ ideal they generate
      redpps: (Linsoln, Fgb) -> Linsoln
        ++ redpps(s,g) returns the simplified form of s after reducing
        ++ modulo a groebner basis g



--               L O C A L    F U N C T I O N S

      B1solve: Linsys -> Linsoln
        ++ B1solve(s) solves the system (s.mat) z = s.vec
        ++ for the variables given by the column indices of s.cols
        ++ in terms of the other variables and the right hand side s.vec
        ++ by assuming that the rank is s.rank,
        ++ that the system is consistent, with the linearly
        ++ independent equations indexed by the given row indices s.rows;
        ++ the coefficients in s.mat involving parameters are treated as
        ++ polynomials.  B1solve(s) returns a particular solution to the
        ++ system and a basis of the homogeneous system (s.mat) z = 0.
      factorset: GR -> L GR
        ++ factorset(p) returns the set of irreducible factors of p.
      maxrank: RankConds -> NNI
        ++ maxrank(r) returns the maximum rank in the list r of regimes
      minrank: RankConds -> NNI
        ++ minrank(r) returns the minimum rank in the list r of regimes
      minset: L L GR -> L L GR
        ++ minset(sl) returns the sublist of sl consisting of the minimal
        ++ lists (with respect to inclusion) in the list sl of lists
      nextSublist: (I, I) -> L L I
        ++ nextSublist(n,k) returns a list of k-subsets of {1, ..., n}.
      overset?: (L GR, L L GR) -> Boolean
        ++ overset?(s,sl) returns true if s properly a sublist of a member
        ++ of sl; otherwise it returns false
      ParCond    : (M GR,NNI) -> Eqns
        ++ ParCond(m,k) returns the list of all k by k subdeterminants in
        ++ the matrix m
      redmat: (M GR, Fgb) -> M GR
        ++ redmat(m,g) returns a matrix whose entries are those of m
        ++ modulo the ideal generated by the groebner basis g
      regime: (Rec,M GR,L GF,L L GR,NNI,NNI,Inputmode) -> Rec3
        ++ regime(y,c, w, p, r, rm, m) returns a regime,
        ++ a list of polynomials specifying the consistency conditions,
        ++ a particular solution and basis representing the general
        ++ solution of the parametric linear system c z = w
        ++ on that regime. The regime returned depends on
        ++ the subdeterminant y.det and the row and column indices.
        ++ The solutions are simplified using the assumption that
        ++ the system has rank r and maximum rank rm. The list p
        ++ represents a list of list of factors of polynomials in
        ++ a groebner basis of the ideal generated by higher order
        ++ subdeterminants, and ius used for the simplification.
        ++ The mode m
        ++ distinguishes the cases when the system is homogeneous,
        ++ or the right hand side is arbitrary, or when there is no
        ++ new right hand side variables.
      sqfree: GR -> GR
        ++ sqfree(p) returns the product of square free factors of p
      inconsistent?: L GR -> Boolean
        ++ inconsistant?(pl) returns true if the system of equations
        ++ p = 0 for p in pl is inconsistent.  It is assumed
        ++ that pl is a groebner basis.
        -- this is needed because of change to
        -- EuclideanGroebnerBasisPackage
      inconsistent?: L PR -> Boolean
        ++ inconsistant?(pl) returns true if the system of equations
        ++ p = 0 for p in pl is inconsistent.  It is assumed
        ++ that pl is a groebner basis.
        -- this is needed because of change to
        -- EuclideanGroebnerBasisPackage

  Definition == add

    inconsistent?(pl:L GR):Boolean ==
      for p in pl repeat
        ground? p => return true
      false
    inconsistent?(pl:L PR):Boolean ==
      for p in pl repeat
        ground? p => return true
      false

    B1solve (sys:Linsys):Linsoln ==
      rss:L I:=sys.rows
      nss:L I:=sys.cols
      k:=sys.rank
      cmat:M GF:=sys.mat
      n:=ncols cmat
      frcols:L I:=setDifference$(L I) (expand$(SEG I) (1..n), nss)
      w:L GF:=sys.vec
      p:V GF:=new(n,0)
      pbas:L V GF:=[]
      if k ~= 0 then
        augmat:M GF:=zero(k,n+1)
        for i in rss for i1 in 1.. repeat
          for j in nss for j1 in 1.. repeat
            augmat(i1,j1):=cmat(i,j)
          for j in frcols for j1 in k+1.. repeat
            augmat(i1,j1):=-cmat(i,j)
          augmat(i1,n+1):=w.i
        augmat:=rowEchelon$(M GF) augmat
        for i in nss for i1 in 1.. repeat p.i:=augmat(i1,n+1)
        for j in frcols for j1 in k+1.. repeat
          pb:V GF:=new(n,0)
          pb.j:=1
          for i in nss for i1 in 1.. repeat
            pb.i:=augmat(i1,j1)
          pbas:=cons(pb,pbas)
      else
        for j in frcols for j1 in k+1.. repeat
          pb:V GF:=new(n,0)
          pb.j:=1
          pbas:=cons(pb,pbas)
      [p,pbas]

    regime (y, coef, w, psbf, rk, rkmax, mode) ==
      -- use the y.det nonzero to simplify the groebner basis
      -- of ideal generated by higher order subdeterminants
      ydetf:L GR:=factorset y.det
      yzero:L GR:=
        rk = rkmax => nil$(L GR)
        psbf:=[setDifference(x, ydetf) for x in psbf]
        groebner$gb [*/x for x in psbf]
      -- simplify coefficients by modulo ideal
      nc:M GF:=dmp2rfi redmat(coef,yzero)
      -- solve the system
      rss:L I:=y.rows;  nss:L I :=y.cols
      sys:Linsys:=[nc,w,rk,rss,nss]$Linsys
      pps:= B1solve(sys)
      pp:=pps.partsol
      frows:L I:=setDifference$(L I) (expand$(SEG I) (1..nrows coef),rss)
      wcd:L PR:= []
      -- case homogeneous rhs
      entry? (mode, [3,6,9,12]$(L I)) =>
               [yzero, ydetf,wcd, redpps(pps, yzero)]$Rec3
      -- case arbitrary rhs, pps not reduced
      for i in frows repeat
          weqn:GF:=+/[nc(i,j)*(pp.j) for j in nss]
          wnum:PR:=numer$GF (w.i - weqn)
          wnum = 0 => "trivially satisfied"
          ground? wnum => return [yzero, ydetf,[1$PR]$(L PR),pps]$Rec3
          wcd:=cons(wnum,wcd)
      entry? (mode, [2,5,8,11]$(L I)) => [yzero, ydetf, wcd, pps]$Rec3
      -- case no new rhs variable
      if not empty? wcd then _
        yzero:=removeDuplicates append(yzero,[pr2dmp pw for pw in wcd])
      test:Rec8:=hasoln (yzero, ydetf)
      not test.sysok => [test.z0, test.n0, [1$PR]$(L PR), pps]$Rec3
      [test.z0, test.n0, [], redpps(pps, test.z0)]$Rec3

    bsolve (coeff, w, h, outname, mode) ==
      r:=nrows coeff
--    n:=ncols coeff
      r ~= #w => error "number of rows unequal on lhs and rhs"
      newfile:FNAME
      rksoln:File Rec3
      count:I:=0
      lrec3:L Rec3:=[]
      filemode:Boolean:= entry? (mode, [7,8,9,10,11,12]$(L I))
      if filemode then
        newfile:=new$FNAME  ("",outname,"regime")
        rksoln:=open$(File Rec3) newfile
      rrcl:RankConds:=
        entry? (mode,[1,2,3,7,8,9]$(L I)) => ParCondList (coeff,0)
        entry? (mode,[4,5,6,10,11,12]$(L I)) => ParCondList (coeff,h)
      rkmax:=maxrank rrcl
      rkmin:=minrank rrcl
      for k in rkmax-rkmin+1..1 by -1 repeat
        rk:=rrcl.k.rank
        pc:Eqns:=rrcl.k.eqns
        psb:Fgb:= (if rk=rkmax then [] else rrcl.(k+1).fgb)
        psbf:L L GR:= [factorset x for x in psb]
        psbf:= minset(psbf)
        for y: Rec in pc repeat
          rec3:Rec3:= regime (y, coeff, w, psbf, rk, rkmax, mode)
          inconsistent? rec3.wcond => "incompatible system"
          if filemode then write!(rksoln, rec3)
          else lrec3:= cons(rec3, lrec3)
          count:=count+1
      if filemode then close! rksoln
      [lrec3, count]$Ranksolns

    factorset y ==
      ground? y => []
      [j.factor for j in factors(factor$mf y)]

    ParCondList (mat, h) ==
      rcl: RankConds:= []
      ps: L GR:=[]
      pc:Eqns:=[]
      npc: Eqns:=[]
      psbf: Fgb:=[]
      done: Boolean := false
      r:=nrows mat
      n:=ncols mat
      maxrk:I:=min(r,n)
      for k: NNI in min(r,n)..h by -1 until done repeat
        pc:= ParCond(mat,k)
        npc:=[]
        (empty? pc) and (k >= 1) => maxrk:= k - 1
        if ground? pc.1.det -- only one is sufficient (neqzro = {})
        then (npc:=pc; done:=true; ps := [1$GR])
        else
          zro:L GR:= (if k = maxrk then [] else rcl.1.fgb)
          covered:Boolean:=false
          for rc: Rec in pc until covered repeat
            p:GR:= redPol$rp (rc.det, zro)
            p = 0 => "incompatible or covered subdeterminant"
            test:=hasoln(zro, [rc.det])
--          zroideal:=ideal(zro)
--          inRadical? (p, zroideal) => "incompatible or covered"
            not test.sysok => "incompatible or covered"
-- The next line is WRONG! cannot replace zro by test.z0
--          zro:=groebner$gb (cons(*/test.n0, test.z0))
            zro:=groebner$gb (cons(p,zro))
            npc:=cons(rc,npc)
            done:= covered:= inconsistent? zro
          ps:=zro
        pcl: Rec2:= construct(k,npc,ps)
        rcl:=cons(pcl,rcl)
      rcl

    redpps(pps, zz) ==
      pv:=pps.partsol
      r:=#pv
      pb:=pps.basis
      n:=#pb + 1
      nummat:M GR:=zero(r,n)
      denmat:M GR:=zero(r,n)
      for i in  1..r repeat
        nummat(i,1):=pr2dmp numer$GF pv.i
        denmat(i,1):=pr2dmp denom$GF pv.i
      for j in 2..n repeat
        for i in 1..r  repeat
          nummat(i,j):=pr2dmp numer$GF (pb.(j-1)).i
          denmat(i,j):=pr2dmp denom$GF (pb.(j-1)).i
      nummat:=redmat(nummat, zz)
      denmat:=redmat(denmat, zz)
      for i in 1..r repeat
        pv.i:=(dmp2rfi nummat(i,1))/(dmp2rfi denmat(i,1))
      for j in 2..n repeat
        pbj:V GF:=new(r,0)
        for i in 1..r repeat
          pbj.i:=(dmp2rfi nummat(i,j))/(dmp2rfi  denmat(i,j))
        pb.(j-1):=pbj
      [pv, pb]

    dmp2rfi (mat:M GR): M GF ==
      r:=nrows mat
      n:=ncols mat
      nmat:M GF:=zero(r,n)
      for i in 1..r repeat
        for j in 1..n repeat
          nmat(i,j):=dmp2rfi mat(i,j)
      nmat

    dmp2rfi (vl: L GR):L GF ==
      [dmp2rfi v for v in vl]

    psolve (mat:M GR, w:L GR): L Rec3 ==
      bsolve(mat, dmp2rfi w, 1, "nofile", 1).rgl
    psolve (mat:M GR, w:L Symbol): L Rec3 ==
      bsolve(mat,  se2rfi w, 1, "nofile", 2).rgl
    psolve (mat:M GR): L Rec3 ==
      bsolve(mat, [0$GF for i in 1..nrows mat], 1, "nofile", 3).rgl

    psolve (mat:M GR, w:L GR, h:PI): L Rec3 ==
      bsolve(mat, dmp2rfi w, h::NNI, "nofile", 4).rgl
    psolve (mat:M GR, w:L Symbol, h:PI): L Rec3 ==
      bsolve(mat, se2rfi w, h::NNI, "nofile", 5).rgl
    psolve (mat:M GR, h:PI): L Rec3 ==
      bsolve(mat, [0$GF for i in 1..nrows mat], h::NNI, "nofile", 6).rgl

    psolve (mat:M GR, w:L GR, outname:S): I ==
      bsolve(mat, dmp2rfi w, 1, outname, 7).rgsz
    psolve (mat:M GR, w:L Symbol, outname:S): I ==
      bsolve(mat, se2rfi w, 1, outname, 8).rgsz
    psolve (mat:M GR, outname:S): I ==
      bsolve(mat, [0$GF for i in 1..nrows mat], 1, outname, 9).rgsz

    nextSublist (n,k) ==
      n <= 0 => []
      k <= 0 => [ nil$(List Integer) ]
      k > n => []
      n = 1 and k = 1 => [[1]]
      mslist: L L I:=[]
      for ms in nextSublist(n-1,k-1) repeat
        mslist:=cons(append(ms,[n]),mslist)
      append(nextSublist(n-1,k), mslist)

    psolve (mat:M GR, w:L GR, h:PI, outname:S): I ==
      bsolve(mat, dmp2rfi w, h::NNI, outname, 10).rgsz
    psolve (mat:M GR, w:L Symbol, h:PI, outname:S): I ==
      bsolve(mat, se2rfi w, h::NNI, outname, 11).rgsz
    psolve (mat:M GR, h:PI, outname:S): I ==
      bsolve(mat,[0$GF for i in 1..nrows mat],h::NNI,outname, 12).rgsz

    hasoln (zro,nzro) ==
      empty? zro => [true, zro, nzro]
      zro:=groebner$gb zro
      inconsistent? zro => [false, zro, nzro]
      empty? nzro =>[true, zro, nzro]
      pnzro:GR:=redPol$rp (*/nzro, zro)
      pnzro = 0 => [false, zro, nzro]
      nzro:=factorset pnzro
      psbf:L L GR:= minset [factorset p for p in zro]
      psbf:= [setDifference(x, nzro) for x in psbf]
      entry? ([], psbf) => [false, zro, nzro]
      zro:=groebner$gb [*/x for x in psbf]
      inconsistent? zro => [false, zro, nzro]
      nzro:=[redPol$rp (p,zro) for p in nzro]
      nzro:=[p for p in nzro | not (ground? p)]
      [true, zro, nzro]



    se2rfi w == [coerce$GF monomial$PR (1$PR, wi, 1) for wi in w]

    pr2dmp p ==
      ground? p => (ground p)::GR
      algCoerceInteractive(p,PR,GR)$(Lisp) pretend GR

    wrregime (lrec3, outname) ==
      newfile:FNAME:=new$FNAME ("",outname,"regime")
      rksoln: File Rec3:=open$(File Rec3) newfile
      count:I:=0  -- number of distinct regimes
--      rec3: Rec3
      for rec3 in lrec3 repeat
          write!(rksoln, rec3)
          count:=count+1
      close!(rksoln)
      count

    dmp2rfi (p:GR):GF ==
      map$plift ((convert #1)@Symbol::GF, #1::PR::GF, p)


    rdregime inname ==
      infilename:=filename$FNAME ("",inname, "regime")
      infile: File Rec3:=open$(File Rec3) (infilename, "input")
      rksoln:L Rec3:=[]
      rec3:Union(Rec3, "failed"):=readIfCan!$(File Rec3) (infile)
      while rec3 case Rec3 repeat
        rksoln:=cons(rec3::Rec3,rksoln) -- replace : to :: for AIX
        rec3:=readIfCan!$(File Rec3) (infile)
      close!(infile)
      rksoln

    maxrank rcl ==
      empty? rcl => 0
      "max"/[j.rank for j in rcl]

    minrank rcl ==
      empty? rcl => 0
      "min"/[j.rank for j in rcl]

    minset lset ==
      empty? lset => lset
      [x for x in lset | not (overset?(x,lset))]

    sqfree p == */[j.factor for j in factors(squareFree p)]


    ParCond (mat, k) ==
      k = 0 => [[1, [], []]$Rec]
      j:NNI:=k::NNI
      DetEqn :Eqns := []
      r:I:= nrows(mat)
      n:I:= ncols(mat)
      k > min(r,n) => error "k exceeds maximum possible rank "
      found:Boolean:=false
      for rss in nextSublist(r, k) until found repeat
        for nss in nextSublist(n, k) until found repeat
          matsub := mat(rss, nss) pretend SM(j, GR)
          detmat := determinant(matsub)
          if detmat ~= 0 then
            found:= (ground? detmat)
            detmat:=sqfree detmat
            neweqn:Rec:=construct(detmat,rss,nss)
            DetEqn:=cons(neweqn, DetEqn)
      found => [first DetEqn]$Eqns
      sort(degree #1.det < degree #2.det, DetEqn)



    overset?(p,qlist) ==
      empty? qlist => false
      or/[part?(brace(q)$Set(GR), brace(p)$Set(GR))$Set(GR) _
                                                for q in qlist]


    redmat (mat,psb) ==
      r:=nrows(mat)
      n:=ncols(mat)
      newmat: M GR:=zero(r,n)
      for i in 1..r repeat
        for j in 1..n repeat
          p:GR:=mat(i,j)
          ground? p => newmat(i,j):=p
          newmat(i,j):=redPol$rp (p,psb)
      newmat

@
<<*>>=
-- See LICENSE.AXIOM for Copyright information

<<package PLEQN ParametricLinearEquations>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}