axiom-source 20170501-3 / usr / share / axiom-20170501 / src / algebra / PLEQN.spad

)abbrev package PLEQN ParametricLinearEquations
++ Author: William Sit, spring 89
++ Description:
++ 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.

ParametricLinearEquations(R,Var,Expon,GR) : SIG == CODE 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
  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)


  SIG ==> 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
      ++ inconsistent?(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
      ++ inconsistent?(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

  CODE ==> 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
      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])
            ^test.sysok => "incompatible or covered"
            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
      for rec3 in lrec3 repeat
          write_!(rksoln, rec3)
          count:=count+1
      close_!(rksoln)
      count

    dmp2rfi (p:GR):GF ==
      map$plift ((v1:Var):GF +-> (convert v1)@Symbol::GF, 
                 (r1:R):GF +-> r1::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((z1:Rec,z2:Rec):Boolean +-> degree z1.det < degree z2.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