1 /*
2 * examples/qchar.C
3 *
4 * Copyright (C) 2007, 2010 A. Urbanska
5 * ========LICENCE========
6 * This file is part of the library LinBox.
7 *
8 * LinBox is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU Lesser General Public
10 * License as published by the Free Software Foundation; either
11 * version 2.1 of the License, or (at your option) any later version.
12 *
13 * This library is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16 * Lesser General Public License for more details.
17 *
18 * You should have received a copy of the GNU Lesser General Public
19 * License along with this library; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
21 * ========LICENCE========
22 */
23
24 /*! @file examples/qchar.C
25 * @example examples/qchar.C
26 * @ingroup examples
27 * @brief Undocumented.
28 */
29
30 #include <givaro/modular.h>
31 #include <linbox/field/gmp-rational.h>
32 #include <linbox/matrix/dense-matrix.h>
33 #include <linbox/matrix/sparse-matrix.h>
34 #include <linbox/solutions/charpoly.h>
35 #include <linbox/solutions/minpoly.h>
36 #include <linbox/ring/polynomial-ring.h>
37
38 using namespace LinBox;
39 using namespace std;
40
41 typedef Givaro::ZRing<Integer> Integers;
42 typedef Integers::Element Integer;
43 typedef Givaro::Modular<double > Field;
44 typedef Field::Element Element;
45 typedef Vector<Integers>::Dense DVector;
46 typedef GMPRationalField Rationals;
47 typedef Rationals::Element Quotient;
48 typedef DenseMatrix<Integers > Blackbox;
49 typedef SparseMatrix<Rationals > RBlackbox;
50
51 //#define _LB_CONT_FR
52
53 typedef UserTimer myTimer;
54
55 myTimer tRationalModular;
56
57 template <class Field, class Polynomial>
printPolynomial(std::ostream & out,const Field & F,const Polynomial & v)58 std::ostream& printPolynomial (std::ostream& out, const Field &F, const Polynomial &v)
59 {
60 for (int i = v.size () - 1; i >= 0; i--) {
61 F.write (out, v[i]);
62 if (i > 0)
63 out << " X^" << i << " + ";
64 out << "\n";
65 }
66 return out;
67 }
68
69 template <class Blackbox, class MyMethod>
70 struct PrecRationalModularCharpoly {
71 const Blackbox &A;
72 const MyMethod &M;
73 const Integer ≺
74
PrecRationalModularCharpolyPrecRationalModularCharpoly75 PrecRationalModularCharpoly(const Integer& detPrec, const Blackbox& b, const MyMethod& n) :
76 prec(detPrec), A(b), M(n)
77 {}
78
79 template<typename Polynomial, typename Field>
operatorPrecRationalModularCharpoly80 Polynomial& operator()(Polynomial& P, const Field& F) const
81 {
82
83 typedef typename Blackbox::template rebind<Field>::other FBlackbox;
84 FBlackbox * Ap;
85
86 tRationalModular.clear();
87 tRationalModular.start();
88 MatrixHom::map(Ap, A, F);
89 tRationalModular.stop();
90 //minpoly(P, *Ap, typename FieldTraits<Field>::categoryTag(), M);
91 charpoly( P, *Ap, typename FieldTraits<Field>::categoryTag(), M);
92 //minpolySymmetric(P, *Ap, typename FieldTraits<Field>::categoryTag(), M);
93 integer p;
94 F.characteristic(p);
95 cout<<"Valence "<<p<<" = P[P.size()-1]" << P[P.size(0)-1] << " R-M time " ;
96 tRationalModular.print(cout);
97 cout << "\n" ;
98
99 typename Field::Element fprec;
100 F.init(fprec, prec);
101 //printPolynomial (std::cout, F, P);
102 delete Ap;
103
104 for (int i=0; i < P.size(); ++i)
105 F.mulin(P[i],fprec);
106 //std::cout<<"prec*Det(A) mod "<<p<<" = P[0]" << P[0] << "\n";
107 return P;
108 }
109 };
110
111 template <class Blackbox, class MyMethod>
112 struct PrecRationalModularMinpoly {
113 const Blackbox &A;
114 const MyMethod &M;
115 const DVector ≺
116
117
PrecRationalModularMinpolyPrecRationalModularMinpoly118 PrecRationalModularMinpoly(const DVector& detPrec, const Blackbox& b, const MyMethod& n) :
119 prec(detPrec), A(b), M(n)
120 {}
121
122 template<typename Polynomial, typename Field>
operatorPrecRationalModularMinpoly123 IterationResult operator()(Polynomial& P, const Field& F) const
124 {
125
126 typedef typename Blackbox::template rebind<Field>::other FBlackbox;
127 FBlackbox * Ap;
128
129 tRationalModular.clear();
130 tRationalModular.start();
131 MatrixHom::map(Ap, A, F);
132 tRationalModular.stop();
133 //minpoly(P, *Ap, typename FieldTraits<Field>::categoryTag(), M);
134 //charpoly( P, *Ap, typename FieldTraits<Field>::categoryTag(), M);
135 minpolySymmetric(P, *Ap, typename FieldTraits<Field>::categoryTag(), M);
136 integer p;
137 F.characteristic(p);
138 //cout<<"Det(A) mod "<<p<<" = P[0]" << P[0] << " R-M time " ;
139 cout<<"Valence "<<p<<" = P[P.size()-1]" << P[P.size()-1] << " R-M time " ;
140 tRationalModular.print(cout);
141 cout << "\n" ;
142
143 typename Field::Element fprec;
144 //F.init(fprec, prec);
145 //printPolynomial (std::cout, F, P);
146 delete Ap;
147
148 for (int i=0; i < P.size()-1; ++i) {
149 F.init(fprec, prec[i]);
150 F.mulin(P[i],fprec);
151 }
152 //std::cout<<"prec*Det(A) mod "<<p<<" = P[0]" << P[0] << "\n";
153 //printPolynomial(cout,F,P);
154 return IterationResult::CONTINUE;
155 }
156 };
157
158
159 void continuedFractionIn(double x,Integer& num, Integer& den, double epsi,const size_t s, Integer den_bound);
160 template <class RMatrix>
161 void generate_precRatMat(string& filename, RMatrix& M, DVector& den, Integer& denPrec);
162 void i_vti_v(RBlackbox& Res, DenseMatrix<Rationals >& M, DVector& den, Integer& detPrec);
163
main(int argc,char ** argv)164 int main (int argc, char** argv)
165 {
166
167 if ((argc < 3) || (argc > 4) ) {
168 cout << "Usage: qchar n matrix" << endl;
169 return -1;
170 }
171
172 size_t n = atoi(argv[1]);
173 cout << "Matrix size " << n<< endl;
174 string filename;
175 filename=argv[2];
176
177 Integers Z;
178 Rationals Q;
179
180 Integer detPrec = 1;
181 Integer detNum;
182 Integer detDen;
183
184 Blackbox DM(Z,n,n);
185 DenseMatrix<Rationals > In(Q,n,n);
186 RBlackbox M(Q,n,n);
187 DVector V(n,1);
188
189 // generate_precRatMat(filename, In, V, detPrec);
190 //i_vti_v(M,In, V, detPrec);
191 generate_precRatMat(filename, M, V, detPrec);
192 cout << "detPrec is " << detPrec << "\n";
193 typedef PolynomialRing<Integers> IntPolRing;
194 IntPolRing::Element c_A;
195
196 Integer max = 1,min=0;
197 for (size_t i=0; i < n; ++i) {
198 for (size_t j=0; j < n; ++j) {
199 Integer a;
200 Q.convert(a,M.getEntry(i,j));
201 // cerr<<"it="<<(*it)<<endl;
202 if (max < a)
203 max = a;
204 if (min > a)
205 min = a;
206 }
207 }
208 if (max<-min)
209 max=-min;
210 else max = max+1;
211 double hadamarcp = (double)n/2.0*(log(double(n))+2*log(double(max))+0.21163275)/log(2.0);
212
213
214 cout << "had" << hadamarcp << "\n";
215 cout << "had2" << (Integer)hadamarcp*detPrec << "\n";
216
217 PrimeIterator<IteratorCategories::HeuristicTag> genprime(FieldTraits<Field>::bestBitSize(M.coldim()));
218 ChineseRemainder< CRABuilderEarlyMultip<Field > > cra(LINBOX_DEFAULT_EARLY_TERMINATION_THRESHOLD);
219 typedef Method::Auto MyMethod;
220 MyMethod Met;
221 //PrecRationalModularCharpoly <RBlackbox ,MyMethod> iteration (detPrec, M, Met);
222 PrecRationalModularMinpoly< RBlackbox ,MyMethod> iteration(V, M, Met);
223 cra (c_A, iteration, genprime);
224 printPolynomial(cout, Z, c_A);
225
226 return 0 ;
227 }
228
229 /* for a rational number num/den construct a continued fraction approximation:
230 - with den bounded by den_bound
231 AND - number of steps bounded by s
232 NEW fraction replaces num/den
233 */
234
continuedFractionIn(Integer & num,Integer & den,double epsi,const size_t s,Integer den_bound)235 void continuedFractionIn(Integer& num, Integer& den, double epsi,const size_t s, Integer den_bound)
236 {
237 //den_bound = 100;
238 Integer a=num;
239 Integer b=den;
240 Integer t;
241 //double y;
242 Integer f[s];
243 f[0] = a/b;// y = f[0];
244 int i=1;
245 while (1) //abs(y-x)>=epsi)
246 {
247 //cout << "a :" << a << " b: " << b << "\n";
248 t = a%b;
249 a = b;
250 b = t;
251 f[i] = a/b;
252 //y = (double)f[i];
253 num = 1;
254 den = f[i];
255 for (int j=i-1; j > 0; --j) {
256 Integer tmp = num;
257 num = den;
258 den = f[j]*den+tmp;
259 //y = (double)f[j]+1/y;
260 }
261 //y = (double)f[0]+1/y;
262 num = den*f[0]+num;
263 if (den >= den_bound) break;
264 ++i;
265 if (i >= s) {
266 break;
267 }
268 }
269 }
270
i_vti_v(RBlackbox & Res,DenseMatrix<Rationals> & M,DVector & den,Integer & detPrec)271 void i_vti_v(RBlackbox& Res, DenseMatrix<Rationals >& M, DVector& den, Integer& detPrec)
272 {
273 detPrec=1;
274 Rationals Q;
275 int n = M.coldim();
276 //DVector denV(n,1);
277 for (int i=0; i < den.size();++i) den[i]=1;
278 for (int i=0; i < n; ++i) {
279 for (int j=0; j <=i ; ++j) {
280 Integer q_num =0;
281 Integer q_den =1;
282 for (int k=0; k < n; ++k) {
283 Integer deno_ik, nume_ik, deno_jk, nume_jk, deno, nume;
284 Q.get_den (deno_ik, M.getEntry(k,i));
285 Q.get_num (nume_ik, M.getEntry(k,i));
286 Q.get_den (deno_jk, M.getEntry(k,j));
287 Q.get_num (nume_jk, M.getEntry(k,j));
288
289 if (i==k) {
290 nume_ik = deno_ik-nume_ik;
291 }
292 else {
293 nume_ik = -nume_ik;
294 }
295
296 if (j==k) {
297 nume_jk = deno_jk-nume_jk;
298 }
299 else {
300 nume_jk = -nume_jk;
301 }
302 //cout << nume_ik << nume_jk;
303
304 deno = deno_ik*deno_jk;
305 nume = nume_ik*nume_jk;
306 //cout << q_num << "/" << q_den << " " ;
307 //cout << nume << "/" << deno << " " ;
308 if (deno==q_den) q_num+=nume;
309 else {
310 if (nume != 0) {
311 Integer g = gcd(q_den, deno);
312 q_num = q_num*deno/g+nume*q_den/g;
313 q_den = q_den*deno/g;
314 }
315 }
316 //cout << q_num << "/" << q_den << " " ;
317 }
318 if (q_num != 0) {
319 Quotient q(q_num,q_den);
320 if (i!=j) {
321 den[i]=lcm(den[i], q_den);
322 den[j]=lcm(den[j], q_den);
323 Res.setEntry(i,j,q);
324 Res.setEntry(j,i,q);
325 cout << i << " " << j << " " << q_num << "/" << q_den << "\n";
326 cout << j << " " << i << " " << q_num << "/" << q_den << "\n";
327 }
328 else {
329 den[i]=lcm(den[i], q_den);
330 Res.setEntry(i,j,q);
331 cout << i << " " << j << " " << q_num << "/" << q_den << "\n";
332 }
333 }
334 }
335 }
336 Integer d = 1;
337 for (int i=0; i < den.size(); ++i) {
338 detPrec *=den[i];
339 d = lcm(d, den[i]);
340 }
341 den[den.size()-1]=d;
342 for (int i=den.size()-2; i >=0; --i) {
343 den[i]=d*den[i+1];
344 }
345
346 }
347
348 template <class RMatrix>
generate_precRatMat(string & filename,RMatrix & M,DVector & den,Integer & denPrec)349 void generate_precRatMat(string& filename, RMatrix& M, DVector& den, Integer& denPrec)
350 {
351 int prec=10;
352 denPrec = 1;
353 ifstream is(filename.c_str(), std::ios::in);
354 int m,n;
355 char c;
356 is >> m >> n;
357 do is >> c; while (isspace (c));
358 if ((m > M.rowdim()) || (n > M.coldim())) {
359 cout << m << " " << n << " " ;
360 cout << "Wrong matrix size\n";
361 return;
362 }
363
364 cout << "Reading matrix\n";
365
366 den.resize(m,1);
367 while (1) {//! @todo temp fix
368 #if 0
369 while (is >> m) //temp fix
370 if (m==0) break;//temp fix
371 #endif
372 is >> m;//temp fix
373 is >> n;
374 #if 0
375 cout << m << n << "\n";
376 long double x;
377 is >> x;
378 #endif
379 char xstr[500];
380 char numstr[500];
381 #if 0
382 cout << x << " ";
383 sprintf(xstr, "%100f", x);
384 cout << xstr << " " ;
385 is >> c; int i=0;
386 while (c!= '\n') {
387 xstr[i]=c;
388 ++i;
389 is >> c;
390 if (i>=199) {
391 xstr[i]=0;
392 break;
393 }
394 }
395 is >> c;
396 #endif
397 is.getline(xstr,500);
398
399 Integer ten=1;
400 if (strchr(xstr, '/') != NULL) {
401 //cout << "xstr " << xstr << "\n";
402 //! @bug non reentrant strtok
403 strcpy(numstr, strtok(xstr, "/"));
404 char tenstr[500];
405 strcpy(tenstr, strtok(NULL, "/"));
406 if (prec < strlen(tenstr)) {
407 int cut = strlen(tenstr)-prec;
408 tenstr[prec]=0;
409 //c = numstr[strlen(numstr)-2];
410 numstr[strlen(numstr)-cut]=0;//temp round down
411 }
412 Integer tmp(tenstr);
413 ten = tmp;
414 //cout << numstr << "/" << ten << "\n";
415 }
416 else
417 {
418
419 // cout << "xstr" << xstr << " " ;
420 int i=0;
421 int j=0;
422 c=xstr[i]; ++i;
423 while (isspace (c)) {
424 if (i >=strlen(xstr)) break;
425 c=xstr[i];
426 ++i;
427 }
428 if (c=='-') {
429 numstr[j]=c;++j;
430 if (i < strlen(xstr)) {
431 c = xstr[i];
432 ++i;
433 }
434 }
435 if (c=='0') {
436 //cout << "zero";
437 if (i < strlen(xstr)) {
438 c = xstr[i];
439 ++i;
440 }
441 } //else cout << " not 0" << c;
442 while (strchr("0123456789",c) != NULL) {
443 //cout << "number";
444 numstr[j]=c; ++j;
445 if (i >=strlen(xstr)) break;
446 c=xstr[i];
447 ++i;
448 }
449 if (c=='.') {
450 if (i < strlen(xstr)) {
451 c = xstr[i];
452 //cout << "c" << c ;
453 ++i;
454 while (strchr("0123456789",c) != NULL) {
455 numstr[j]=c;++j;
456 if (i >=strlen(xstr)) break;
457 c=xstr[i];
458 ++i;
459 }
460 }
461 else {
462 cout << "wrong number format at .";
463 break;
464 }
465 }
466 numstr[j]=0;
467 //cout << "num" << numstr << " " ;
468
469 size_t dplaces=0;
470 int exp=0;
471 j=0;
472 c = xstr[j]; ++j;
473 while (c != '.') {
474 if (j >= strlen(xstr)) break;
475 c = xstr[j];
476 ++j;
477 }
478 if (j < strlen(xstr)) {
479 c = xstr[j]; ++j;
480 while (c != 'e') {
481 ++dplaces;
482 ten *= 10;
483 if (j >= strlen(xstr)) break;
484 c = xstr[j];
485 ++j;
486 }
487 }
488 if (j < strlen(xstr)) {
489 sscanf(xstr+j, "%d", &exp);
490 j=j-2;c = xstr[j];
491 while (c=='0') {
492 ten/=10;
493 --j;
494 c = xstr[j];
495 }
496 }
497
498 dplaces -=exp;
499 if (exp > 0) {
500 for (int i=0; i < exp; ++i) ten /=10;
501 }
502 else if (exp < 0) {
503 for (int i=0; i < exp; ++i) ten *=10;
504 }
505 //double nume = x * (double)ten;
506 }
507 Integer num (numstr);
508 Integer g;
509
510 #ifdef _LB_CONT_FR
511 double epsi = 1/(double)ten*10;
512 size_t s=10;
513 continuedFractionIn(num, ten, epsi, s, ten);
514 #endif
515 //cout << m << " " << n << " " << num << "/" << ten << "\n";
516 g = gcd(num,ten);
517 //g =1;
518 Quotient q(num/g,ten/g);
519 if ((m==0)&&(n==0)&&(num==0)) break;//temp fix
520 //M.setEntry(m-1,n-1,q);//temp fix
521 M.setEntry(m,n,q);
522 den[m-1] = lcm(den[m-1],ten/g);
523 }
524 Integer d = 1;
525 for (int i=0; i < den.size(); ++i) {
526 denPrec *=den[i];
527 d = lcm(d, den[i]);
528 }
529 den[den.size()-1]=d;
530 cout << d << "\n";
531 for (int i=den.size()-2; i >=0; --i) {
532 den[i]=d*den[i+1];
533 cout << den[i] << "\n";
534 }
535 }
536
537 #undef _LB_CONT_FR
538
539 // Local Variables:
540 // mode: C++
541 // tab-width: 4
542 // indent-tabs-mode: nil
543 // c-basic-offset: 4
544 // End:
545 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
546