1 /* linbox/blackbox/rational-reconstruction-base.h 2 * Copyright (C) 2009 Anna Marszalek 3 * 4 * Written by Anna Marszalek <aniau@astronet.pl> 5 * 6 * ========LICENCE======== 7 * This file is part of the library LinBox. 8 * 9 * LinBox is free software: you can redistribute it and/or modify 10 * it under the terms of the GNU Lesser General Public 11 * License as published by the Free Software Foundation; either 12 * version 2.1 of the License, or (at your option) any later version. 13 * 14 * This library is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 * Lesser General Public License for more details. 18 * 19 * You should have received a copy of the GNU Lesser General Public 20 * License along with this library; if not, write to the Free Software 21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22 * ========LICENCE======== 23 */ 24 25 #ifndef __LINBOX_charpoly_rational_H 26 #define __LINBOX_charpoly_rational_H 27 28 #include "linbox/util/commentator.h" 29 #include "linbox/util/timer.h" 30 #include "linbox/ring/modular.h" 31 32 //#include "linbox/field/gmp-rational.h" 33 #include "linbox/blackbox/rational-matrix-factory.h" 34 #include "linbox/algorithms/cra-builder-early-multip.h" 35 #include "linbox/algorithms/cra-domain.h" 36 //#include "linbox/algorithms/rational-cra.h" 37 #include "linbox/algorithms/rational-reconstruction-base.h" 38 #include "linbox/algorithms/classic-rational-reconstruction.h" 39 #include "linbox/solutions/charpoly.h" 40 #include "linbox/blackbox/compose.h" 41 #include "linbox/blackbox/diagonal.h" 42 43 namespace LinBox 44 { 45 //typedef GMPRationalField Rationals; 46 //typedef Rationals::Element Quotient; 47 48 /* 49 * Computes the characteristic polynomial of a rational dense matrix 50 */ 51 52 template<class T1, class T2> 53 struct MyModularCharpoly{ 54 T1* t1; 55 T2* t2; 56 57 int switcher; 58 59 MyModularCharpoly(T1* s1, T2* s2, int s = 1) {t1=s1; t2=s2;switcher = s;} 60 setSwitcherMyModularCharpoly61 int setSwitcher(int s) {return switcher = s;} 62 63 template<typename Polynomial, typename Field> operatorMyModularCharpoly64 IterationResult operator()(Polynomial& P, const Field& F) const 65 { 66 if (switcher ==1) { 67 return t1->operator()(P,F); 68 } 69 else { 70 return t2->operator()(P,F); 71 } 72 } 73 }; 74 75 template <class Blackbox, class MyMethod> 76 struct MyRationalModularCharpoly { 77 const Blackbox &A; 78 const MyMethod &M; 79 const std::vector<Integer> &mul;//multiplicative prec; 80 MyRationalModularCharpolyMyRationalModularCharpoly81 MyRationalModularCharpoly(const Blackbox& b, const MyMethod& n, 82 const std::vector<Integer >& p) : 83 A(b), M(n), mul(p) 84 {} MyRationalModularCharpolyMyRationalModularCharpoly85 MyRationalModularCharpoly(MyRationalModularCharpoly& C) : 86 // MyRationalModularCharpoly(C.A,C.M,C.mul) 87 A(C.A),M(C.M),mul(C.mul) 88 {} 89 90 template<typename Polynomial, typename Field> operatorMyRationalModularCharpoly91 IterationResult operator()(Polynomial& P, const Field& F) const 92 { 93 typedef typename Blackbox::template rebind<Field>::other FBlackbox; 94 FBlackbox * Ap; 95 MatrixHom::map(Ap, A, F); 96 charpoly( P, *Ap, typename FieldTraits<Field>::categoryTag(), M); 97 typename std::vector<Integer >::const_iterator it = mul.begin(); 98 typename Polynomial::iterator it_p = P.begin(); 99 for (;it_p !=P.end(); ++it, ++it_p) { 100 typename Field::Element e; 101 F.init(e, *it); 102 F.mulin(*it_p,e); 103 } 104 105 delete Ap; 106 return IterationResult::CONTINUE; 107 } 108 }; 109 110 template <class Blackbox, class MyMethod> 111 struct MyIntegerModularCharpoly { 112 const Blackbox &A; 113 const MyMethod &M; 114 const std::vector<typename Blackbox::Field::Element> &vD;//diagonal div. prec; 115 const std::vector<typename Blackbox::Field::Element > &mul;//multiplicative prec; 116 MyIntegerModularCharpolyMyIntegerModularCharpoly117 MyIntegerModularCharpoly(const Blackbox& b, const MyMethod& n, 118 const std::vector<typename Blackbox::Field::Element>& ve, 119 const std::vector<typename Blackbox::Field::Element >& p) : 120 A(b), M(n), vD(ve), mul(p) {} 121 MyIntegerModularCharpolyMyIntegerModularCharpoly122 MyIntegerModularCharpoly(MyIntegerModularCharpoly& C) : 123 // MyIntegerModularCharpoly(C.A,C.M,C.vD,C.mul) 124 A(C.A),M(C.M),vD(C.vD),mul(C.mul) 125 {} 126 127 template<typename Polynomial, typename Field> operatorMyIntegerModularCharpoly128 IterationResult operator()(Polynomial& P, const Field& F) const 129 { 130 typedef typename Blackbox::template rebind<Field>::other FBlackbox; 131 FBlackbox * Ap; 132 MatrixHom::map(Ap, A, F); 133 134 typename std::vector<typename Blackbox::Field::Element>::const_iterator it; 135 136 int i=0; 137 for (it = vD.begin(); it != vD.end(); ++it,++i) { 138 typename Field::Element t,tt; 139 F.init(t,*it); 140 F.invin(t); 141 for (int j=0; j < A.coldim(); ++j) { 142 F.mulin(Ap->refEntry(i,j),t); 143 } 144 } 145 146 charpoly( P, *Ap, typename FieldTraits<Field>::categoryTag(), M); 147 typename std::vector<typename Blackbox::Field::Element >::const_iterator it2 = mul.begin(); 148 typename Polynomial::iterator it_p = P.begin(); 149 for (;it_p !=P.end(); ++it2, ++it_p) { 150 typename Field::Element e; 151 F.init(e, *it2); 152 F.mulin(*it_p,e); 153 } 154 155 delete Ap; 156 return IterationResult::CONTINUE; 157 } 158 }; 159 160 template <class Rationals, template <class> class Vector, class MyMethod > 161 Vector<typename Rationals::Element>& rational_charpoly (Vector<typename Rationals::Element> &p, 162 const BlasMatrix<Rationals > &A, 163 const MyMethod &Met= Method::Auto()) 164 { 165 166 typedef typename Rationals::Element Quotient; 167 168 commentator().start ("Rational Charpoly", "Rminpoly"); 169 170 typedef Givaro::ModularBalanced<double> Field; 171 PrimeIterator<IteratorCategories::HeuristicTag> genprime(FieldTraits<Field>::bestBitSize(A.coldim())); 172 173 std::vector<Integer> F(A.rowdim()+1,1); 174 std::vector<Integer> M(A.rowdim()+1,1); 175 std::vector<Integer> Di(A.rowdim()); 176 177 RationalMatrixFactory<Givaro::ZRing<Integer>,Rationals,BlasMatrix<Rationals > > FA(&A); 178 Integer da=1, di=1; Integer D=1; 179 FA.denominator(da); 180 181 for (int i=(int)M.size()-2; i >= 0 ; --i) { 182 //c[m]=1, c[0]=det(A); 183 FA.denominator(di,i); 184 D *=di; 185 Di[(size_t)i]=di; 186 M[(size_t)i] = M[(size_t)i+1]*da; 187 } 188 for (int i=0; i < (int) M.size() ; ++i ) { 189 gcd(M[(size_t)i],M[(size_t)i],D); 190 } 191 192 Givaro::ZRing<Integer> Z; 193 BlasMatrix<Givaro::ZRing<Integer> > Atilde(Z,A.rowdim(), A.coldim()); 194 FA.makeAtilde(Atilde); 195 196 ChineseRemainder< CRABuilderEarlyMultip<Field> > cra(LINBOX_DEFAULT_EARLY_TERMINATION_THRESHOLD); 197 MyRationalModularCharpoly<BlasMatrix<Rationals > , MyMethod> iteration1(A, Met, M); 198 MyIntegerModularCharpoly<BlasMatrix<Givaro::ZRing<Integer> >, MyMethod> iteration2(Atilde, Met, Di, M); 199 MyModularCharpoly<MyRationalModularCharpoly<BlasMatrix<Rationals > , MyMethod>, 200 MyIntegerModularCharpoly<BlasMatrix<Givaro::ZRing<Integer> >, MyMethod> > iteration(&iteration1,&iteration2); 201 202 RReconstruction<Givaro::ZRing<Integer>, ClassicMaxQRationalReconstruction<Givaro::ZRing<Integer> > > RR; 203 204 std::vector<Integer> PP; // use of integer due to non genericity of cra. PG 2005-08-04 205 UserTimer t1,t2; 206 t1.clear(); 207 t2.clear(); 208 t1.start(); 209 cra(2,PP,iteration1,genprime); 210 t1.stop(); 211 t2.start(); 212 cra(2,PP,iteration2,genprime); 213 t2.stop(); 214 215 216 217 if (t1.time() < t2.time()) { 218 //cout << "ratim"; 219 iteration.setSwitcher(1); 220 } 221 else { 222 //cout << "intim"; 223 iteration.setSwitcher(2); 224 } 225 226 int k=4; 227 while (! cra(k,PP, iteration, genprime)) { 228 k *=2; 229 Integer m; //Integer r; Integer a,b; 230 cra.getModulus(m); 231 cra.result(PP);//need to divide 232 for (int i=0; i < (int)PP.size(); ++i) { 233 Integer D_1; 234 inv(D_1,M[(size_t)i],m); 235 PP[(size_t)i] = (PP[(size_t)i]*D_1) % m; 236 } 237 Integer den,den1; 238 std::vector<Integer> num(A.rowdim()+1); 239 std::vector<Integer> num1(A.rowdim()+1); 240 if (RR.reconstructRational(num,den,PP,m,-1)) {//performs reconstruction strating form c[m], use c[(size_t)i] as prec for c[(size_t)i-1] 241 cra(1,PP,iteration,genprime); 242 cra.getModulus(m); 243 for (int i=0; i < (int)PP.size(); ++i) { 244 Integer D_1; 245 inv(D_1,M[(size_t)i],m); 246 PP[(size_t)i] = (PP[(size_t)i]*D_1) % m; 247 } 248 if (RR.reconstructRational(num1,den1,PP,m,-1)) { 249 bool terminated = true; 250 if (den==den1) { 251 for (int i=0; i < (int)num.size(); ++i) { 252 if (num[(size_t)i] != num1[(size_t)i]) { 253 terminated =false; 254 break; 255 } 256 } 257 } 258 else { 259 terminated = false; 260 } 261 //set p 262 if (terminated) { 263 size_t i =0; 264 integer t,tt,ttt; 265 integer err; 266 // size_t max_err = 0; 267 Quotient qerr; 268 p.resize(PP.size()); 269 typename Vector <typename Rationals::Element>::iterator it; 270 Rationals Q; 271 for (it= p.begin(); it != p.end(); ++it, ++i) { 272 A.field().init(*it, num[(size_t)i],den); 273 Q.get_den(t,*it); 274 if (it != p.begin()) Q.get_den(tt,*(it-1)); 275 else tt = 1; 276 Q.init(qerr,t,tt); 277 278 } 279 return p; 280 // break; 281 } 282 } 283 } 284 } 285 286 cra.result(PP); 287 288 size_t i =0; 289 integer t,tt; 290 integer err; 291 // size_t max_res=0;int max_i; 292 // double rel; 293 // size_t max_resu=0; int max_iu; 294 // size_t max_err = 0; 295 Quotient qerr; 296 p.resize(PP.size()); 297 298 typename Vector <typename Rationals::Element>::iterator it; 299 300 Rationals Q; 301 for (it= p.begin(); it != p.end(); ++it, ++i) { 302 A.field().init(*it, PP[(size_t)i],M[(size_t)i]); 303 Q.get_den(t, *it); 304 Q.get_num(tt,*it); 305 err = M[(size_t)i]/t; 306 // size_t resi = err.bitsize() + tt.bitsize() -1; 307 // size_t resu = t.bitsize() + tt.bitsize() -1; 308 // if (resi > max_res) {max_res = resi; max_i=i;} 309 // if (resu > max_resu) {max_resu = resu; max_iu =i;} 310 //size_t resu = t.bitsize() + tt.bitsize() -1; 311 //if (err.bitsize() > max_err) max_err = err.bitsize(); 312 } 313 314 // max_res=0; 315 for (it= p.begin()+1; it != p.end(); ++it) { 316 //A.field().init(*it, PP[(size_t)i],M[(size_t)i]); 317 Q.get_den(t, *it); 318 Q.get_den(tt, *(it-1)); 319 Q.init(qerr,t,tt); 320 Q.get_num(tt, *it); 321 // size_t resi = Q.bitsize(t,qerr) + tt.bitsize() -2; 322 // if (resi > max_res) {max_res = resi; max_i=i;} 323 //if (err.bitsize() > max_err) max_err = err.bitsize(); 324 } 325 326 commentator().stop ("done", NULL, "Iminpoly"); 327 328 return p; 329 330 } 331 332 } 333 334 #endif //__LINBOX_charpoly_rational_H 335 336 // Local Variables: 337 // mode: C++ 338 // tab-width: 4 339 // indent-tabs-mode: nil 340 // c-basic-offset: 4 341 // End: 342 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 343