aboutsummaryrefslogtreecommitdiff
path: root/src/algebra/pleqn.spad.pamphlet
blob: f5653815373ec609e586c00bad8795839bf02d59 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
\documentclass{article}
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra pleqn.spad}
\author{William Sit}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
\begin{verbatim}
++ This package completely solves a parametric linear system of equations
++ by decomposing the set of all parametric values for which the linear
++ system is consistent into a union of quasi-algebraic  sets (which need
++ not be irredundant, but most of the time is). Each quasi-algebraic
++ set is described by a list of polynomials that vanish on the set, and
++ a list of polynomials that vanish at no point of the set.
++ For each quasi-algebraic set, the solution of the linear system
++ is given, as a particular solution and  a basis of the homogeneous
++ system.
++ The parametric linear system should be given in matrix form, with
++ a coefficient matrix and a right hand side vector. The entries
++ of the coefficient matrix and right hand side vector should be
++ polynomials in the parametric variables, over a Euclidean domain
++ of characteristic zero.
++
++ If the system is homogeneous, the right hand side need not be given.
++ The right hand side can also be  replaced by an indeterminate vector,
++ in which case, the conditions required for consistency will also be
++ given.
++
++ The package has other facilities for saving results to external
++ files, as well as solving the system for a specified minimum rank.
++ Altogether there are 12 mode maps for psolve, as explained below.

-- modified to conform with new runtime system 06/04/90
-- updated with comments for MB, 02/16/94
-- also cleaned up some unnecessary arguments in regime routine
--
-- MB: In order to allow the rhs to be indeterminate, while working
-- mainly with the parametric variables on the lhs (less number of
-- variables), certain conversions of internal representation from
-- GR to Polynomial R and Fraction Polynomial R are done.  At the time
-- of implementation, I thought that this may be more efficient.  I
-- have not done any comparison since that requires rewriting the
-- package.  My own application needs to call this package quite often,
-- and most computations involves only polynomials in the parametric
-- variables.

-- The 12 modes of psolve are probably not all necessary.  Again, I
-- was thinking that if there are many regimes and many ranks, then
-- the output is quite big, and it may be nice to be able to save it
-- and read the results in later to continue computing rather than
-- recomputing.  Because of the combinatorial nature of the algorithm
-- (computing all subdeterminants!), it does not take a very big matrix
-- to get into many regimes.  But I now have second thoughts of this
-- design, since most of the time, the results are just intermediate,
-- passed to be further processed.  On the other hand, there is probably
-- no penalty in leaving the options as is.
\end{verbatim}
\section{package PLEQN ParametricLinearEquations}
<<package PLEQN ParametricLinearEquations>>=
)abbrev package PLEQN ParametricLinearEquations
++  Author: William Sit, spring 89
ParametricLinearEquations(R,Var,Expon,GR):
            Declaration == Definition where

  R         :  Join(EuclideanDomain, CharacteristicZero)
  -- Warning: does not work if R is a field! because of Fraction R
  Var       :  Join(OrderedSet,ConvertibleTo (Symbol))
  Expon     :  OrderedAbelianMonoidSup
  GR        :  PolynomialCategory(R,Expon,Var)
  F        == Fraction R
  FILE ==> FileCategory
  FNAME ==> FileName
  GB  ==> EuclideanGroebnerBasisPackage
--  GBINTERN ==> GroebnerInternalPackage
  I   ==> Integer
  L   ==> List
  M   ==> Matrix
  NNI ==> NonNegativeInteger
  OUT ==> OutputForm
  P   ==> Polynomial
  PI  ==> PositiveInteger
  SEG ==> Segment
  SM  ==> SquareMatrix
  S   ==> String
  V   ==> Vector
  mf    ==> MultivariateFactorize(Var,Expon,R,GR)
  rp    ==> GB(R,Expon,Var,GR)
  gb    ==> GB(R,Expon,Var,GR)
  PR    ==> P R
  GF    ==> Fraction PR
  plift ==> PolynomialCategoryLifting(Expon,Var,R,GR,GF)
  Inputmode ==> Integer
  groebner ==> euclideanGroebner
  redPol  ==> euclideanNormalForm

