Initial population.

1 files changed, 654 insertions, 0 deletions

diff --git a/src/algebra/pleqn.spad.pamphlet b/src/algebra/pleqn.spad.pamphlet
new file mode 100644
index 00000000..8b0569f3
--- /dev/null
+++ b/src/algebra/pleqn.spad.pamphlet
@@ -0,0 +1,654 @@
+\documentclass{article}
+\usepackage{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 ==
+      i,j,i1,j1:I
+      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) ==
+      i,j:I
+      -- 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
+      y:Rec
+      k:NNI
+      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 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:=[]
+      rc: Rec
+      done: Boolean := false
+      r:=nrows mat
+      n:=ncols mat
+      maxrk:I:=min(r,n)
+      k:NNI
+      for k 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 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"
+            ^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 | ^(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 | ^(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/[(brace$(Set GR) q) <$(Set GR) (brace$(Set GR) p) _
+                                                for q in qlist]
+
+
+    redmat (mat,psb) ==
+      i,j:I
+      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}