forked from fricas/fricas
-
Notifications
You must be signed in to change notification settings - Fork 0
/
muldep.spad
163 lines (143 loc) · 6.24 KB
/
muldep.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
)abbrev package GCDBAS GcdBasis
++ Gcd basis
++ Author: Waldek Hebisch
++ Date Created: 2009
++ Description:
++ Gcd basis provides functions to find structure of multiplicative
++ group (and semigroup) generated by elements of Gcd domain.
-- Currently the package uses naive quadratic algorithm: for each a_i it
-- computes gcd of a_i with all elements of current basis. Nontrivial gcd is
-- used to split elements.
GcdBasis(R : GcdDomain) : Exports == Implementation where
ResRec ==> Record(basis : Vector R, transform : Matrix Integer)
Exports ==> with
gcdBasis : Vector R -> Vector R
++ gcdBasis(v) returns basis part of \spad{gcdDecomposition(v)}.
gcdDecomposition : Vector R -> ResRec
++ gcdDecomposition(v) returns \spad{[b, t]} such that
++ elements of \spad{b} are relatively prime and that
++ \spad{v(i) = product(b(j)^(t(j, i)), j=1..n)}
gcdDecomposition : Vector Fraction R -> ResRec
++ gcdDecomposition(v) returns \spad{[b, t]} such that
++ elements of \spad{b} are relatively prime and that
++ \spad{v(i) = product(b(j)^(t(j, i)), j=1..n)}
-- multiplicativeBasis: Vector Fraction R -> ResRec
Implementation ==> add
gcdBasis(v) == gcdDecomposition(v).basis
SplitRec ==> Record(fac1 : R, fac2 : R, commonfac : R)
SplitRes ==> Union(SplitRec, "failed")
splitNums(a : R, b : R) : SplitRes ==
cf : R := gcd(a, b)
unit?(cf) => "failed"
a1 : R := (a exquo cf)::R
b1 : R := (b exquo cf)::R
[a1, b1, cf]
-- gcdDecomposition(vector([1/666, 11*7*9*5*128, 13*7*243*125*8]))
gcdDecomposition(v : Vector Fraction R) ==
n := #v
nv : Vector R := new(2*n, 0) -- $(Vector R)
for i in 1..n repeat
nv(i) := numer(v(i))
nv(i+n) := denom(v(i))
pr := gcdDecomposition(nv)
cb := pr.basis
ct0 := pr.transform
m := #cb
ct : Matrix Integer := new(m, n, 0)
for i in 1..m repeat
for j in 1..n repeat
ct(i, j) := ct0(i, j) - ct0(i, j + n)
[cb, ct]
-- gcdDecomposition(vector([11*7*9*5*128, 13*7*243*125*8]))
-- gcdDecomposition(vector([x^2-1, x^5-1]))
gcdDecomposition(v : Vector R) ==
cb : FlexibleArray R := empty()
n := #v
ct : FlexibleArray Vector Integer := empty()
for i in 1..n repeat
a := v(i)
i0 := #cb
for j in 1..i0 while not(unit?(a)) repeat
b := cb(j)
unit?(b) => iterate
pquo := a exquo b
if pquo case R then
while pquo case R repeat
ct(j)(i) := ct(j)(i) + 1
a := pquo
pquo := a exquo b
unit?(a) => break
sr1 := splitNums(a, b)
if sr1 case SplitRec then
sr : SplitRec := sr1
a1 := sr.fac1
b1 := sr.fac2
cf := sr.commonfac
if unit?(gcd(a1, cf)) and unit?(gcd(b1, cf)) then
concat!(cb, b1)
concat!(ct, copy(ct(j)))
cb(j) := cf
ct(j)(i) := ct(j)(i) + 1
a := a1
else
dr1 := gcdDecomposition(vector([a1, b1, cf]))
cb1 : Vector R := dr1.basis
ct1 := dr1.transform
firstAdded := true
a := 1::R
n1 := #cb1
ov := ct(j)
for k in 1..n1 repeat
c := ct1(k, 2) + ct1(k, 3)
if c > 0 then
if firstAdded then
firstAdded := false
cb(j) := cb1(k)
m := j
else
concat!(cb, cb1(k))
concat!(ct, empty())
m := #cb
nv := new(n, 0)$(Vector Integer)
for l in 1..n repeat
nv(l) := c*ov(l)
nv(i) := nv(i) + ct1(k, 1) + ct1(k, 3)
ct(m) := nv
else
for l in 1..(ct1(k, 1)) repeat
a := a*cb1(k)
if not(a = 1) then
concat!(cb, a)
nv := new(n, 0)$(Vector Integer)
nv(i) := 1
concat!(ct, nv)
[vector(parts cb), matrix([parts(ct(i)) for i in 1..#ct])]
)abbrev package MULDEP MultiplicativeDependence
MultiplicativeDependence() : Exports == Implementation where
FI ==> Fraction Integer
VZ ==> Vector Integer
RU ==> Union(Vector FI, "failed")
Exports ==> with
logDependenceQ : (List FI, FI) -> RU
++ logDependenceQ([q1, ..., qn], q0) finds rational
++ constants \spad{c1, ...cn} such that
++ \spad{q1^c1*...*qn^cn=u*q0} where \spad{u} is a unit
Implementation ==> add
import from GcdBasis(Integer)
solveOverQ(m : Matrix FI , v : Vector FI) : RU ==
particularSolution(m, v)$(LinearSystemMatrixPackage(FI,
Vector FI, Vector FI, Matrix FI))
logDependenceQ(lq : List FI, q : FI) : RU ==
empty?(lq) => "failed"
n := #lq
v := vector(concat(lq, [q]))$(Vector FI)
dr := gcdDecomposition(v)
bas := dr.basis
k0 := #bas
li := [i for i in 1..k0 | not(unit?(bas(i)))]
empty?(li) => "failed"
tm := dr.transform
-- k := nrows(tm)
m := matrix([[(tm(i, j))::FI for j in 1..n] for i in li])$(Matrix FI)
qv := vector([(tm(i, n+1))::FI for i in li])$(Vector FI)
solveOverQ(m, qv)