-- MB: The following macros are data structures to store mostly
-- intermediate results
-- Rec stores a subdeterminant with corresponding row and column indices
-- Fgb is a Groebner basis for the ideal generated by the subdeterminants
--     of a given rank.
-- Linsys specifies a linearly independent system of a given system
--     assuming a given rank, using given row and column indices
-- Linsoln stores the solution to the parametric linear system as a basis
--     and a particular solution (for a given regime)
-- Rec2 stores the rank, and a list of subdeterminants of that rank,
--     and a Groebner basis for the ideal they generate.
-- Rec3 stores a regime and the corresponding solution; the regime is
--     given by a list of equations (eqzro) and one inequation (neqzro)
--     describing the quasi-algebraic set which is the regime; the
--     additional consistency conditions due to the rhs is given by wcond.
-- Ranksolns stores a list of regimes and their solutions, and the number
--     of regimes all together.
-- Rec8 (temporary) stores a quasi-algebraic set with an indication
-- whether it is empty (sysok = false) or not (sysok = true).

-- I think psolve should be renamed parametricSolve, or even
-- parametricLinearSolve.  On the other hand, may be just solve will do.
-- Please feel free to change it to conform with system conventions.
-- Most psolve routines return a list of regimes and solutions,
-- except those that output to file when the number of regimes is
-- returned instead.
-- This version has been tested on the pc version 1.608 March 13, 1992

  Rec   ==> Record(det:GR,rows:L I,cols:L I)
  Eqns  ==> L Rec
  Fgb   ==> L GR  -- groebner basis
  Linsoln   ==> Record(partsol:V GF,basis:L V GF)
  Linsys    ==> Record(mat:M GF,vec:L GF,rank:NNI,rows:L I,cols:L I)
  Rec2  ==> Record(rank:NNI,eqns:Eqns,fgb:Fgb)
  RankConds ==> L Rec2
  Rec3  ==> Record(eqzro:L GR, neqzro:L GR,wcond:L PR, bsoln:Linsoln)
  Ranksolns ==> Record(rgl:L Rec3,rgsz:I)
  Rec8 ==> Record(sysok:Boolean, z0:L GR, n0:L GR)


  Declaration == with
      psolve: (M GR, L GR) -> L Rec3
        ++ psolve(c,w) solves c z = w for all possible ranks
        ++ of the matrix c and given right hand side vector w
        -- this is mode 1
      psolve: (M GR, L Symbol) -> L Rec3
        ++ psolve(c,w) solves c z = w for all possible ranks
        ++ of the matrix c and indeterminate right hand side w
        -- this is mode 2
      psolve:  M GR        -> L Rec3
        ++ psolve(c) solves the homogeneous linear system
        ++ c z = 0 for all possible ranks of the matrix c
        -- this is mode 3
      psolve: (M GR, L GR, PI) -> L Rec3
        ++ psolve(c,w,k) solves c z = w for all possible ranks >= k
        ++ of the matrix c and given right hand side vector w
        -- this is mode 4
      psolve: (M GR, L Symbol, PI) -> L Rec3
        ++ psolve(c,w,k) solves c z = w for all possible ranks >= k
        ++ of the matrix c and indeterminate right hand side w
        -- this is mode 5
      psolve: (M GR,       PI) -> L Rec3
        ++ psolve(c) solves the homogeneous linear system
        ++ c z = 0 for all possible ranks >= k of the matrix c
        -- this is mode 6
      psolve: (M GR, L GR, S) -> I
        ++ psolve(c,w,s) solves c z = w for all possible ranks
        ++ of the matrix c and given right hand side vector w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 7
      psolve: (M GR, L Symbol, S) -> I
        ++ psolve(c,w,s) solves c z = w for all possible ranks
        ++ of the matrix c and indeterminate right hand side w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 8
      psolve: (M GR,       S) -> I
        ++ psolve(c,s) solves c z = 0 for all possible ranks
        ++ of the matrix c and given right hand side vector w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 9
      psolve: (M GR, L GR, PI, S) -> I
        ++ psolve(c,w,k,s) solves c z = w for all possible ranks >= k
        ++ of the matrix c and given right hand side w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 10
      psolve: (M GR, L Symbol, PI, S) -> I
        ++ psolve(c,w,k,s) solves c z = w for all possible ranks >= k
        ++ of the matrix c and indeterminate right hand side w,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 11
      psolve: (M GR,       PI, S) -> I
        ++ psolve(c,k,s) solves c z = 0 for all possible ranks >= k
        ++ of the matrix c,
        ++ writes the results to a file named s, and returns the
        ++ number of regimes
        -- this is mode 12

      wrregime   : (L Rec3, S) -> I
        ++ wrregime(l,s) writes a list of regimes to a file named s
        ++ and returns the number of regimes written
      rdregime   : S -> L Rec3
        ++ rdregime(s) reads in a list from a file with name s

        -- for internal use only --
      -- these are exported so my other packages can use them

      bsolve: (M GR, L GF, NNI, S, Inputmode) -> Ranksolns
        ++ bsolve(c, w, r, s, m) returns a list of regimes and
        ++ solutions of the system c z = w for ranks at least r;
        ++ depending on the mode m chosen, it writes the output to
        ++ a file given by the string s.
      dmp2rfi: GR -> GF
        ++ dmp2rfi(p) converts p to target domain
      dmp2rfi: M GR -> M GF
        ++ dmp2rfi(m) converts m to target domain
      dmp2rfi: L GR -> L GF
        ++ dmp2rfi(l) converts l to target domain
      se2rfi:  L Symbol -> L GF
        ++ se2rfi(l) converts l to target domain
      pr2dmp: PR -> GR
        ++ pr2dmp(p) converts p to target domain
      hasoln: (Fgb, L GR) -> Rec8
        ++ hasoln(g, l) tests whether the quasi-algebraic set
        ++ defined by p = 0 for p in g and q ~= 0 for q in l
        ++ is empty or not and returns a simplified definition
        ++ of the quasi-algebraic set
        -- this is now done in QALGSET package
      ParCondList: (M GR,NNI) -> RankConds
        ++ ParCondList(c,r) computes a list of subdeterminants of each
        ++ rank >= r of the matrix c and returns
        ++ a groebner basis for the
        ++ ideal they generate
      redpps: (Linsoln, Fgb) -> Linsoln
        ++ redpps(s,g) returns the simplified form of s after reducing
        ++ modulo a groebner basis g



