1)abbrev package PLEQN ParametricLinearEquations 2++ Author: William Sit, spring 89 3++ Description: 4++ This package completely solves a parametric linear system of equations 5++ by decomposing the set of all parametric values for which the linear 6++ system is consistent into a union of quasi-algebraic sets (which need 7++ not be irredundant, but most of the time is). Each quasi-algebraic 8++ set is described by a list of polynomials that vanish on the set, and 9++ a list of polynomials that vanish at no point of the set. 10++ For each quasi-algebraic set, the solution of the linear system 11++ is given, as a particular solution and a basis of the homogeneous 12++ system. 13++ The parametric linear system should be given in matrix form, with 14++ a coefficient matrix and a right hand side vector. The entries 15++ of the coefficient matrix and right hand side vector should be 16++ polynomials in the parametric variables, over a Euclidean domain 17++ of characteristic zero. 18++ 19++ If the system is homogeneous, the right hand side need not be given. 20++ The right hand side can also be replaced by an indeterminate vector, 21++ in which case, the conditions required for consistency will also be 22++ given. 23-- MB: In order to allow the rhs to be indeterminate, while working 24-- mainly with the parametric variables on the lhs (less number of 25-- variables), certain conversions of internal representation from 26-- GR to Polynomial R and Fraction Polynomial R are done. At the time 27-- of implementation, I thought that this may be more efficient. I 28-- have not done any comparison since that requires rewriting the 29-- package. My own application needs to call this package quite often, 30-- and most computations involves only polynomials in the parametric 31-- variables. 32 33-- The 12 modes of psolve are probably not all necessary. Again, I 34-- was thinking that if there are many regimes and many ranks, then 35-- the output is quite big, and it may be nice to be able to save it 36-- and read the results in later to continue computing rather than 37-- recomputing. Because of the combinatorial nature of the algorithm 38-- (computing all subdeterminants!), it does not take a very big matrix 39-- to get into many regimes. But I now have second thoughts of this 40-- design, since most of the time, the results are just intermediate, 41-- passed to be further processed. On the other hand, there is probably 42-- no penalty in leaving the options as is. 43ParametricLinearEquations(R, Var, Expon, GR): 44 Declaration == Definition where 45 46 R : Join(EuclideanDomain, PolynomialFactorizationExplicit, 47 CharacteristicZero) 48 -- Warning: does not work if R is a field! because of Fraction R 49 Var : Join(OrderedSet, ConvertibleTo (Symbol)) 50 Expon : OrderedAbelianMonoidSup 51 GR : PolynomialCategory(R, Expon, Var) 52 F == Fraction R 53 FILE ==> FileCategory 54 FNAME ==> FileName 55 GB ==> EuclideanGroebnerBasisPackage 56-- GBINTERN ==> GroebnerInternalPackage 57 I ==> Integer 58 L ==> List 59 M ==> Matrix 60 NNI ==> NonNegativeInteger 61 OUT ==> OutputForm 62 P ==> Polynomial 63 PI ==> PositiveInteger 64 SEG ==> Segment 65 SM ==> SquareMatrix 66 S ==> String 67 V ==> Vector 68 mf ==> MultivariateFactorize(Var, Expon, R, GR) 69 rp ==> GB(R, Expon, Var, GR) 70 gb ==> GB(R, Expon, Var, GR) 71 PR ==> P R 72 GF ==> Fraction PR 73 plift ==> PolynomialCategoryLifting(Expon, Var, R, GR, GF) 74 Inputmode ==> Integer 75 groebner ==> euclideanGroebner 76 redPol ==> euclideanNormalForm 77 78-- MB: The following macros are data structures to store mostly 79-- intermediate results 80-- Rec stores a subdeterminant with corresponding row and column indices 81-- Fgb is a Groebner basis for the ideal generated by the subdeterminants 82-- of a given rank. 83-- Linsys specifies a linearly independent system of a given system 84-- assuming a given rank, using given row and column indices 85-- Linsoln stores the solution to the parametric linear system as a basis 86-- and a particular solution (for a given regime) 87-- Rec2 stores the rank, and a list of subdeterminants of that rank, 88-- and a Groebner basis for the ideal they generate. 89-- Rec3 stores a regime and the corresponding solution; the regime is 90-- given by a list of equations (eqzro) and one inequation (neqzro) 91-- describing the quasi-algebraic set which is the regime; the 92-- additional consistency conditions due to the rhs is given by wcond. 93-- Ranksolns stores a list of regimes and their solutions, and the number 94-- of regimes all together. 95-- Rec8 (temporary) stores a quasi-algebraic set with an indication 96-- whether it is empty (sysok = false) or not (sysok = true). 97 98-- I think psolve should be renamed parametricSolve, or even 99-- parametricLinearSolve. On the other hand, may be just solve will do. 100-- Please feel free to change it to conform with system conventions. 101-- Most psolve routines return a list of regimes and solutions, 102-- except those that output to file when the number of regimes is 103-- returned instead. 104-- This version has been tested on the pc version 1.608 March 13, 1992 105 106 Rec ==> Record(det : GR, rows : L I, cols : L I) 107 Eqns ==> L Rec 108 Fgb ==> L GR -- groebner basis 109 Linsoln ==> Record(partsol : V GF, basis : L V GF) 110 Linsys ==> Record(mat : M GF, vec : L GF, rank : NNI, rows : L I, cols : L I) 111 Rec2 ==> Record(rank : NNI, eqns : Eqns, fgb : Fgb) 112 RankConds ==> L Rec2 113 Rec3 ==> Record(eqzro : L GR, neqzro : L GR, wcond : L PR, bsoln : Linsoln) 114 Ranksolns ==> Record(rgl : L Rec3, rgsz : I) 115 Rec8 ==> Record(sysok : Boolean, z0 : L GR, n0 : L GR) 116 117 118 Declaration == with 119 psolve : (M GR, L GR) -> L Rec3 120 ++ psolve(c, w) solves c z = w for all possible ranks 121 ++ of the matrix c and given right hand side vector w 122 -- this is mode 1 123 psolve : (M GR, L Symbol) -> L Rec3 124 ++ psolve(c, w) solves c z = w for all possible ranks 125 ++ of the matrix c and indeterminate right hand side w 126 -- this is mode 2 127 psolve : M GR -> L Rec3 128 ++ psolve(c) solves the homogeneous linear system 129 ++ c z = 0 for all possible ranks of the matrix c 130 -- this is mode 3 131 psolve : (M GR, L GR, PI) -> L Rec3 132 ++ psolve(c, w, k) solves c z = w for all possible ranks >= k 133 ++ of the matrix c and given right hand side vector w 134 -- this is mode 4 135 psolve : (M GR, L Symbol, PI) -> L Rec3 136 ++ psolve(c, w, k) solves c z = w for all possible ranks >= k 137 ++ of the matrix c and indeterminate right hand side w 138 -- this is mode 5 139 psolve : (M GR, PI) -> L Rec3 140 ++ psolve(c) solves the homogeneous linear system 141 ++ c z = 0 for all possible ranks >= k of the matrix c 142 -- this is mode 6 143 psolve : (M GR, L GR, S) -> I 144 ++ psolve(c, w, s) solves c z = w for all possible ranks 145 ++ of the matrix c and given right hand side vector w, 146 ++ writes the results to a file named s, and returns the 147 ++ number of regimes 148 -- this is mode 7 149 psolve : (M GR, L Symbol, S) -> I 150 ++ psolve(c, w, s) solves c z = w for all possible ranks 151 ++ of the matrix c and indeterminate right hand side w, 152 ++ writes the results to a file named s, and returns the 153 ++ number of regimes 154 -- this is mode 8 155 psolve : (M GR, S) -> I 156 ++ psolve(c, s) solves c z = 0 for all possible ranks 157 ++ of the matrix c and given right hand side vector w, 158 ++ writes the results to a file named s, and returns the 159 ++ number of regimes 160 -- this is mode 9 161 psolve : (M GR, L GR, PI, S) -> I 162 ++ psolve(c, w, k, s) solves c z = w for all possible ranks >= k 163 ++ of the matrix c and given right hand side w, 164 ++ writes the results to a file named s, and returns the 165 ++ number of regimes 166 -- this is mode 10 167 psolve : (M GR, L Symbol, PI, S) -> I 168 ++ psolve(c, w, k, s) solves c z = w for all possible ranks >= k 169 ++ of the matrix c and indeterminate right hand side w, 170 ++ writes the results to a file named s, and returns the 171 ++ number of regimes 172 -- this is mode 11 173 psolve : (M GR, PI, S) -> I 174 ++ psolve(c, k, s) solves c z = 0 for all possible ranks >= k 175 ++ of the matrix c, 176 ++ writes the results to a file named s, and returns the 177 ++ number of regimes 178 -- this is mode 12 179 180 wrregime : (L Rec3, S) -> I 181 ++ wrregime(l, s) writes a list of regimes to a file named s 182 ++ and returns the number of regimes written 183 rdregime : S -> L Rec3 184 ++ rdregime(s) reads in a list from a file with name s 185 186 -- for internal use only -- 187 -- these are exported so my other packages can use them 188 189 bsolve : (M GR, L GF, NNI, S, Inputmode) -> Ranksolns 190 ++ bsolve(c, w, r, s, m) returns a list of regimes and 191 ++ solutions of the system c z = w for ranks at least r; 192 ++ depending on the mode m chosen, it writes the output to 193 ++ a file given by the string s. 194 dmp2rfi : GR -> GF 195 ++ dmp2rfi(p) converts p to target domain 196 dmp2rfi : M GR -> M GF 197 ++ dmp2rfi(m) converts m to target domain 198 dmp2rfi : L GR -> L GF 199 ++ dmp2rfi(l) converts l to target domain 200 se2rfi : L Symbol -> L GF 201 ++ se2rfi(l) converts l to target domain 202 pr2dmp : PR -> GR 203 ++ pr2dmp(p) converts p to target domain 204 hasoln : (Fgb, L GR) -> Rec8 205 ++ hasoln(g, l) tests whether the quasi-algebraic set 206 ++ defined by p = 0 for p in g and q ~= 0 for q in l 207 ++ is empty or not and returns a simplified definition 208 ++ of the quasi-algebraic set 209 -- this is now done in QALGSET package 210 ParCondList : (M GR, NNI) -> RankConds 211 ++ ParCondList(c, r) computes a list of subdeterminants of each 212 ++ rank >= r of the matrix c and returns 213 ++ a groebner basis for the 214 ++ ideal they generate 215 redpps : (Linsoln, Fgb) -> Linsoln 216 ++ redpps(s, g) returns the simplified form of s after reducing 217 ++ modulo a groebner basis g 218 219 220 221-- L O C A L F U N C T I O N S 222 223 B1solve : Linsys -> Linsoln 224 ++ B1solve(s) solves the system (s.mat) z = s.vec 225 ++ for the variables given by the column indices of s.cols 226 ++ in terms of the other variables and the right hand side s.vec 227 ++ by assuming that the rank is s.rank, 228 ++ that the system is consistent, with the linearly 229 ++ independent equations indexed by the given row indices s.rows; 230 ++ the coefficients in s.mat involving parameters are treated as 231 ++ polynomials. B1solve(s) returns a particular solution to the 232 ++ system and a basis of the homogeneous system (s.mat) z = 0. 233 factorset : GR -> L GR 234 ++ factorset(p) returns the set of irreducible factors of p. 235 maxrank : RankConds -> NNI 236 ++ maxrank(r) returns the maximum rank in the list r of regimes 237 minrank : RankConds -> NNI 238 ++ minrank(r) returns the minimum rank in the list r of regimes 239 minset : L L GR -> L L GR 240 ++ minset(sl) returns the sublist of sl consisting of the minimal 241 ++ lists (with respect to inclusion) in the list sl of lists 242 nextSublist : (I, I) -> L L I 243 ++ nextSublist(n, k) returns a list of k-subsets of {1, ..., n}. 244 overset? : (L GR, L L GR) -> Boolean 245 ++ overset?(s, sl) returns true if s properly a sublist of a member 246 ++ of sl; otherwise it returns false 247 ParCond : (M GR, NNI) -> Eqns 248 ++ ParCond(m, k) returns the list of all k by k subdeterminants in 249 ++ the matrix m 250 redmat : (M GR, Fgb) -> M GR 251 ++ redmat(m, g) returns a matrix whose entries are those of m 252 ++ modulo the ideal generated by the groebner basis g 253 regime : (Rec, M GR, L GF, L L GR, NNI, NNI, Inputmode) -> Rec3 254 ++ regime(y, c, w, p, r, rm, m) returns a regime, 255 ++ a list of polynomials specifying the consistency conditions, 256 ++ a particular solution and basis representing the general 257 ++ solution of the parametric linear system c z = w 258 ++ on that regime. The regime returned depends on 259 ++ the subdeterminant y.det and the row and column indices. 260 ++ The solutions are simplified using the assumption that 261 ++ the system has rank r and maximum rank rm. The list p 262 ++ represents a list of list of factors of polynomials in 263 ++ a groebner basis of the ideal generated by higher order 264 ++ subdeterminants, and ius used for the simplification. 265 ++ The mode m 266 ++ distinguishes the cases when the system is homogeneous, 267 ++ or the right hand side is arbitrary, or when there is no 268 ++ new right hand side variables. 269 sqfree : GR -> GR 270 ++ sqfree(p) returns the product of square free factors of p 271 inconsistent? : L GR -> Boolean 272 ++ inconsistent?(pl) returns true if the system of equations 273 ++ p = 0 for p in pl is inconsistent. It is assumed 274 ++ that pl is a groebner basis. 275 -- this is needed because of change to 276 -- EuclideanGroebnerBasisPackage 277 inconsistent? : L PR -> Boolean 278 ++ inconsistent?(pl) returns true if the system of equations 279 ++ p = 0 for p in pl is inconsistent. It is assumed 280 ++ that pl is a groebner basis. 281 -- this is needed because of change to 282 -- EuclideanGroebnerBasisPackage 283 284 Definition == add 285 286 inconsistent?(pl : L GR) : Boolean == 287 for p in pl repeat 288 ground? p => return true 289 false 290 inconsistent?(pl : L PR) : Boolean == 291 for p in pl repeat 292 ground? p => return true 293 false 294 295 B1solve (sys : Linsys) : Linsoln == 296 i, j, i1, j1 : I 297 rss : L I := sys.rows 298 nss : L I := sys.cols 299 k := sys.rank 300 cmat : M GF := sys.mat 301 n := ncols cmat 302 frcols : L I := setDifference$(L I) (expand$(SEG I) (1..n), nss) 303 w : L GF := sys.vec 304 p : V GF := new(n, 0) 305 pbas : L V GF := [] 306 if k ~= 0 then 307 augmat : M GF := zero(k, n+1) 308 for i in rss for i1 in 1.. repeat 309 for j in nss for j1 in 1.. repeat 310 augmat(i1, j1) := cmat(i, j) 311 for j in frcols for j1 in k+1.. repeat 312 augmat(i1, j1) := -cmat(i, j) 313 augmat(i1, n+1) := w.i 314 augmat := rowEchelon$(M GF) augmat 315 for i in nss for i1 in 1.. repeat p.i := augmat(i1, n+1) 316 for j in frcols for j1 in k+1.. repeat 317 pb : V GF := new(n, 0) 318 pb.j := 1 319 for i in nss for i1 in 1.. repeat 320 pb.i := augmat(i1, j1) 321 pbas := cons(pb, pbas) 322 else 323 for j in frcols for j1 in k+1.. repeat 324 pb : V GF := new(n, 0) 325 pb.j := 1 326 pbas := cons(pb, pbas) 327 [p, pbas] 328 329 regime (y, coef, w, psbf, rk, rkmax, mode) == 330 i, j : I 331 -- use the y.det nonzero to simplify the groebner basis 332 -- of ideal generated by higher order subdeterminants 333 ydetf : L GR := factorset y.det 334 yzero : L GR := 335 rk = rkmax => []$(L GR) 336 psbf := [setDifference(x, ydetf) for x in psbf] 337 groebner$gb [*/x for x in psbf] 338 -- simplify coefficients by modulo ideal 339 nc : M GF := dmp2rfi redmat(coef, yzero) 340 -- solve the system 341 rss : L I := y.rows; nss : L I := y.cols 342 sys : Linsys := [nc, w, rk, rss, nss]$Linsys 343 pps := B1solve(sys) 344 pp := pps.partsol 345 frows : L I := setDifference$(L I) (expand$(SEG I) (1..nrows coef), rss) 346 wcd : L PR := [] 347 -- case homogeneous rhs 348 entry? (mode, [3, 6, 9, 12]$(L I)) => 349 [yzero, ydetf, wcd, redpps(pps, yzero)]$Rec3 350 -- case arbitrary rhs, pps not reduced 351 for i in frows repeat 352 weqn : GF := +/[nc(i, j)*(pp.j) for j in nss] 353 wnum : PR := numer$GF (w.i - weqn) 354 wnum = 0 => "trivially satisfied" 355 ground? wnum => return [yzero, ydetf, [1$PR]$(L PR), pps]$Rec3 356 wcd := cons(wnum, wcd) 357 entry? (mode, [2, 5, 8, 11]$(L I)) => [yzero, ydetf, wcd, pps]$Rec3 358 -- case no new rhs variable 359 if not empty? wcd then _ 360 yzero := removeDuplicates append(yzero, [pr2dmp pw for pw in wcd]) 361 test : Rec8 := hasoln (yzero, ydetf) 362 not test.sysok => [test.z0, test.n0, [1$PR]$(L PR), pps]$Rec3 363 [test.z0, test.n0, [], redpps(pps, test.z0)]$Rec3 364 365 bsolve (coeff, w, h, outname, mode) == 366 r := nrows coeff 367-- n := ncols coeff 368 r ~= #w => error "number of rows unequal on lhs and rhs" 369 newfile : FNAME 370 rksoln : File Rec3 371 count : I := 0 372 lrec3 : L Rec3 := [] 373 filemode : Boolean := entry? (mode, [7, 8, 9, 10, 11, 12]$(L I)) 374 if filemode then 375 newfile := new$FNAME ("",outname,"regime") 376 rksoln := open$(File Rec3) newfile 377 y : Rec 378 k : NNI 379 rrcl : RankConds := 380 entry? (mode, [1, 2, 3, 7, 8, 9]$(L I)) => ParCondList (coeff, 0) 381 entry? (mode, [4, 5, 6, 10, 11, 12]$(L I)) => ParCondList (coeff, h) 382 rkmax := maxrank rrcl 383 rkmin := minrank rrcl 384 for k in rkmax-rkmin+1..1 by -1 repeat 385 rk := rrcl.k.rank 386 pc : Eqns := rrcl.k.eqns 387 psb : Fgb := (if rk = rkmax then [] else rrcl.(k+1).fgb) 388 psbf : L L GR := [factorset x for x in psb] 389 psbf := minset(psbf) 390 for y in pc repeat 391 rec3 : Rec3 := regime (y, coeff, w, psbf, rk, rkmax, mode) 392 inconsistent? rec3.wcond => "incompatible system" 393 if filemode then write!(rksoln, rec3) 394 else lrec3 := cons(rec3, lrec3) 395 count := count+1 396 if filemode then close! rksoln 397 [lrec3, count]$Ranksolns 398 399 factorset y == 400 ground? y => [] 401 [j.factor for j in factorList(factor$mf y)] 402 403 ParCondList (mat, h) == 404 rcl : RankConds := [] 405 ps : L GR := [] 406 pc : Eqns := [] 407 npc : Eqns := [] 408 psbf : Fgb := [] 409 rc : Rec 410 done : Boolean := false 411 r := nrows mat 412 n := ncols mat 413 maxrk : I := min(r, n) 414 k : NNI 415 for k in min(r, n)..h by -1 until done repeat 416 pc := ParCond(mat, k) 417 npc := [] 418 (empty? pc) and (k >= 1) => maxrk := k - 1 419 if ground? pc.1.det -- only one is sufficient (neqzro = {}) 420 then (npc := pc; done := true; ps := [1$GR]) 421 else 422 zro : L GR := (if k = maxrk then [] else rcl.1.fgb) 423 covered : Boolean := false 424 for rc in pc until covered repeat 425 p : GR := redPol$rp (rc.det, zro) 426 p = 0 => "incompatible or covered subdeterminant" 427 test := hasoln(zro, [rc.det]) 428-- zroideal := ideal(zro) 429-- inRadical? (p, zroideal) => "incompatible or covered" 430 not test.sysok => "incompatible or covered" 431-- The next line is WRONG! cannot replace zro by test.z0 432-- zro := groebner$gb (cons(*/test.n0, test.z0)) 433 zro := groebner$gb (cons(p, zro)) 434 npc := cons(rc, npc) 435 done := covered := inconsistent? zro 436 ps := zro 437 pcl : Rec2 := construct(k, npc, ps) 438 rcl := cons(pcl, rcl) 439 rcl 440 441 redpps(pps, zz) == 442 pv := pps.partsol 443 r := #pv 444 pb := pps.basis 445 n := #pb + 1 446 nummat : M GR := zero(r, n) 447 denmat : M GR := zero(r, n) 448 for i in 1..r repeat 449 nummat(i, 1) := pr2dmp numer$GF pv.i 450 denmat(i, 1) := pr2dmp denom$GF pv.i 451 for j in 2..n repeat 452 for i in 1..r repeat 453 nummat(i, j) := pr2dmp numer$GF (pb.(j-1)).i 454 denmat(i, j) := pr2dmp denom$GF (pb.(j-1)).i 455 nummat := redmat(nummat, zz) 456 denmat := redmat(denmat, zz) 457 for i in 1..r repeat 458 pv.i := (dmp2rfi nummat(i, 1))/(dmp2rfi denmat(i, 1)) 459 for j in 2..n repeat 460 pbj : V GF := new(r, 0) 461 for i in 1..r repeat 462 pbj.i := (dmp2rfi nummat(i, j))/(dmp2rfi denmat(i, j)) 463 pb.(j-1) := pbj 464 [pv, pb] 465 466 dmp2rfi (mat : M GR) : M GF == 467 r := nrows mat 468 n := ncols mat 469 nmat : M GF := zero(r, n) 470 for i in 1..r repeat 471 for j in 1..n repeat 472 nmat(i, j) := dmp2rfi mat(i, j) 473 nmat 474 475 dmp2rfi (vl : L GR) : L GF == 476 [dmp2rfi v for v in vl] 477 478 psolve (mat : M GR, w : L GR) : L Rec3 == 479 bsolve(mat, dmp2rfi w, 1, "nofile", 1).rgl 480 psolve (mat : M GR, w : L Symbol) : L Rec3 == 481 bsolve(mat, se2rfi w, 1, "nofile", 2).rgl 482 psolve (mat : M GR) : L Rec3 == 483 bsolve(mat, [0$GF for i in 1..nrows mat], 1, "nofile", 3).rgl 484 485 psolve (mat : M GR, w : L GR, h : PI) : L Rec3 == 486 bsolve(mat, dmp2rfi w, h::NNI, "nofile", 4).rgl 487 psolve (mat : M GR, w : L Symbol, h : PI) : L Rec3 == 488 bsolve(mat, se2rfi w, h::NNI, "nofile", 5).rgl 489 psolve (mat : M GR, h : PI) : L Rec3 == 490 bsolve(mat, [0$GF for i in 1..nrows mat], h::NNI, "nofile", 6).rgl 491 492 psolve (mat : M GR, w : L GR, outname : S) : I == 493 bsolve(mat, dmp2rfi w, 1, outname, 7).rgsz 494 psolve (mat : M GR, w : L Symbol, outname : S) : I == 495 bsolve(mat, se2rfi w, 1, outname, 8).rgsz 496 psolve (mat : M GR, outname : S) : I == 497 bsolve(mat, [0$GF for i in 1..nrows mat], 1, outname, 9).rgsz 498 499 nextSublist (n, k) == 500 n <= 0 => [] 501 k <= 0 => [ []$(List Integer) ] 502 k > n => [] 503 n = 1 and k = 1 => [[1]] 504 mslist : L L I := [] 505 for ms in nextSublist(n-1, k-1) repeat 506 mslist := cons(append(ms, [n]), mslist) 507 append(nextSublist(n-1, k), mslist) 508 509 psolve (mat : M GR, w : L GR, h : PI, outname : S) : I == 510 bsolve(mat, dmp2rfi w, h::NNI, outname, 10).rgsz 511 psolve (mat : M GR, w : L Symbol, h : PI, outname : S) : I == 512 bsolve(mat, se2rfi w, h::NNI, outname, 11).rgsz 513 psolve (mat : M GR, h : PI, outname : S) : I == 514 bsolve(mat, [0$GF for i in 1..nrows mat], h::NNI, outname, 12).rgsz 515 516 hasoln (zro, nzro) == 517 empty? zro => [true, zro, nzro] 518 zro := groebner$gb zro 519 inconsistent? zro => [false, zro, nzro] 520 empty? nzro =>[true, zro, nzro] 521 pnzro : GR := redPol$rp (*/nzro, zro) 522 pnzro = 0 => [false, zro, nzro] 523 nzro := factorset pnzro 524 psbf : L L GR := minset [factorset p for p in zro] 525 psbf := [setDifference(x, nzro) for x in psbf] 526 entry? ([], psbf) => [false, zro, nzro] 527 zro := groebner$gb [*/x for x in psbf] 528 inconsistent? zro => [false, zro, nzro] 529 nzro := [redPol$rp (p, zro) for p in nzro] 530 nzro := [p for p in nzro | not (ground? p)] 531 [true, zro, nzro] 532 533 534 535 se2rfi w == [coerce$GF monomial$PR (1$PR, wi, 1) for wi in w] 536 537 pr2dmp p == 538 ground? p => (ground p)::GR 539 algCoerceInteractive(p, PR, GR)$(Lisp) pretend GR 540 541 wrregime (lrec3, outname) == 542 newfile:FNAME := new$FNAME ("",outname,"regime") 543 rksoln : File Rec3 := open$(File Rec3) newfile 544 count : I := 0 -- number of distinct regimes 545-- rec3: Rec3 546 for rec3 in lrec3 repeat 547 write!(rksoln, rec3) 548 count := count+1 549 close!(rksoln) 550 count 551 552 dmp2rfi (p : GR) : GF == 553 map$plift ((v1 : Var) : GF +-> (convert v1)@Symbol::GF, _ 554 (r1 : R) : GF +-> r1::PR::GF, p) 555 556 557 rdregime inname == 558 infilename := filename$FNAME ("",inname, "regime") 559 infile : File Rec3 := open$(File Rec3) (infilename, "input") 560 rksoln : L Rec3 := [] 561 rec3 : Union(Rec3, "failed") := readIfCan!$(File Rec3) (infile) 562 while rec3 case Rec3 repeat 563 rksoln := cons(rec3::Rec3, rksoln) -- replace : to :: for AIX 564 rec3 := readIfCan!$(File Rec3) (infile) 565 close!(infile) 566 rksoln 567 568 maxrank rcl == 569 empty? rcl => 0 570 "max"/[j.rank for j in rcl] 571 572 minrank rcl == 573 empty? rcl => 0 574 "min"/[j.rank for j in rcl] 575 576 minset lset == 577 empty? lset => lset 578 [x for x in lset | not (overset?(x, lset))] 579 580 sqfree p == */[j.factor for j in factorList(squareFree p)] 581 582 583 ParCond (mat, k) == 584 k = 0 => [[1, [], []]$Rec] 585 j : NNI := k::NNI 586 DetEqn : Eqns := [] 587 r : I := nrows(mat) 588 n : I := ncols(mat) 589 k > min(r,n) => error "k exceeds maximum possible rank " 590 found : Boolean := false 591 for rss in nextSublist(r, k) until found repeat 592 for nss in nextSublist(n, k) until found repeat 593 matsub := mat(rss, nss) pretend SM(j, GR) 594 detmat := determinant(matsub) 595 if detmat ~= 0 then 596 found := (ground? detmat) 597 detmat := sqfree detmat 598 neweqn : Rec := construct(detmat, rss, nss) 599 DetEqn := cons(neweqn, DetEqn) 600 found => [first DetEqn]$Eqns 601 sort((z1 : Rec, z2 : Rec) : Boolean +-> degree z1.det < degree z2.det, DetEqn) 602 603 604 605 overset?(p, qlist) == 606 empty? qlist => false 607 or/[(set$(Set GR) q) <$(Set GR) (set$(Set GR) p) _ 608 for q in qlist] 609 610 611 redmat (mat, psb) == 612 i, j : I 613 r := nrows(mat) 614 n := ncols(mat) 615 newmat : M GR := zero(r, n) 616 for i in 1..r repeat 617 for j in 1..n repeat 618 p : GR := mat(i, j) 619 ground? p => newmat(i, j) := p 620 newmat(i, j) := redPol$rp (p, psb) 621 newmat 622 623-- See LICENSE.AXIOM for Copyright information 624