1 /* linbox/solutions/charpoly.h 2 * Copyright (C) 2005 Clement Pernet 3 * 4 * Written by Clement Pernet <clement.pernet@imag.fr> 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_H 26 #define __LINBOX_charpoly_H 27 28 #include "linbox/solutions/methods.h" 29 #include "linbox/util/debug.h" 30 #include "linbox/field/field-traits.h" 31 #include "linbox/matrix/dense-matrix.h" 32 #include "linbox/matrix/matrix-domain.h" 33 #include "linbox/randiter/random-prime.h" 34 #include "linbox/algorithms/bbcharpoly.h" 35 // Namespace in which all LinBox library code resides 36 37 namespace LinBox 38 { 39 40 // for specialization with respect to the DomainCategory 41 template< class Blackbox, class Polynomial, class MyMethod, class DomainCategory> 42 Polynomial& charpoly ( Polynomial &P, 43 const Blackbox &A, 44 const DomainCategory &tag, 45 const MyMethod &M); 46 47 48 /* //error handler for rational domain 49 template <class Blackbox, class Polynomial> 50 Polynomial &charpoly (Polynomial& P, 51 const Blackbox& A, 52 const RingCategories::RationalTag& tag, 53 const Method::Auto& M) 54 { 55 throw LinboxError("LinBox ERROR: charpoly is not yet defined over a rational domain"); 56 } 57 */ 58 /** \brief ...using an optional Method parameter 59 \param P - the output characteristic polynomial. If the polynomial 60 is of degree d, this random access container has size d+1, the 0-th 61 entry is the constant coefficient and the d-th is 1 since the charpoly 62 is monic. 63 \param A - a blackbox matrix 64 Optional \param M - the method object. Generally, the default 65 object suffices and the algorithm used is determined by the class of M. 66 Basic methods are Method::Blackbox, Method::Elimination, and 67 Method::Auto (the default). 68 See methods.h for more options. 69 \return a reference to P. 70 */ 71 template <class Blackbox, class Polynomial, class MyMethod> charpoly(Polynomial & P,const Blackbox & A,const MyMethod & M)72 Polynomial& charpoly (Polynomial & P, 73 const Blackbox & A, 74 const MyMethod & M){ 75 return charpoly ( P, A, typename FieldTraits<typename Blackbox::Field>::categoryTag(), M); 76 } 77 78 79 /// \brief ...using default method 80 template<class Blackbox, class Polynomial> charpoly(Polynomial & P,const Blackbox & A)81 Polynomial& charpoly (Polynomial & P, 82 const Blackbox & A) 83 { 84 return charpoly (P, A, Method::Auto()); 85 } 86 87 // The charpoly with Auto Method 88 template <class Blackbox, class Polynomial> charpoly(Polynomial & P,const Blackbox & A,const RingCategories::ModularTag & tag,const Method::Auto & M)89 Polynomial& charpoly (Polynomial & P, 90 const Blackbox & A, 91 const RingCategories::ModularTag & tag, 92 const Method::Auto & M) 93 { 94 //return charpoly(P, A, tag, Method::Blackbox(M)); 95 return charpoly(P, A, tag, Method::DenseElimination(M)); 96 } 97 98 // The charpoly with Auto Method 99 template <class Domain, class Polynomial> charpoly(Polynomial & P,const SparseMatrix<Domain> & A,const RingCategories::ModularTag & tag,const Method::Auto & M)100 Polynomial& charpoly (Polynomial & P, 101 const SparseMatrix<Domain> & A, 102 const RingCategories::ModularTag & tag, 103 const Method::Auto & M) 104 { 105 return charpoly(P, A, tag, Method::Blackbox(M)); 106 } 107 108 // The charpoly with Auto Method 109 template<class Domain, class Polynomial> charpoly(Polynomial & P,const BlasMatrix<Domain> & A,const RingCategories::ModularTag & tag,const Method::Auto & M)110 Polynomial& charpoly (Polynomial & P, 111 const BlasMatrix<Domain> & A, 112 const RingCategories::ModularTag & tag, 113 const Method::Auto & M) 114 { 115 return charpoly(P, A, tag, Method::DenseElimination(M)); 116 } 117 118 // The charpoly with Elimination Method 119 template<class Blackbox, class Polynomial> charpoly(Polynomial & P,const Blackbox & A,const RingCategories::ModularTag & tag,const Method::Elimination & M)120 Polynomial& charpoly (Polynomial & P, 121 const Blackbox & A, 122 const RingCategories::ModularTag & tag, 123 const Method::Elimination & M) 124 { 125 return charpoly(P, A, tag, Method::DenseElimination(M)); 126 } 127 128 129 /** @brief Compute the characteristic polynomial over \f$\mathbf{Z}_p\f$. 130 * 131 * Compute the characteristic polynomial of a matrix using dense 132 * elimination methods 133 134 * @param P Polynomial where to store the result 135 * @param A Blackbox representing the matrix 136 * @param tag 137 * @param M 138 */ 139 template <class Blackbox, class Polynomial > charpoly(Polynomial & P,const Blackbox & A,const RingCategories::ModularTag & tag,const Method::DenseElimination & M)140 Polynomial& charpoly (Polynomial & P, 141 const Blackbox & A, 142 const RingCategories::ModularTag & tag, 143 const Method::DenseElimination & M) 144 { 145 if (A.coldim() != A.rowdim()) 146 throw LinboxError("LinBox ERROR: matrix must be square for characteristic polynomial computation\n"); 147 148 BlasMatrix< typename Blackbox::Field > BBB (A); 149 BlasMatrixDomain< typename Blackbox::Field > BMD (BBB.field()); 150 //BlasVector<typename Blackbox::Field,Polynomial> P2(A.field(),P); 151 BMD.charpoly (P, BBB); 152 return P;//= P2.getRep() ; 153 } 154 155 156 } 157 158 #include "linbox/algorithms/matrix-hom.h" 159 160 #include "linbox/algorithms/rational-cra-var-prec.h" 161 #include "linbox/algorithms/cra-builder-var-prec-early-multip.h" 162 #include "linbox/algorithms/charpoly-rational.h" 163 164 namespace LinBox 165 { 166 template <class Blackbox, class MyMethod> 167 struct IntegerModularCharpoly { 168 const Blackbox &A; 169 const MyMethod &M; 170 IntegerModularCharpolyIntegerModularCharpoly171 IntegerModularCharpoly(const Blackbox& b, const MyMethod& n) : 172 A(b), M(n) 173 {} 174 175 template<typename Field, class Polynomial> operatorIntegerModularCharpoly176 IterationResult operator()(Polynomial& P, const Field& F) const 177 { 178 typedef typename Blackbox::template rebind<Field>::other FBlackbox; 179 FBlackbox Ap(A, F); 180 charpoly (P, Ap, typename FieldTraits<Field>::categoryTag(), M); 181 return IterationResult::CONTINUE; 182 // std::cerr << "Charpoly(A) mod "<<F.characteristic()<<" = "<<P; 183 // integer p; 184 // F.characteristic(p); 185 // std::cerr<<"Charpoly(A) mod "<<p<<" = "<<P; 186 } 187 }; 188 189 template <class Blackbox, class Polynomial> charpoly(Polynomial & P,const Blackbox & A,const RingCategories::IntegerTag & tag,const Method::Auto & M)190 Polynomial& charpoly (Polynomial & P, 191 const Blackbox & A, 192 const RingCategories::IntegerTag & tag, 193 const Method::Auto & M) 194 { 195 commentator().start ("Integer Charpoly", "Icharpoly"); 196 if (useBlackboxMethod(A)) 197 charpoly(P, A, tag, Method::Blackbox(M) ); 198 else 199 charpoly(P, A, tag, Method::DenseElimination(M) ); 200 commentator().stop ("done", NULL, "Icharpoly"); 201 return P; 202 } 203 204 } 205 206 #if defined(__LINBOX_HAVE_NTL) 207 208 #include "linbox/algorithms/cia.h" 209 namespace LinBox 210 { 211 /** @brief Compute the characteristic polynomial over {\bf Z} 212 * 213 * Compute the characteristic polynomial of a matrix using dense 214 * elimination methods 215 216 * @param P Polynomial where to store the result 217 * @param A \ref Black-Box representing the matrix 218 */ 219 220 221 template <class Blackbox, class Polynomial> charpoly(Polynomial & P,const Blackbox & A,const RingCategories::IntegerTag & tag,const Method::DenseElimination & M)222 Polynomial& charpoly (Polynomial & P, 223 const Blackbox & A, 224 const RingCategories::IntegerTag & tag, 225 const Method::DenseElimination & M) 226 { 227 if (A.coldim() != A.rowdim()) 228 throw LinboxError("LinBox ERROR: matrix must be square for characteristic polynomial computation\n"); 229 return cia (P, A, M); 230 } 231 232 /** Compute the characteristic polynomial over {\bf Z} 233 * 234 * Compute the characteristic polynomial of a matrix, represented via 235 * a blackBox. 236 * 237 * @param P Polynomial where to store the result 238 * @param A \ref Black-Box representing the matrix 239 */ 240 template <class Blackbox, class Polynomial/*, class Categorytag*/ > charpoly(Polynomial & P,const Blackbox & A,const RingCategories::IntegerTag & tag,const Method::Blackbox & M)241 Polynomial& charpoly (Polynomial & P, 242 const Blackbox & A, 243 const RingCategories::IntegerTag & tag, 244 const Method::Blackbox & M) 245 { 246 if (A.coldim() != A.rowdim()) 247 throw LinboxError("LinBox ERROR: matrix must be square for characteristic polynomial computation\n"); 248 return BBcharpoly::blackboxcharpoly (P, A, tag, M); 249 } 250 251 } 252 253 #else // no NTL 254 255 #include "linbox/ring/modular.h" 256 #include "linbox/algorithms/cra-domain.h" 257 #include "linbox/algorithms/cra-builder-full-multip.h" 258 #include "linbox/algorithms/cra-builder-early-multip.h" 259 #include "linbox/algorithms/matrix-hom.h" 260 261 namespace LinBox 262 { 263 264 template <class Matrix, class Polynomial, class Method> charpoly(Polynomial & P,const Matrix & A,const RingCategories::IntegerTag & tag,const Method & M)265 Polynomial& charpoly (Polynomial & P, 266 const Matrix & A, 267 const RingCategories::IntegerTag & tag, 268 const Method & M) 269 { 270 if (A.coldim() != A.rowdim()) 271 throw LinboxError("LinBox ERROR: matrix must be square for characteristic polynomial computation\n"); 272 273 commentator().start ("Integer Charpoly : chinese remaindering", "IntCharpoly"); 274 275 typedef Givaro::ModularBalanced<double> Field; 276 PrimeIterator<IteratorCategories::HeuristicTag> genprime(FieldTraits<Field>::bestBitSize(A.coldim())); 277 278 // @todo: use a value for the switch provided by the method and not by a macro 279 #ifdef __LINBOX_HEURISTIC_CRA 280 ChineseRemainder< CRABuilderEarlyMultip<Field > > cra(LINBOX_DEFAULT_EARLY_TERMINATION_THRESHOLD); 281 #else 282 double hbound = FastCharPolyHadamardBound(A); 283 ChineseRemainder< CRABuilderFullMultip<Field > > cra(hbound); 284 #endif 285 IntegerModularCharpoly<Matrix, Method> iteration(A, M); 286 cra.operator() (P, iteration, genprime); 287 commentator().stop ("done", NULL, "IbbCharpoly"); 288 #ifdef _LB_CRATIMING 289 cra.reportTimes(std::clog) << std::endl; 290 #endif 291 return P; 292 } 293 294 } 295 296 #endif 297 298 namespace LinBox 299 { 300 /** Compute the characteristic polynomial over \f$\mathbf{Z}_p\f$. 301 * 302 * Compute the characteristic polynomial of a matrix, represented via 303 * a blackBox. 304 * 305 * @param P Polynomial where to store the result 306 * @param A Blackbox representing the matrix 307 * @param tag 308 * @param M 309 */ 310 template <class Blackbox, class Polynomial/*, class Categorytag*/ > charpoly(Polynomial & P,const Blackbox & A,const RingCategories::ModularTag & tag,const Method::Blackbox & M)311 Polynomial& charpoly (Polynomial & P, 312 const Blackbox & A, 313 const RingCategories::ModularTag & tag, 314 const Method::Blackbox & M) 315 { 316 if (A.coldim() != A.rowdim()) 317 throw LinboxError("LinBox ERROR: matrix must be square for characteristic polynomial computation\n"); 318 319 return BBcharpoly::blackboxcharpoly (P, A, tag, M); 320 } 321 322 template < class Blackbox, class Polynomial, class MyMethod> charpoly(Polynomial & P,const Blackbox & A,const RingCategories::RationalTag & tag,const MyMethod & M)323 Polynomial &charpoly (Polynomial& P, const Blackbox& A, 324 const RingCategories::RationalTag& tag, const MyMethod& M) 325 { 326 commentator().start ("Rational Charpoly", "Rcharpoly"); 327 328 typedef Givaro::ModularBalanced<double> Field; 329 PrimeIterator<IteratorCategories::HeuristicTag> genprime(FieldTraits<Field>::bestBitSize(A.coldim())); 330 RationalChineseRemainderVarPrec< CRABuilderVarPrecEarlyMultip<Field > > rra(LINBOX_DEFAULT_EARLY_TERMINATION_THRESHOLD); 331 IntegerModularCharpoly<Blackbox,MyMethod> iteration(A, M); 332 333 Givaro::ZRing<Integer> Z; 334 DensePolynomial<Givaro::ZRing<Integer> > PP(Z); // use of integer due to non genericity of cra. PG 2005-08-04 335 Integer den; 336 rra(dynamic_cast<decltype(PP)::Domain_t::Storage_t&>(PP),den, iteration, genprime); 337 size_t i =0; 338 P.resize(PP.size()); 339 for (typename Polynomial::iterator it= P.begin(); it != P.end(); ++it, ++i) 340 A.field().init(*it, PP[i],den); 341 342 commentator().stop ("done", NULL, "Rcharpoly"); 343 344 return P; 345 } 346 347 } // end of LinBox namespace 348 #endif // __LINBOX_charpoly_H 349 350 // Local Variables: 351 // mode: C++ 352 // tab-width: 4 353 // indent-tabs-mode: nil 354 // c-basic-offset: 4 355 // End: 356 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 357