--               L O C A L    F U N C T I O N S

      B1solve: Linsys -> Linsoln
        ++ B1solve(s) solves the system (s.mat) z = s.vec
        ++ for the variables given by the column indices of s.cols
        ++ in terms of the other variables and the right hand side s.vec
        ++ by assuming that the rank is s.rank,
        ++ that the system is consistent, with the linearly
        ++ independent equations indexed by the given row indices s.rows;
        ++ the coefficients in s.mat involving parameters are treated as
        ++ polynomials.  B1solve(s) returns a particular solution to the
        ++ system and a basis of the homogeneous system (s.mat) z = 0.
      factorset: GR -> L GR
        ++ factorset(p) returns the set of irreducible factors of p.
      maxrank: RankConds -> NNI
        ++ maxrank(r) returns the maximum rank in the list r of regimes
      minrank: RankConds -> NNI
        ++ minrank(r) returns the minimum rank in the list r of regimes
      minset: L L GR -> L L GR
        ++ minset(sl) returns the sublist of sl consisting of the minimal
        ++ lists (with respect to inclusion) in the list sl of lists
      nextSublist: (I, I) -> L L I
        ++ nextSublist(n,k) returns a list of k-subsets of {1, ..., n}.
      overset?: (L GR, L L GR) -> Boolean
        ++ overset?(s,sl) returns true if s properly a sublist of a member
        ++ of sl; otherwise it returns false
      ParCond    : (M GR,NNI) -> Eqns
        ++ ParCond(m,k) returns the list of all k by k subdeterminants in
        ++ the matrix m
      redmat: (M GR, Fgb) -> M GR
        ++ redmat(m,g) returns a matrix whose entries are those of m
        ++ modulo the ideal generated by the groebner basis g
      regime: (Rec,M GR,L GF,L L GR,NNI,NNI,Inputmode) -> Rec3
        ++ regime(y,c, w, p, r, rm, m) returns a regime,
        ++ a list of polynomials specifying the consistency conditions,
        ++ a particular solution and basis representing the general
        ++ solution of the parametric linear system c z = w
        ++ on that regime. The regime returned depends on
        ++ the subdeterminant y.det and the row and column indices.
        ++ The solutions are simplified using the assumption that
        ++ the system has rank r and maximum rank rm. The list p
        ++ represents a list of list of factors of polynomials in
        ++ a groebner basis of the ideal generated by higher order
        ++ subdeterminants, and ius used for the simplification.
        ++ The mode m
        ++ distinguishes the cases when the system is homogeneous,
        ++ or the right hand side is arbitrary, or when there is no
        ++ new right hand side variables.
      sqfree: GR -> GR
        ++ sqfree(p) returns the product of square free factors of p
      inconsistent?: L GR -> Boolean
        ++ inconsistant?(pl) returns true if the system of equations
        ++ p = 0 for p in pl is inconsistent.  It is assumed
        ++ that pl is a groebner basis.
        -- this is needed because of change to
        -- EuclideanGroebnerBasisPackage
      inconsistent?: L PR -> Boolean
        ++ inconsistant?(pl) returns true if the system of equations
        ++ p = 0 for p in pl is inconsistent.  It is assumed
        ++ that pl is a groebner basis.
        -- this is needed because of change to
        -- EuclideanGroebnerBasisPackage

  Definition == add

    inconsistent?(pl:L GR):Boolean ==
      for p in pl repeat
        ground? p => return true
      false
    inconsistent?(pl:L PR):Boolean ==
      for p in pl repeat
        ground? p => return true
      false

    B1solve (sys:Linsys):Linsoln ==
      i,j,i1,j1:I
      rss:L I:=sys.rows
      nss:L I:=sys.cols
      k:=sys.rank
      cmat:M GF:=sys.mat
      n:=ncols cmat
      frcols:L I:=setDifference$(L I) (expand$(SEG I) (1..n), nss)
      w:L GF:=sys.vec
      p:V GF:=new(n,0)
      pbas:L V GF:=[]
      if k ~= 0 then
        augmat:M GF:=zero(k,n+1)
        for i in rss for i1 in 1.. repeat
          for j in nss for j1 in 1.. repeat
            augmat(i1,j1):=cmat(i,j)
          for j in frcols for j1 in k+1.. repeat
            augmat(i1,j1):=-cmat(i,j)
          augmat(i1,n+1):=w.i
        augmat:=rowEchelon$(M GF) augmat
        for i in nss for i1 in 1.. repeat p.i:=augmat(i1,n+1)
        for j in frcols for j1 in k+1.. repeat
          pb:V GF:=new(n,0)
          pb.j:=1
          for i in nss for i1 in 1.. repeat
            pb.i:=augmat(i1,j1)
          pbas:=cons(pb,pbas)
      else
        for j in frcols for j1 in k+1.. repeat
          pb:V GF:=new(n,0)
          pb.j:=1
          pbas:=cons(pb,pbas)
      [p,pbas]

    regime (y, coef, w, psbf, rk, rkmax, mode) ==
      i,j:I
      -- use the y.det nonzero to simplify the groebner basis
      -- of ideal generated by higher order subdeterminants
      ydetf:L GR:=factorset y.det
      yzero:L GR:=
        rk = rkmax => nil$(L GR)
        psbf:=[setDifference(x, ydetf) for x in psbf]
        groebner$gb [*/x for x in psbf]
      -- simplify coefficients by modulo ideal
      nc:M GF:=dmp2rfi redmat(coef,yzero)
      -- solve the system
      rss:L I:=y.rows;  nss:L I :=y.cols
      sys:Linsys:=[nc,w,rk,rss,nss]$Linsys
      pps:= B1solve(sys)
      pp:=pps.partsol
      frows:L I:=setDifference$(L I) (expand$(SEG I) (1..nrows coef),rss)
      wcd:L PR:= []
      -- case homogeneous rhs
      entry? (mode, [3,6,9,12]$(L I)) =>
               [yzero, ydetf,wcd, redpps(pps, yzero)]$Rec3
      -- case arbitrary rhs, pps not reduced
      for i in frows repeat
          weqn:GF:=+/[nc(i,j)*(pp.j) for j in nss]
          wnum:PR:=numer$GF (w.i - weqn)
          wnum = 0 => "trivially satisfied"
          ground? wnum => return [yzero, ydetf,[1$PR]$(L PR),pps]$Rec3
          wcd:=cons(wnum,wcd)
      entry? (mode, [2,5,8,11]$(L I)) => [yzero, ydetf, wcd, pps]$Rec3
      -- case no new rhs variable
      if not empty? wcd then _
        yzero:=removeDuplicates append(yzero,[pr2dmp pw for pw in wcd])
      test:Rec8:=hasoln (yzero, ydetf)
      not test.sysok => [test.z0, test.n0, [1$PR]$(L PR), pps]$Rec3
      [test.z0, test.n0, [], redpps(pps, test.z0)]$Rec3

    bsolve (coeff, w, h, outname, mode) ==
      r:=nrows coeff
