-
Notifications
You must be signed in to change notification settings - Fork 46
/
skpol.spad
558 lines (473 loc) · 19.6 KB
/
skpol.spad
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
)abbrev category MORECAT MultivariateSkewPolynomialCategory
MultivariateSkewPolynomialCategory(R : Ring, E : OrderedAbelianMonoidSup,
Var : OrderedSet) : Category == MaybeSkewPolynomialCategory(R, E, Var)
)abbrev category SSKPOLC SolvableSkewPolynomialCategory
++ Description:
++ This is the category of polynomials in noncommutative
++ variables over noncommutative rings. We do not
++ assume that variables and elements of the base
++ ring commute. We assume that the polynomial ring
++ is of solvable type, so noncommutative version
++ of Buchberger algorithm works.
SolvableSkewPolynomialCategory(R : Ring,
Expon : OrderedAbelianMonoidSup) : Category == Join(Ring, LeftModule(R)) with
leadingCoefficient : % -> R
++ leadingCoefficient(p) returns the coefficient of the highest
++ degree term of p.
leadingMonomial : % -> %
++ leadingMonomial(p) returns the monomial of p with the highest degree.
degree : % -> Expon
++ degree(p) returns the maximum of the exponents of the terms of p.
reductum : % -> %
++ reductum(u) returns u minus its leading monomial
++ returns zero if handed the zero element.
monomial : (R, Expon) -> %
++ monomial(r, e) makes a term from a coefficient r and an exponent e.
)abbrev package NGROEB NGroebnerPackage
++ Description:
++ This is package computes noncommutative Groebner basis.
++ Based on commutative version. Note that this package
++ accepts rings as base domain, however computed basis
++ is over left fraction field. Computations are done
++ in fraction free way (coefficients stay in base ring).
NGroebnerPackage(Dom, Expon, Dpol) : T == C where
Dom : LeftOreRing
Expon : OrderedAbelianMonoidSup
Dpol : SolvableSkewPolynomialCategory(Dom, Expon)
NNI ==> NonNegativeInteger
critPair ==> Record( lcmfij : Expon, totdeg : NonNegativeInteger,
poli : Dpol, polj : Dpol )
sugarPol ==> Record( totdeg : NonNegativeInteger, pol : Dpol)
T == with
groebner : List(Dpol) -> List(Dpol)
++ groebner(lp) computes a groebner basis for a polynomial ideal
++ generated by the list of polynomials lp.
redPol : (Dpol, List(Dpol)) -> Dpol
++ normalForm(poly, gb) reduces the polynomial poly modulo the
++ precomputed groebner basis gb giving up to a constant factor
++ a canonical representative of the residue class.
hMonic : Dpol -> Dpol
++ hMonic(p) tries to remove content from p
virtualDegree : Dpol -> NonNegativeInteger
++ virtualDegree
sPol : critPair -> Dpol
++ sPol
C == add
import from OutputForm
credPol : (Dpol, List(Dpol)) -> Dpol
critM : (Expon, Expon) -> Boolean
redPo : (Dpol, List(Dpol) ) -> Record(poly : Dpol, mult : Dom)
-- hMonic: Dpol -> Dpol
updatF : (Dpol, NNI, List(sugarPol) ) -> List(sugarPol)
-- sPol: critPair -> Dpol
updatD : (List(critPair), List(critPair)) -> List(critPair)
critpOrder : (critPair, critPair) -> Boolean
makeCrit : (sugarPol, Dpol, NonNegativeInteger) -> critPair
-- virtualDegree : Dpol -> NonNegativeInteger
lcmCo : (Dom, Dom) -> Record(co1 : Dom, co2 : Dom)
lcmCo(c1, c2) ==
lcc := lcmCoef(c1, c2)
[lcc.coeff1, lcc.coeff2]
------ Definition of intermediate functions
if Dpol has totalDegree : Dpol -> NonNegativeInteger then
virtualDegree p == totalDegree p
else
virtualDegree p == 0
------ ordering of critpairs
critpOrder(cp1, cp2) ==
cp1.totdeg < cp2.totdeg => true
cp2.totdeg < cp1.totdeg => false
cp1.lcmfij < cp2.lcmfij
------ creating a critical pair
makeCrit(sp1, p2, totdeg2) ==
p1 := sp1.pol
deg := sup(degree(p1), degree(p2))
e1 := subtractIfCan(deg, degree(p1))::Expon
e2 := subtractIfCan(deg, degree(p2))::Expon
tdeg := max(sp1.totdeg + virtualDegree(monomial(1, e1)),
totdeg2 + virtualDegree(monomial(1, e2)))
[deg, tdeg, p1, p2]$critPair
------ calculate basis
gbasis(Pol : List(Dpol)) : List(Dpol) ==
D, D1 : List(critPair)
--------- create D and Pol
Pol1 := sort((z1, z2) +-> degree z1 > degree z2, Pol)
basPols := updatF(hMonic(first Pol1), virtualDegree(first Pol1), [])
Pol1 := rest(Pol1)
D := []
while not(empty?(Pol1)) repeat
h := hMonic(first(Pol1))
Pol1 := rest(Pol1)
toth := virtualDegree h
D1 := [makeCrit(x, h, toth) for x in basPols]
D := updatD(sort(critpOrder, D1), D)
basPols := updatF(h, toth, basPols)
D := sort(critpOrder, D)
redPols := [x.pol for x in basPols]
while not(empty?(D)) repeat
D0 := first D
s := hMonic(sPol(D0))
D := rest(D)
h := hMonic(redPol(s, redPols))
h = 0 => iterate
degree(h) = 0 =>
D := []
basPols := updatF(h, 0, [])
break
D1 := [makeCrit(x, h, D0.totdeg) for x in basPols]
D := updatD(sort(critpOrder, D1), D)
basPols := updatF(h, D0.totdeg, basPols)
redPols := concat(redPols, h)
Pol := [x.pol for x in basPols]
Pol
---- calculate minimal basis for ordered F
minGbasis(F : List(Dpol)) : List(Dpol) ==
empty?(F) => []
newbas := minGbasis rest F
cons(hMonic credPol( first(F), newbas), newbas)
groebner(F) ==
minGbasis(sort((x, y) +-> degree x > degree y, gbasis(F)))
--- calculate S-polynomial of a critical pair
sPol(p : critPair)==
Tij := p.lcmfij
fi := p.poli
fj := p.polj
fi := monomial(1, subtractIfCan(Tij, degree fi)::Expon)*fi
fj := monomial(1, subtractIfCan(Tij, degree fj)::Expon)*fj
cc := lcmCo(leadingCoefficient fi, leadingCoefficient fj)
cc.co1*fi - cc.co2*fj
redPo(s : Dpol, F : List(Dpol)) ==
m : Dom := 1
Fh := F
while not(s = 0 or empty?(F)) repeat
f1 := first(F)
s1 := degree(s)
e : Union(Expon, "failed")
d1 := degree(f1)
(e := subtractIfCan(s1, d1)) case Expon =>
f2 := monomial(1, e)*f1
cc := lcmCo(leadingCoefficient s, leadingCoefficient f2)
s := cc.co1*s - cc.co2*f2
if degree(s) = s1 then
print(message("no progress in reduction"))
print(s1::OutputForm)
print(d1::OutputForm)
print(e::OutputForm)
print(degree(f2)::OutputForm)
print(cc::OutputForm)
print(s::OutputForm)
error "no progress in reduction"
m := m*cc.co1
F := Fh
F := rest F
[s, m]
redPol(s : Dpol, F : List(Dpol)) ==
credPol(redPo(s, F).poly, F)
if Dpol has primitivePart : Dpol -> Dpol then
hMonic(p) == primitivePart p
else
hMonic(p) == p
--- reduce all terms of h mod F (iterative version )
credPol(h : Dpol, F : List(Dpol) ) ==
empty?(F) => h
h0 : Dpol := monomial(leadingCoefficient h, degree h)
while (h := reductum h) ~= 0 repeat
hred := redPo(h, F)
h := hred.poly
h0 := (hred.mult)*h0 + monomial(leadingCoefficient(h), degree h)
h0
--- concat ordered critical pair lists D1 and D2
updatD(D1 : List(critPair), D2 : List(critPair)) ==
empty?(D1) => D2
empty?(D2) => D1
dl1 := first(D1)
dl2 := first(D2)
critpOrder(dl1, dl2) => cons(dl1, updatD(D1.rest, D2))
cons(dl2, updatD(D1, D2.rest))
--- crit M - true, if lexp#2 multiple of lexp#1
critM(e1 : Expon, e2 : Expon) ==
en : Union(Expon, "failed")
(en := subtractIfCan(e2, e1)) case Expon
--- concat F and h and erase polys in F with lexp
--- being multiple of lexp h
--- Can do this because such poly can be recoverd
--- from critical pair
updatF(h : Dpol, deg : NNI, F : List(sugarPol)) ==
empty?(F) => [[deg, h]]
f1:= first(F)
critM(degree(h), degree(f1.pol)) => updatF(h, deg, rest(F))
cons(f1, updatF(h, deg, rest(F)))
)abbrev domain SKSMP SparseMultivariateSkewPolynomial
++ Description:
++ SparseMultivariateSkewPolynomial(R, Var, sigma, delta) defines
++ a mutivariate Ore ring over R in variables from V. \spad{sigma(v)}
++ gives automorphism of R corresponding to variable \spad{v} and
++ \spad{delta(v)} gives corresponding derivative.
++
SparseMultivariateSkewPolynomial(R : Ring, Var : OrderedSet,
sigma : Var -> Automorphism R , delta : Var -> (R -> R)
) : Exports == Implementation where
Exports ==> MultivariateSkewPolynomialCategory(R, IndexedExponents(Var),
Var) with
D : Var -> %
++ D(v) returns operator corresponding to derivative with
++ respect to v in R.
if Var has variable : Symbol -> Var then
Delta : Symbol -> %
++ Delta(s) returns operator corresponding to derivative with
++ respect to s in R.
Smp ==> SparseMultivariatePolynomial(R, Var)
Implementation ==> Smp add
D(v : Var) : % == monomial(1, v, 1)
if Var has with (variable : Symbol -> Union(Var, "failed")) then
Delta(s : Symbol) : % ==
vu := variable(s)
vu case "failed" => error "missing variable"
D(vu@Var)
Sup := SparseUnivariatePolynomial(%)
Supp := SparseUnivariatePolynomial(Smp)
-- Fake domain so that we can use UnivariateSkewPolynomialCategoryOps
id : Automorphism(%) := morphism((x +-> x)@(% -> %))
Upol := SparseUnivariateSkewPolynomial(%, id, (x : %) : % +-> 0)
lift_map(f : R -> R) : % -> % ==
x +-> map(f, x)
lift_morphism(m : Automorphism R) : Automorphism % ==
morphism(lift_map((x : R) : R +-> m(x)),
lift_map((x : R) : R +-> inv(m)(x)))
import from UnivariateSkewPolynomialCategoryOps(%, Upol)
x : % * y : % ==
ground?(x) => retract(x)@R*y
v := mainVariable(x)::Var
xu := univariate(x pretend Smp, v)$Smp pretend Upol
yu := univariate(y pretend Smp, v)$Smp pretend Upol
ru := times(xu, yu, lift_morphism(sigma(v)), lift_map(delta(v)))
multivariate(ru pretend Supp, v)$Smp pretend %
x : % * r : R == x*(r::%)
if R has IntegralDomain then
x : % exquo y : % ==
ground?(y) => x exquo retract(y)@R
ground?(x) => "failed"
maxd := totalDegree(x)
ly := leadingMonomial(y) pretend Smp
res : % := 0
while x ~= 0 repeat
lx1 := leadingMonomial(x)
totalDegree(lx1) > maxd => return "failed"
lx := lx1 pretend Smp
cu := lx exquo ly
cu case "failed" => return "failed"
cc := (cu@Smp) pretend %
res := res + cc
x := x - cc*y
res
x : % ^ n : PositiveInteger ==
res := x
for i in 2..n repeat res := res*x
res
x : % ^ n : NonNegativeInteger ==
n = 0 => 1
x^(qcoerce(n)@PositiveInteger)
coerce(x : %) : OutputForm ==
ground?(x) => (retract(x)@R)::OutputForm
v := mainVariable(x)::Var
xu := univariate(x pretend Smp, v) pretend Upol
outputForm(xu, sub('D::OutputForm, v::OutputForm))
)abbrev package PDOHLP PartialDifferentialOperatorHelper
PartialDifferentialOperatorHelper(R : PartialDifferentialRing(Var),
Var : OrderedSet) : with
id_map : Var -> Automorphism(R)
diff_map : Var -> (R -> R)
== add
id_map(v) == 1$Automorphism(R)
diff_map(v) == (x : R) : R +-> differentiate(x, v)
)abbrev domain PDO PartialDifferentialOperator
++ Description:
++ PartialDifferentialOperator(R, V) defines a ring of partial
++ differential operators in variables from V and with coefficients
++ in a partial differential ring R. Multiplication of operators
++ corresponds to composition of operators.
++
PartialDifferentialOperator(R : PartialDifferentialRing(Var),
Var : OrderedSet) : Exports == Implementation where
HLP ==> PartialDifferentialOperatorHelper(R, Var)
Exports ==> MultivariateSkewPolynomialCategory(R,
IndexedExponents(Var), Var) with
D : Var -> %
++ D(v) returns the operator corresponding to derivative with
++ respect to v in R.
adjoint : % -> %
++ adjoint(p) returns the adjoint of operator p.
if R has LeftOreRing then
SolvableSkewPolynomialCategory(R, IndexedExponents(Var))
Implementation ==> SparseMultivariateSkewPolynomial(R, Var, id_map$HLP,
diff_map$HLP) add
Smp ==> SparseMultivariatePolynomial(R, Var)
Sup := SparseUnivariatePolynomial(%)
adjoint(x : %) : % ==
ground?(x) => x
v := mainVariable(x)::Var
xu := univariate(x pretend Smp, v)$Smp pretend Sup
res : % := 0
xu := map(adjoint, xu)$Sup
while xu ~= 0 repeat
d := degree(xu)
sign := (odd?(d) => -1 ; 1)
res := res + sign*monomial(1, v, d)*leadingCoefficient(xu)
xu := reductum(xu)
res
)abbrev package OREMAT OrePolynomialMatrixOperations
OrePolynomialMatrixOperations(F, LO) : Exports == Implementation where
F : Field
LO : UnivariateSkewPolynomialCategory(F)
M ==> Matrix LO
V ==> Vector F
Param_Rec_F ==> Record(ratpart : F, coeffs : Vector F)
Param_Rec_V ==> Record(ratpart : V, coeffs : Vector F)
FPL ==> Record(particular : List Param_Rec_F, basis : List F)
VPL ==> Record(particular : List Param_Rec_V, basis : List V)
Exports ==> with
rowEchelon : M -> M
++ \spad{rowEchelon(m)} returns the row echelon form of the matrix m.
rowEchelon : (M, F, List V) -> Record(mat : M, vecs : List V)
++ \spad{rowEchelon(m, c, lv)} returns [m2, lv2] such that
++ m2 is the row echelon form of the matrix m. lv2
++ is transformed lv using c as parameter to \spad{apply}.
solve : (M, F, List V, (LO, List F) -> FPL) -> Union(VPL, "failed")
++ solve(m, c, lv, solf) returns "failed" of
++ [[vp1, ..., vp_m], [b1, ..., bl]]
++ such that bi-s are basis of solutions of homogeneous system
++ m bi = 0. Each vpi = [r, [c1, ..., cn]] is a particular solution
++ of a parametric matrix equation m r = \sum ci vi where
++ lv = [v1, ..., vn]. solf is scalar solver, c is a
++ parameter to \spad{apply} (needed for action of LO on F).
++ "failed" means that system is underdetermined.
Implementation ==> add
import from LinearCombinationUtilities(F, SparseUnivariatePolynomial F)
rowEchelon(m : M) == rowEchelon(m, 0, []).mat
rowEchelon(m, cc, lv) ==
mm := copy m
lw := [copy(vv) for vv in lv]
nc := ncols(mm)
nr := nrows(mm)
for w in lw repeat
if #w ~= nc then
error "rowEchelon: elements of lv must have size = ncols(m)"
i : Integer := 1
for j in 1..nc repeat
if i > nr then break
pivk : Integer := -1
dpiv : Integer := -1
mij : LO
for k in i..nr repeat
if not zero?(mkj := mm(k, j)) and (dp1 := degree(mkj);
(dpiv = -1 or dp1 < dpiv)) then
mij := mkj
pivk := k
dpiv := dp1
dpiv = - 1 => iterate
mm := swapRows!(mm, i, pivk)
for w in lw repeat
swap!(w, i, pivk)
for k in (i + 1)..nr | not zero?(mkj := mm(k, j)) repeat
(mij, c, d, u, v) := right_ext_ext_GCD(mij, mkj)
for k1 in (j + 1)..nc repeat
el1 := qelt(mm, i, k1)
el2 := qelt(mm, k, k1)
qsetelt!(mm, i, k1, c*el1 + d*el2)
qsetelt!(mm, k, k1, u*el1 + v*el2)
for w in lw repeat
elf1 := qelt(w, i)
elf2 := qelt(w, k)
qsetelt!(w, i, apply(c, cc, elf1) + apply(d, cc, elf2))
qsetelt!(w, k, apply(u, cc, elf1) + apply(v, cc, elf2))
qsetelt!(mm, i, j, mij)
qsetelt!(mm, k, j, 0)
for k in 1..(i - 1) repeat
zero?(mkj := mm(k, j)) => iterate
(q, r) := rightDivide(mkj, mij)
qsetelt!(mm, k, j, r)
for k1 in (j + 1)..nc repeat
qsetelt!(mm, k, k1, qelt(mm, k, k1) - q*qelt(mm, i, k1))
for w in lw repeat
qsetelt!(w, k, qelt(w, k) - apply(q, cc, qelt(w, i)))
i := i + 1
[mm, lw]
triangular_solve(m : M, cc : F, lv : List V,
solf : (LO, List F) -> FPL) : Union(VPL, "failed") ==
nc := ncols(m)
nr := nrows(m)
nc > nr => "failed"
n := #lv
bvl : List V := []
cb := [new(n, 0)$Vector(F) for i in 1..n]
cba := [new(nc, 0)$Vector(F) for bv in cb]
for i in 1..n for bv in cb repeat
bv(i) := 1
for i in nc..1 by -1 repeat
m(i, i) = 0 => return "failed"
lvi0 := [v(i) for v in lv]
lvi1 := [lin_comb(bv, lvi0) for bv in cb]
nlvi := #lvi1
nbas := #bvl
lvi1 := concat(lvi1, [0 for bv in bvl])
lvi2 : List(F) := []
for bv in concat(cba, bvl) for vvi in lvi1 repeat
for j in i+1..nc repeat
vvi := vvi - apply(m(i, j), cc, bv(j))
lvi2 := cons(vvi, lvi2)
lvi := reverse!(lvi2)
resi := solf(m(i, i), lvi)
nbvl : List V := []
for be in resi.basis repeat
bv := new(nc, 0)
bv(i) := be
nbvl := cons(bv, nbvl)
prl := resi.particular
nrl := #prl
ncvl := [pr.coeffs for pr in prl]
nsl1 := [pr.ratpart for pr in prl]
cm1 := matrix([parts(ncv) for ncv in ncvl])$Matrix(F)
cm2 := rowEchelon(horizConcat(cm1, scalarMatrix(nrl, 1)))
nsl2 : List(F) := []
noff := nlvi + nbas + 1
for j in 1..nrl repeat
ss : F := 0
for k in noff.. for slk in nsl1 repeat
ss := ss + cm2(j, k)*slk
nsl2 := cons(ss, nsl2)
nrl1 : NonNegativeInteger := 0
for j in nrl..1 by -1 for slj in nsl2 repeat
for k in 1..nlvi repeat
if cm2(j, k) ~= 0 then
nrl1 := j
break
nrl1 ~= 0 => break
nbcv := new(nbas, 0)$Vector(F)
for k in 1..nbas repeat
nbcv(k) := cm2(j, k + nlvi)
bv := lin_comb(nbcv, bvl)
bv(i) := slj
nbvl := cons(bv, nbvl)
nsl2 := reverse!(nsl2)
ncb : List(Vector(F)) := []
ncba : List(Vector(F)) := []
for j in 1..nrl1 for slj in nsl2 repeat
nbcv1 := new(nbas, 0)$Vector(F)
for k in 1..nbas repeat
nbcv1(k) := cm2(j, k + nlvi)
nbcv2 := new(nlvi, 0)$Vector(F)
for k in 1..nlvi repeat
nbcv2(k) := cm2(j, k)
bv := lin_comb(nbcv2, cb)
ncb := cons(bv, ncb)
na := lin_comb(nbcv2, cba)
na := lin_comb!(nbcv1, na, bvl)
na(i) := slj
ncba := cons(na, ncba)
cb := reverse!(ncb)
cba := reverse!(ncba)
bvl := nbvl
[[[ba, bv] for ba in cba for bv in cb], bvl]
solve(m, cc, lv, solf) ==
rec := rowEchelon(m, cc, lv)
triangular_solve(rec.mat, cc, rec.vecs, solf)