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.
-- 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.
ParametricLinearEquations(R, Var, Expon, GR):
            Declaration == Definition where

  R         :  Join(EuclideanDomain, PolynomialFactorizationExplicit,
                    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
        ++ 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

  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 => []$(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 factorList(factor$mf y)]

    ParCondList (mat, h) ==
      rcl : RankConds := []
      ps : L GR := []
      pc : Eqns := []
      npc : Eqns := []
      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"
            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 => [ []$(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 ((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 | not (overset?(x, lset))]

    sqfree p == */[j.factor for j in factorList(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/[(set$(Set GR) q) <$(Set GR) (set$(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.txt for Copyright information