1 #include <linbox/linbox-config.h>
2
3 #include <iostream>
4 #include <string>
5 #include <vector>
6 #include <list>
7
8 #include <linbox/util/timer.h>
9
10 #include <linbox/ring/modular.h>
11 #include <linbox/matrix/sparse-matrix.h>
12 #include <linbox/algorithms/smith-form-sparseelim-local.h>
13
14 #include <linbox/util/timer.h>
15
16 #include <linbox/ring/local2_32.h>
17 #include <linbox/ring/pir-modular-int32.h>
18 #include <linbox/ring/givaro-poly.h>
19 #include <linbox/ring/givaro-poly-local.h>
20 #include <linbox/ring/givaro-poly-quotient.h>
21 #include <linbox/randiter/givaro-poly.h>
22 #include <linbox/algorithms/smith-form-local.h>
23 #include <linbox/algorithms/smith-form-iliopoulos2.h>
24 #include <linbox/algorithms/smith-form-adaptive.h>
25
26 using namespace LinBox;
27 using namespace std;
28
29 typedef Givaro::Modular<double> Field;
30 typedef typename Field::Element Element;
31 typedef Givaro::Poly1Dom<Field,Givaro::Dense> PolyDom;
32 typedef GivaroPoly<Field> PolyRing;
33 typedef GivaroPolyQuotient<PolyDom> QuotRing;
34 typedef GivaroPolyRandIter<PolyRing> PolyRandIter;
35 typedef typename PolyRing::Element PolyElement;
36 typedef MatrixDomain<PolyRing> MatDom;
37 typedef typename MatDom::OwnMatrix Mat;
38 typedef MatrixDomain<QuotRing> QuotMatDom;
39 typedef typename QuotMatDom::OwnMatrix QMat;
40
41 typedef GivaroPolyLocal<PolyDom> LocalRing;
42 typedef MatrixDomain<LocalRing> LocalMatDom;
43 typedef LocalMatDom::OwnMatrix LMat;
44
45 typedef IliopoulosDomain<PolyRing> IliopoulosDom;
46 typedef SmithFormLocal<QuotRing> LocalDom;
47 typedef DenseVector<PolyRing> FactorVector;
48
49 int VERBOSE = 0;
50
local(Mat & M,PolyElement & f,QuotRing & R)51 void local(Mat &M, PolyElement &f, QuotRing &R) {
52 size_t n = M.coldim();
53 QMat A(R, n, n);
54 for (size_t i = 0; i < n; i++) {
55 for (size_t j = 0; j < n; j++) {
56 PolyElement tmp;
57 M.getEntry(tmp, i, j);
58 A.setEntry(i, j, tmp);
59 }
60 }
61
62 LocalDom sfl;
63 list<PolyElement> l;
64
65 Timer timer;
66 timer.start();
67 sfl(l, A, R);
68 timer.stop();
69 cout << "Local Time: " << timer << endl;
70
71 if (VERBOSE) {
72 list<PolyElement>::const_iterator iterator;
73 for (iterator = l.begin(); iterator != l.end(); ++iterator) {
74 R.write(cout, *iterator) << endl;
75 }
76 cout << endl;
77 }
78 }
79
local2(Mat & M,LocalRing & R)80 void local2(Mat &M, LocalRing &R) {
81 size_t n = M.coldim();
82 LMat A(R, n, n);
83 for (size_t i = 0; i < n; i++) {
84 for (size_t j = 0; j < n; j++) {
85 PolyElement tmp;
86 M.getEntry(tmp, i, j);
87
88 LocalRing::Element ltmp;
89 R.init(ltmp, tmp);
90
91 A.setEntry(i, j, ltmp);
92 }
93 }
94
95 SmithFormLocal<LocalRing> sfl;
96 list<LocalRing::Element> l;
97
98 Timer timer;
99 timer.start();
100 sfl(l, A, R);
101 timer.stop();
102 cout << "Local2 Time: " << timer << endl;
103
104 if (VERBOSE) {
105 list<LocalRing::Element>::const_iterator iterator;
106 for (iterator = l.begin(); iterator != l.end(); ++iterator) {
107 R.write(cout, *iterator) << endl;
108 }
109 cout << endl;
110 }
111 }
112
ilio(Mat & M,PolyElement & f,MatDom & MD,PolyRing & R)113 void ilio(Mat &M, PolyElement &f, MatDom &MD, PolyRing &R) {
114 IliopoulosDom sfi(R);
115
116 size_t n = M.coldim();
117 Mat A(MD.field(), n, n);
118 MD.copy(A, M);
119
120 FactorVector l;
121 l.resize(n);
122
123 Timer timer;
124 timer.start();
125 sfi.smithForm(l, A, f);
126 timer.stop();
127 cout << "Iliopolous Time: " << timer << endl;
128
129 if (VERBOSE) {
130 for (size_t i = 0; i < l.size(); i++) {
131 R.write(cout, l[i]) << endl;
132 }
133 cout << endl;
134 }
135 }
136
generateM(Mat & M,MatDom & MD,PolyRing & R,PolyElement & f,PolyElement & g,size_t n,size_t e)137 void generateM(
138 Mat &M,
139 MatDom &MD,
140 PolyRing &R,
141 PolyElement &f,
142 PolyElement &g,
143 size_t n,
144 size_t e
145 ) {
146 integer max;
147 R.convert(max, f);
148 max *= 10;
149
150 Mat L(R, n, n);
151 for (size_t i = 0; i < n; i++) {
152 for (size_t j = 0; j < i; j++) {
153 PolyElement e;
154 R.init(e, rand() % max);
155 R.modin(e, f);
156 L.setEntry(i, j, e);
157 }
158 }
159 for (size_t i = 0; i < n; i++) {
160 L.setEntry(i, i, R.one);
161 }
162
163 Mat T(R, n, n);
164 for (size_t i = 0; i < n; i++) {
165 for (size_t j = i+1; j < n; j++) {
166 PolyElement e;
167 R.init(e, rand() % max);
168 R.modin(e, f);
169 T.setEntry(i, j, e);
170 }
171 }
172 for (size_t i = 0; i < n; i++) {
173 T.setEntry(i, i, R.one);
174 }
175
176 Mat D(R, n, n);
177 for (size_t i = 0; i < n; i++) {
178 size_t ei = rand() % e;
179 PolyElement di;
180 R.assign(di, R.one);
181
182 for (size_t j = 0; j < ei; j++) {
183 R.mulin(di, g);
184 }
185
186 D.setEntry(i, i, di);
187 if (VERBOSE > 1) {
188 R.write(cout << i << ", " << i << ": ", di) << endl;
189 }
190 }
191
192 MD.mul(M, L, D);
193 MD.mulin(M, T);
194
195 if (VERBOSE > 1) {
196 cout << endl;
197 for (size_t i = 0; i < n; i++) {
198 for (size_t j = 0; j < n; j++) {
199 PolyElement e;
200 M.getEntry(e, i, j);
201 R.write(cout << i << ", " << j << ": ", e) << endl;
202 }
203 }
204 }
205 }
206
main(int argc,char * argv[])207 int main(int argc, char* argv[]) {
208 size_t p = 3;
209 size_t n = 50;
210 size_t e = 6;
211 int gn = 13;
212
213 static Argument args[] = {
214 { 'p', "-p P", "Set the field GF(p)", TYPE_INT, &p},
215 { 'n', "-n N", "Set the matrix dimensions", TYPE_INT, &n},
216 { 'e', "-e E", "Set the max exponent", TYPE_INT, &e},
217 { 'g', "-g G", "Set irreducible such that f(p) = G", TYPE_INT, &gn},
218 { 'v', "-v V", "Set verbosity", TYPE_INT, &VERBOSE},
219 END_OF_ARGUMENTS
220 };
221
222 parseArguments(argc, argv, args);
223
224 Field F(p);
225 PolyRing R(F, "x");
226 PolyRandIter PRI(R, 0, 10);
227 MatDom MD(R);
228
229 PolyElement g, f;
230 R.init(g, gn);
231 R.write(cout << "irred: ", g) << endl;
232
233 R.pow(f, g, e);
234 R.write(cout << "quotient: ", f) << endl;
235
236 QuotRing QR(R, f);
237 LocalRing LR(R, g, e);
238
239 Mat M(R, n, n);
240 generateM(M, MD, R, f, g, n, e);
241
242 cout << endl;
243
244 ilio(M, f, MD, R);
245
246 //local(M, f, QR);
247
248 local2(M, LR);
249 }
250
251 // Local Variables:
252 // mode: C++
253 // tab-width: 4
254 // indent-tabs-mode: nil
255 // c-basic-offset: 4
256 // End:
257 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
258