\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} <>= )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 <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}