--    n:=ncols coeff
      r ~= #w => error "number of rows unequal on lhs and rhs"
      newfile:FNAME
      rksoln:File Rec3
      count:I:=0
      lrec3:L Rec3:=[]
      filemode:Boolean:= entry? (mode, [7,8,9,10,11,12]$(L I))
      if filemode then
        newfile:=new$FNAME  ("",outname,"regime")
        rksoln:=open$(File Rec3) newfile
      y:Rec
      k:NNI
      rrcl:RankConds:=
        entry? (mode,[1,2,3,7,8,9]$(L I)) => ParCondList (coeff,0)
        entry? (mode,[4,5,6,10,11,12]$(L I)) => ParCondList (coeff,h)
      rkmax:=maxrank rrcl
      rkmin:=minrank rrcl
      for k in rkmax-rkmin+1..1 by -1 repeat
        rk:=rrcl.k.rank
        pc:Eqns:=rrcl.k.eqns
        psb:Fgb:= (if rk=rkmax then [] else rrcl.(k+1).fgb)
        psbf:L L GR:= [factorset x for x in psb]
        psbf:= minset(psbf)
        for y in pc repeat
          rec3:Rec3:= regime (y, coeff, w, psbf, rk, rkmax, mode)
          inconsistent? rec3.wcond => "incompatible system"
          if filemode then write_!(rksoln, rec3)
          else lrec3:= cons(rec3, lrec3)
          count:=count+1
      if filemode then close_! rksoln
      [lrec3, count]$Ranksolns

    factorset y ==
      ground? y => []
      [j.factor for j in factors(factor$mf y)]

    ParCondList (mat, h) ==
      rcl: RankConds:= []
      ps: L GR:=[]
      pc:Eqns:=[]
      npc: Eqns:=[]
      psbf: Fgb:=[]
      rc: Rec
      done: Boolean := false
      r:=nrows mat
      n:=ncols mat
      maxrk:I:=min(r,n)
      k:NNI
      for k in min(r,n)..h by -1 until done repeat
        pc:= ParCond(mat,k)
        npc:=[]
        (empty? pc) and (k >= 1) => maxrk:= k - 1
        if ground? pc.1.det -- only one is sufficient (neqzro = {})
        then (npc:=pc; done:=true; ps := [1$GR])
        else
          zro:L GR:= (if k = maxrk then [] else rcl.1.fgb)
          covered:Boolean:=false
          for rc in pc until covered repeat
            p:GR:= redPol$rp (rc.det, zro)
            p = 0 => "incompatible or covered subdeterminant"
            test:=hasoln(zro, [rc.det])
