1 /* Copyright (C) LinBox 2 * 3 * author: Zhendong Wan 4 * 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 26 27 #ifndef __LINBOX_minpoly_integer_H 28 #define __LINBOX_minpoly_integer_H 29 30 #include <iostream> 31 #include <cmath> 32 33 /*! \file algorithms/minpoly-integer.h 34 * Compute the minpoly of a matrix over an integer ring using modular arithmetic 35 * @todo better filter out repeated primes 36 */ 37 38 39 40 #include "linbox/field/field-traits.h" 41 #include "linbox/algorithms/matrix-hom.h" 42 #include "linbox/vector/vector-domain.h" 43 #include "linbox/randiter/random-prime.h" 44 //#include "linbox/solutions/minpoly.h" 45 #include "linbox/util/commentator.h" 46 #include <fflas-ffpack/ffpack/ffpack.h> 47 #include "linbox/algorithms/cra-builder-early-multip.h" 48 49 namespace LinBox 50 { 51 52 /* compute the minpoly of a matrix over the Integer ring 53 * via modular method over Field. 54 */ 55 template <class _Integer, class _Field> 56 class MinPoly { 57 public: 58 typedef _Field Field; 59 //typedef _Integer Integer; 60 typedef typename Field::Element Element; 61 62 /* Blackbox case */ 63 template<class Poly, class IMatrix> 64 static Poly& minPoly(Poly& y, const IMatrix& M); 65 66 template <class IMatrix> 67 static int minPolyDegree (const IMatrix& M, int n_try = 1); 68 69 template<class Poly, class IMatrix> 70 static Poly& minPoly(Poly& y, const IMatrix& M, int degree); 71 72 template<class Poly, class IMatrix> 73 static Poly& minPolyNonSymmetric(Poly& y, const IMatrix& M, int degree); 74 75 template<class Poly, class IMatrix> 76 static Poly& minPolySymmetric(Poly& y, const IMatrix& M, int degree); 77 78 template <class IMatrix> 79 static bool isSymmetric(const IMatrix& M, int n_try = 1); 80 }; 81 82 template <class _Integer, class _Field> 83 class MinPolyBlas { 84 public: 85 typedef _Field Field; 86 typedef _Integer Integer; 87 typedef typename Field::Element Element; 88 89 template <class Poly, class Ring> 90 static Poly& minPolyBlas (Poly& y, const BlasMatrix<Ring>& M); 91 92 template <class Poly, class Ring> 93 static Poly& minPolyBlas (Poly& y, const BlasMatrix<Ring>& M, int degree); 94 95 template <class Ring> 96 static int minPolyDegreeBlas (const BlasMatrix<Ring>& M, int n_try = 1); 97 }; 98 99 template<class _Integer, class _Field> 100 template<class Poly, class IMatrix> minPoly(Poly & y,const IMatrix & M)101 Poly& MinPoly<_Integer, _Field>::minPoly(Poly& y, const IMatrix& M) 102 { 103 int degree = minPolyDegree (M); 104 minPoly(y, M, degree); 105 return y; 106 } 107 108 template <class _Integer, class _Field> 109 template <class IMatrix> minPolyDegree(const IMatrix & M,int n_try)110 int MinPoly<_Integer, _Field>::minPolyDegree (const IMatrix& M, int n_try) 111 { 112 int degree = 0; 113 typedef typename IMatrix::template rebind<Field>::other FBlackbox; 114 typedef std::vector<Element> FPoly; 115 FPoly fp; 116 PrimeIterator<IteratorCategories::HeuristicTag> primeg(FieldTraits<Field>::bestBitSize(M.coldim())); 117 for (int i = 0; i < n_try; ++ i) { 118 ++primeg; 119 Field F(*primeg); 120 FBlackbox fbb(M, F); 121 minpoly (fp, fbb); 122 if (degree < ((int) fp.size() - 1)) degree = fp.size() -1; 123 } 124 return degree; 125 } 126 127 template <class _Integer, class _Field> 128 template<class Poly, class IMatrix> minPoly(Poly & y,const IMatrix & M,int degree)129 Poly& MinPoly<_Integer, _Field>::minPoly(Poly& y, const IMatrix& M, int degree) 130 { 131 if (isSymmetric(M)) { 132 //std::cout << "Symmetric:\n"; 133 minPolySymmetric(y, M, degree); 134 } 135 else { 136 //std::cout << "NonSymmetric:\n"; 137 minPolyNonSymmetric(y, M, degree); 138 } 139 return y; 140 } 141 142 template <class _Integer, class _Field> 143 template<class Poly, class IMatrix> minPolyNonSymmetric(Poly & y,const IMatrix & M,int degree)144 Poly& MinPoly<_Integer, _Field>::minPolyNonSymmetric(Poly& y, const IMatrix& M, int degree) 145 { 146 147 typedef typename IMatrix::template rebind<Field>::other FBlackbox; 148 typedef std::vector<Element> FPoly; 149 150 PrimeIterator<IteratorCategories::HeuristicTag> primeg(FieldTraits<Field>::bestBitSize(M.coldim())); 151 152 FPoly fp (degree + 1); 153 // typename FPoly::iterator fp_p; 154 y.resize (degree + 1); 155 156 CRABuilderEarlyMultip< _Field > cra(LINBOX_DEFAULT_EARLY_TERMINATION_THRESHOLD); 157 do { 158 ++primeg; 159 Field F(*primeg); 160 FBlackbox fbb(M, F); 161 minpoly (fp, fbb); 162 cra.initialize(F, fp); 163 } while( (int)fp.size() - 1 != degree); // Test for Bad primes 164 165 while(! cra.terminated()) { 166 ++primeg; while(cra.noncoprime(*primeg)) ++primeg; 167 Field F(*primeg); 168 FBlackbox fbb(M, F); 169 minpoly (fp, fbb); 170 if ((int)fp.size() - 1 != degree) { 171 commentator().report (Commentator::LEVEL_IMPORTANT, 172 INTERNAL_DESCRIPTION) << "Bad prime.\n"; 173 continue; 174 } 175 cra.progress(F, fp); 176 } 177 178 cra. result (y); 179 // commentator().report (Commentator::LEVEL_IMPORTANT, INTERNAL_DESCRIPTION) << "Number of primes needed: " << cra. steps() << std::endl; 180 return y; 181 } 182 183 template <class _Integer, class _Field> 184 template<class Poly, class IMatrix> minPolySymmetric(Poly & y,const IMatrix & M,int degree)185 Poly& MinPoly<_Integer, _Field>::minPolySymmetric(Poly& y, const IMatrix& M, int degree) 186 { 187 188 typedef typename IMatrix::template rebind<Field>::other FBlackbox; 189 typedef std::vector<Element> FPoly; 190 191 PrimeIterator<IteratorCategories::HeuristicTag> primeg(FieldTraits<Field>::bestBitSize(M.coldim())); 192 193 FPoly fp (degree + 1); 194 // typename FPoly::iterator fp_p; 195 y.resize (degree + 1); 196 197 CRABuilderEarlyMultip< _Field > cra(LINBOX_DEFAULT_EARLY_TERMINATION_THRESHOLD); 198 do { 199 ++primeg; 200 Field F(*primeg); 201 FBlackbox fbb(M,F); 202 minpolySymmetric (fp, fbb); 203 cra.initialize(F, fp); 204 } while( (int)fp.size() - 1 != degree); // Test for Bad primes 205 206 while(! cra.terminated()) { 207 ++primeg; while(cra.noncoprime(*primeg)) ++primeg; 208 Field F(*primeg); 209 FBlackbox fbb(M,F); 210 minpolySymmetric (fp, fbb); 211 if ((int)fp.size() - 1 != degree) { 212 commentator().report (Commentator::LEVEL_IMPORTANT, 213 INTERNAL_DESCRIPTION) << "Bad prime.\n"; 214 continue; 215 } 216 cra.progress(F, fp); 217 } 218 cra. result (y); 219 //std::cout << "Number of primes needed: " << cra. steps() << std::endl; 220 return y; 221 } 222 223 224 template <class _Integer, class _Field> 225 template <class IMatrix> isSymmetric(const IMatrix & M,int n_try)226 bool MinPoly<_Integer, _Field>::isSymmetric(const IMatrix& M, int n_try) 227 { 228 typedef typename IMatrix::Field Ring; 229 typedef typename Ring::Element Element; 230 Ring R(M. field()); int order = M. rowdim(); 231 std::vector<Element> x(order), mx (order), xm (order); 232 typename std::vector<Element>::iterator x_p; 233 VectorDomain<Ring> RVD (R); 234 if (M. rowdim() != M. coldim()) return false; 235 236 for (int i = 0; i < n_try; ++ i) { 237 for (x_p = x. begin(); x_p != x. end(); ++ x_p) 238 R. init (*x_p, rand()); 239 M. apply (mx, x); 240 M. applyTranspose (xm, x); 241 if (!RVD.areEqual(mx, xm)) return false; 242 } 243 return true; 244 } 245 246 template <class _Integer, class _Field> 247 template <class Poly, class Ring> minPolyBlas(Poly & y,const BlasMatrix<Ring> & M)248 Poly& MinPolyBlas<_Integer, _Field>::minPolyBlas (Poly& y, const BlasMatrix<Ring>& M) 249 { 250 int degree = minPolyDegreeBlas (M); 251 minPolyBlas (y, M, degree); 252 return y; 253 } 254 255 template <class _Integer, class _Field> 256 template <class Poly, class Ring> minPolyBlas(Poly & y,const BlasMatrix<Ring> & M,int degree)257 Poly& MinPolyBlas<_Integer, _Field>::minPolyBlas (Poly& y, const BlasMatrix<Ring>& M, int degree) 258 { 259 260 y. resize (degree + 1); 261 size_t n = M. rowdim(); 262 PrimeIterator<IteratorCategories::HeuristicTag> primeg(FieldTraits<Field>::bestBitSize(M.coldim())); 263 Element* FA = new Element [n*n]; 264 Element* X = new Element [n*(n+1)]; 265 size_t* Perm = new size_t[n]; 266 Element* p; 267 typename BlasMatrix<Ring>::ConstIterator raw_p; 268 std::vector<Element> poly (degree + 1); 269 // typename std::vector<Element>::iterator poly_ptr; 270 271 CRABuilderEarlyMultip< _Field > cra(LINBOX_DEFAULT_EARLY_TERMINATION_THRESHOLD); 272 do { 273 ++primeg; while(cra.noncoprime(*primeg)) ++primeg; 274 Field F(*primeg); 275 for (p = FA, raw_p = M. Begin(); 276 p != FA + (n*n); ++ p, ++ raw_p) 277 278 F. init (*p, *raw_p); 279 280 FFPACK::MinPoly((typename _Field::Father_t) F, poly, n, FA, n, X, n, Perm); 281 282 cra.initialize(F, poly); 283 } while( poly. size() != degree + 1) ; // Test for Bad primes 284 285 while (! cra. terminated()) { 286 ++primeg; while(cra.noncoprime(*primeg)) ++primeg; 287 Field F(*primeg); 288 for (p = FA, raw_p = M. Begin(); 289 p != FA + (n*n); ++ p, ++ raw_p) 290 291 F. init (*p, *raw_p); 292 293 FFPACK::MinPoly((typename _Field::Father_t) F, poly, n, FA, n, X, n, Perm); 294 295 if(poly. size() != degree + 1) { 296 commentator().report (Commentator::LEVEL_IMPORTANT, 297 INTERNAL_DESCRIPTION) << "Bad prime.\n"; 298 continue; 299 } 300 cra.progress(F, poly); 301 } 302 cra. result(y); 303 //std::cout << "Number of primes needed: " << cra. steps() << std::endl; 304 delete[] FA; delete[] X; delete[] Perm; 305 306 return y; 307 } 308 309 310 template <class _Integer, class _Field> 311 template <class Ring> minPolyDegreeBlas(const BlasMatrix<Ring> & M,int n_try)312 int MinPolyBlas<_Integer, _Field>::minPolyDegreeBlas (const BlasMatrix<Ring>& M, int n_try) 313 { 314 size_t n = M. rowdim(); 315 int degree = 0; 316 Element* FA = new Element [n*n]; 317 Element* X = new Element [n*(n+1)]; 318 size_t* Perm = new size_t[n]; 319 Element* p; 320 std::vector<Element> Poly; 321 322 PrimeIterator<IteratorCategories::HeuristicTag> primeg(FieldTraits<Field>::bestBitSize(M.coldim())); 323 324 typename BlasMatrix<Ring>::ConstIterator raw_p; 325 for (int i = 0; i < n_try; ++ i) { 326 ++primeg; 327 Field F(*primeg); 328 for (p = FA, raw_p = M. Begin(); 329 p!= FA + (n*n); ++ p, ++ raw_p) 330 F. init (*p, *raw_p); 331 332 FFPACK::MinPoly((typename _Field::Father_t) F, Poly, n, FA, n, X, n, Perm); 333 334 if (degree < Poly. size() - 1) 335 degree = Poly. size() -1; 336 } 337 delete[] FA; delete[] X; delete[] Perm; 338 339 return degree; 340 } 341 } // LinBox 342 343 #endif //__LINBOX_minpoly_integer_H 344 345 // Local Variables: 346 // mode: C++ 347 // tab-width: 4 348 // indent-tabs-mode: nil 349 // c-basic-offset: 4 350 // End: 351 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 352