--          zroideal:=ideal(zro)
--          inRadical? (p, zroideal) => "incompatible or covered"
            ^test.sysok => "incompatible or covered"
-- The next line is WRONG! cannot replace zro by test.z0
--          zro:=groebner$gb (cons(*/test.n0, test.z0))
            zro:=groebner$gb (cons(p,zro))
            npc:=cons(rc,npc)
            done:= covered:= inconsistent? zro
          ps:=zro
        pcl: Rec2:= construct(k,npc,ps)
        rcl:=cons(pcl,rcl)
      rcl

    redpps(pps, zz) ==
      pv:=pps.partsol
      r:=#pv
      pb:=pps.basis
      n:=#pb + 1
      nummat:M GR:=zero(r,n)
      denmat:M GR:=zero(r,n)
      for i in  1..r repeat
        nummat(i,1):=pr2dmp numer$GF pv.i
        denmat(i,1):=pr2dmp denom$GF pv.i
      for j in 2..n repeat
        for i in 1..r  repeat
          nummat(i,j):=pr2dmp numer$GF (pb.(j-1)).i
          denmat(i,j):=pr2dmp denom$GF (pb.(j-1)).i
      nummat:=redmat(nummat, zz)
      denmat:=redmat(denmat, zz)
      for i in 1..r repeat
        pv.i:=(dmp2rfi nummat(i,1))/(dmp2rfi denmat(i,1))
      for j in 2..n repeat
        pbj:V GF:=new(r,0)
        for i in 1..r repeat
          pbj.i:=(dmp2rfi nummat(i,j))/(dmp2rfi  denmat(i,j))
        pb.(j-1):=pbj
      [pv, pb]

    dmp2rfi (mat:M GR): M GF ==
      r:=nrows mat
      n:=ncols mat
      nmat:M GF:=zero(r,n)
      for i in 1..r repeat
        for j in 1..n repeat
          nmat(i,j):=dmp2rfi mat(i,j)
      nmat

    dmp2rfi (vl: L GR):L GF ==
      [dmp2rfi v for v in vl]

    psolve (mat:M GR, w:L GR): L Rec3 ==
      bsolve(mat, dmp2rfi w, 1, "nofile", 1).rgl
    psolve (mat:M GR, w:L Symbol): L Rec3 ==
      bsolve(mat,  se2rfi w, 1, "nofile", 2).rgl
    psolve (mat:M GR): L Rec3 ==
      bsolve(mat, [0$GF for i in 1..nrows mat], 1, "nofile", 3).rgl

    psolve (mat:M GR, w:L GR, h:PI): L Rec3 ==
      bsolve(mat, dmp2rfi w, h::NNI, "nofile", 4).rgl
    psolve (mat:M GR, w:L Symbol, h:PI): L Rec3 ==
      bsolve(mat, se2rfi w, h::NNI, "nofile", 5).rgl
    psolve (mat:M GR, h:PI): L Rec3 ==
      bsolve(mat, [0$GF for i in 1..nrows mat], h::NNI, "nofile", 6).rgl

    psolve (mat:M GR, w:L GR, outname:S): I ==
      bsolve(mat, dmp2rfi w, 1, outname, 7).rgsz
    psolve (mat:M GR, w:L Symbol, outname:S): I ==
      bsolve(mat, se2rfi w, 1, outname, 8).rgsz
    psolve (mat:M GR, outname:S): I ==
      bsolve(mat, [0$GF for i in 1..nrows mat], 1, outname, 9).rgsz

    nextSublist (n,k) ==
      n <= 0 => []
      k <= 0 => [ nil$(List Integer) ]
      k > n => []
      n = 1 and k = 1 => [[1]]
      mslist: L L I:=[]
      for ms in nextSublist(n-1,k-1) repeat
        mslist:=cons(append(ms,[n]),mslist)
      append(nextSublist(n-1,k), mslist)

    psolve (mat:M GR, w:L GR, h:PI, outname:S): I ==
      bsolve(mat, dmp2rfi w, h::NNI, outname, 10).rgsz
    psolve (mat:M GR, w:L Symbol, h:PI, outname:S): I ==
      bsolve(mat, se2rfi w, h::NNI, outname, 11).rgsz
    psolve (mat:M GR, h:PI, outname:S): I ==
      bsolve(mat,[0$GF for i in 1..nrows mat],h::NNI,outname, 12).rgsz

    hasoln (zro,nzro) ==
      empty? zro => [true, zro, nzro]
      zro:=groebner$gb zro
      inconsistent? zro => [false, zro, nzro]
      empty? nzro =>[true, zro, nzro]
      pnzro:GR:=redPol$rp (*/nzro, zro)
      pnzro = 0 => [false, zro, nzro]
      nzro:=factorset pnzro
      psbf:L L GR:= minset [factorset p for p in zro]
      psbf:= [setDifference(x, nzro) for x in psbf]
      entry? ([], psbf) => [false, zro, nzro]
      zro:=groebner$gb [*/x for x in psbf]
      inconsistent? zro => [false, zro, nzro]
      nzro:=[redPol$rp (p,zro) for p in nzro]
      nzro:=[p for p in nzro | ^(ground? p)]
      [true, zro, nzro]



    se2rfi w == [coerce$GF monomial$PR (1$PR, wi, 1) for wi in w]

    pr2dmp p ==
      ground? p => (ground p)::GR
      algCoerceInteractive(p,PR,GR)$(Lisp) pretend GR

    wrregime (lrec3, outname) ==
      newfile:FNAME:=new$FNAME ("",outname,"regime")
      rksoln: File Rec3:=open$(File Rec3) newfile
      count:I:=0  -- number of distinct regimes
--      rec3: Rec3
      for rec3 in lrec3 repeat
          write_!(rksoln, rec3)
          count:=count+1
      close_!(rksoln)
      count

    dmp2rfi (p:GR):GF ==
      map$plift ((convert #1)@Symbol::GF, #1::PR::GF, p)


    rdregime inname ==
      infilename:=filename$FNAME ("",inname, "regime")
      infile: File Rec3:=open$(File Rec3) (infilename, "input")
      rksoln:L Rec3:=[]
      rec3:Union(Rec3, "failed"):=readIfCan_!$(File Rec3) (infile)
      while rec3 case Rec3 repeat
        rksoln:=cons(rec3::Rec3,rksoln) -- replace : to :: for AIX
        rec3:=readIfCan_!$(File Rec3) (infile)
      close_!(infile)
      rksoln

    maxrank rcl ==
      empty? rcl => 0
      "max"/[j.rank for j in rcl]

    minrank rcl ==
      empty? rcl => 0
      "min"/[j.rank for j in rcl]

    minset lset ==
      empty? lset => lset
      [x for x in lset | ^(overset?(x,lset))]

    sqfree p == */[j.factor for j in factors(squareFree p)]


    ParCond (mat, k) ==
      k = 0 => [[1, [], []]$Rec]
      j:NNI:=k::NNI
      DetEqn :Eqns := []
      r:I:= nrows(mat)
      n:I:= ncols(mat)
      k > min(r,n) => error "k exceeds maximum possible rank "
      found:Boolean:=false
      for rss in nextSublist(r, k) until found repeat
        for nss in nextSublist(n, k) until found repeat
          matsub := mat(rss, nss) pretend SM(j, GR)
          detmat := determinant(matsub)
          if detmat ~= 0 then
            found:= (ground? detmat)
            detmat:=sqfree detmat
            neweqn:Rec:=construct(detmat,rss,nss)
            DetEqn:=cons(neweqn, DetEqn)
      found => [first DetEqn]$Eqns
      sort(degree #1.det < degree #2.det, DetEqn)



    overset?(p,qlist) ==
      empty? qlist => false
      or/[part?(brace(q)$Set(GR), brace(p)$Set(GR))$Set(GR) _
                                                for q in qlist]


    redmat (mat,psb) ==
      i,j:I
      r:=nrows(mat)
      n:=ncols(mat)
      newmat: M GR:=zero(r,n)
      for i in 1..r repeat
        for j in 1..n repeat
          p:GR:=mat(i,j)
          ground? p => newmat(i,j):=p
          newmat(i,j):=redPol$rp (p,psb)
      newmat

@
<<*>>=
-- See LICENSE.AXIOM for Copyright information

<<package PLEQN ParametricLinearEquations>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}