1 /* 2 * Copyright (C) 2013 Pascal Giorgi 3 * Romain Lebreton 4 * 5 * Written by Pascal Giorgi <pascal.giorgi@lirmm.fr> 6 * Romain Lebreton <lebreton@lirmm.fr> 7 * 8 * ========LICENCE======== 9 * This file is part of the library LinBox. 10 * 11 * LinBox is free software: you can redistribute it and/or modify 12 * it under the terms of the GNU Lesser General Public 13 * License as published by the Free Software Foundation; either 14 * version 2.1 of the License, or (at your option) any later version. 15 * 16 * This library is distributed in the hope that it will be useful, 17 * but WITHOUT ANY WARRANTY; without even the implied warranty of 18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 19 * Lesser General Public License for more details. 20 * 21 * You should have received a copy of the GNU Lesser General Public 22 * License along with this library; if not, write to the Free Software 23 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 24 * ========LICENCE======== 25 */ 26 #ifndef __LINBOX_polynomial_matrix_H 27 #define __LINBOX_polynomial_matrix_H 28 #include <vector> 29 #include "linbox/vector/subvector.h" 30 #include "linbox/matrix/dense-matrix.h" 31 #include "linbox/field/hom.h" 32 #include "fflas-ffpack/utils/align-allocator.h" 33 #include "givaro/modular.h" 34 #include <algorithm> 35 36 #ifdef TRACK_MEMORY_MATPOL 37 uint64_t max_memory=0, cur_memory=0; 38 #define ADD_MEM(x) {cur_memory+=x; max_memory=std::max(max_memory,cur_memory);} 39 #define DEL_MEM(x) {cur_memory-=x;} 40 #define STR_MEMINFO std::right<<"\033[31m [ MEM: cur="<<cur_memory/1000000.<<" Mo --- max="<<max_memory/1000000.<<" Mo \033[0m]" 41 #define PRINT_MEMINFO std::cerr<<"\033[31m[ MEM: cur="<<cur_memory/1000000.<<" Mo --- max="<<max_memory/1000000.<<" Mo ]\033[0m"<<std::endl; 42 #else 43 #define ADD_MEM(X) ; 44 #define DEL_MEM(X) ; 45 #endif 46 47 #define COPY_BLOCKSIZE 32 48 49 namespace LinBox{ 50 51 enum PMType {polfirst, matfirst}; 52 enum PMStorage {plain, view, const_view}; 53 54 // Generic handler class for Polynomial Matrix 55 template<size_t type, size_t storage, class Field> 56 class PolynomialMatrix; 57 element_storage(const Field & F)58 template<typename Field> uint64_t element_storage(const Field& F) { integer p;F.characteristic(p); return length(p);} element_storage(const Givaro::Modular<Givaro::Integer> & F)59 template<> uint64_t element_storage(const Givaro::Modular<Givaro::Integer> &F) { integer p;F.characteristic(p); return length(p)+sizeof(Givaro::Integer);} 60 61 // Class for Polynomial Matrix stored as a Matrix of Polynomials 62 template<class _Field> 63 class PolynomialMatrix<PMType::polfirst,PMStorage::plain,_Field> { 64 public: 65 typedef _Field Field; 66 typedef typename Field::Element Element; 67 typedef BlasMatrix<Field> Matrix; 68 //typedef typename std::vector<Element>::iterator Iterator; 69 //typedef typename std::vector<Element>::const_iterator ConstIterator; 70 typedef std::vector<Element,AlignedAllocator<Element, Alignment::DEFAULT>> VECT; 71 typedef typename VECT::iterator Iterator; 72 typedef typename VECT::const_iterator ConstIterator; 73 //typedef vector<Element> Polynomial; 74 typedef Subvector<Iterator,ConstIterator> Polynomial; 75 typedef PolynomialMatrix<PMType::polfirst,PMStorage::plain,_Field> Self_t; 76 typedef PolynomialMatrix<PMType::matfirst,PMStorage::plain,_Field> Other_t; 77 78 //PolynomialMatrix() {} 79 80 // construct a polynomial matrix in f[x]^(m x n) of degree (s-1) 81 PolynomialMatrix(const Field& f, size_t r, size_t c, size_t s, size_t stor=0) : 82 _store((stor?stor:s)), _repview(r*c),_rep(r*c*_store,f.zero), _row(r), _col(c), _size(s), _fld(&f) { 83 for (size_t i=0;i<_row;i++) 84 for (size_t j=0;j<_col;j++) 85 _repview[i*_col+j]= Polynomial(_rep.begin()+(i*_col+j)*_store,_size); 86 //integer p; 87 //std::cout<<"MatrixP allocating : "<<r*c*s*length(f.characteristic(p))/1000000.<<"Mo"<<std::endl; 88 //std::cout<<"(ALLOC) PolynomialMatrix<polfirst> at "<<this<<" : "<<r<<"x"<<c<<" - size= "<<s<<" ==> "<<MB(realmeminfo())<<" Mo "<<STR_MEMINFO<<std::endl; 89 ADD_MEM(realmeminfo()); 90 } 91 92 PolynomialMatrix(const Self_t&) = delete; 93 ~PolynomialMatrix()94 ~PolynomialMatrix(){ 95 DEL_MEM(realmeminfo()); 96 _rep.clear(); 97 //std::cout<<"(FREE) PolynomialMatrix<polfirst> at "<<this<<" : "<<_row<<"x"<<_col<<" - size= "<<_store<<" ==> "<<MB(realmeminfo())<<" Mo "<<STR_MEMINFO<<std::endl; 98 99 //integer p; 100 //std::cout<<"MatrixP Desallocating : "<<_row*_col*_store*length(_fld->characteristic(p))/1000000.<<"Mo"<<std::endl; 101 102 } 103 clear()104 void clear(){ 105 _rep.resize(0); 106 } 107 108 // set the matrix A to the matrix of degree k in the polynomial matrix setMatrix(Matrix & A,size_t k)109 void setMatrix(Matrix& A, size_t k) const { 110 typename Matrix::Iterator it=A.Begin(); 111 for(size_t i=0;i<_row*_col;i++,it++) 112 *it = get(i,k); 113 } 114 // retrieve the matrix of degree k in the polynomial matrix 115 Matrix operator[](size_t k)const { 116 Matrix A(field(), _row, _col); 117 setMatrix(A,k); 118 return A; 119 } 120 // retrieve the polynomial at entry (i,j) in the matrix operator()121 inline Polynomial& operator()(size_t i, size_t j) {return operator()(i*_col+j);} operator()122 inline const Polynomial& operator()(size_t i, size_t j)const {return operator()(i*_col+j);} 123 // retrieve the polynomial at the position i in the storage of the matrix operator()124 inline Polynomial& operator()(size_t i) {return _repview[i];} operator()125 inline const Polynomial& operator()(size_t i)const {return _repview[i];} 126 127 // resize the polynomial length of the polynomial matrix resize(size_t s)128 void resize(size_t s){ 129 if (s==_store) return; 130 //std::cout<<"MATPOL RESIZING : "<<_store<<" --> "<<s<<std::endl; 131 if (s>_store){ 132 _rep.resize(s*_row*_col); 133 size_t k=s*_row*_col-1; 134 for(size_t i=0;i<_row*_col;i++){ 135 size_t j=s; 136 for(;j>=_size;j--,k--) 137 _rep[k]=_fld->zero; 138 for(;j>size_t(-1);j--,k--) 139 _rep[k]=_rep[i*_store+j]; 140 } 141 } 142 else { 143 size_t k=0; 144 for(size_t i=0;i<_row*_col;i++) 145 for (size_t j=0;j<s;j++,k++) 146 _rep[k]=_rep[i*_store+j]; 147 _rep.resize(s*_row*_col); 148 //_rep.shrink_to_fit(); 149 } 150 integer p;_fld->characteristic(p); size_t bb=p.bitsize(); if(bb>64) bb+=128; bb/=8; 151 152 #ifdef MEM_PMBASIS 153 size_t mem=realmeminfo(); 154 ADD_MEM(realmeminfo()); 155 DEL_MEM(mem); 156 #endif 157 _store=s; 158 setsize(s); 159 } 160 changesize(size_t s)161 void changesize(size_t s){ 162 if (s <=_store){ 163 for (size_t i=0;i<_row*_col;i++) 164 _repview[i]= Polynomial(_rep.begin()+i*_store,s); 165 _size=s; 166 } 167 else { 168 std::cerr<<"ERROR: in PolynomialMatrix<polfirs> -> setsize with too large size"<<std::endl; 169 throw LinboxError("LinBox ERROR: PolynomialMatrix setsize \n"); 170 } 171 172 } setsize(size_t s)173 void setsize(size_t s){ 174 if (s <=_store){ 175 for(size_t i=0;i<_row*_col;i++){ 176 for(size_t j=_size ;j<s;j++) 177 _rep[i*_store+j]=_fld->zero; 178 } 179 for (size_t i=0;i<_row*_col;i++) 180 _repview[i]= Polynomial(_rep.begin()+i*_store,s); 181 _size=s; 182 } 183 else { 184 std::cerr<<"ERROR: in PolynomialMatrix<polfirs> -> setsize with too large size"<<std::endl; 185 throw LinboxError("LinBox ERROR: PolynomialMatrix setsize \n"); 186 } 187 } 188 189 190 // copy elt from M[beg..end], _size must be >= j-i copy(const Self_t & M,size_t beg,size_t end)191 void copy(const Self_t& M, size_t beg, size_t end){ 192 //cout<<"copying.....polfirst to polfirst.....same field"<<endl; 193 for (size_t k=0;k<_row*_col;k++){ 194 size_t j=0; 195 for (size_t i=beg;i<=end;i++){ 196 ref(k,j)=M.get(k,i); 197 j++; 198 } 199 } 200 } 201 template<typename OtherField> copy(const PolynomialMatrix<PMType::polfirst,PMStorage::plain,OtherField> & M,size_t beg,size_t end)202 void copy(const PolynomialMatrix<PMType::polfirst,PMStorage::plain,OtherField> & M, size_t beg, size_t end){ 203 //cout<<"copying.....polfirst to polfirst.....other field"<<endl; 204 Hom<OtherField,Field> hom(M.field(),field()) ; 205 for (size_t k=0;k<_row*_col;k++){ 206 size_t j=0; 207 for (size_t i=beg;i<=end;i++,j++){ 208 hom.image(ref(k,j),M.get(k,i)); 209 } 210 } 211 } 212 213 // copy elt from M[beg..end], _size must be >= end-beg+1 214 // M is stored as a Polynomial of Matrices 215 template <size_t storage> copy(const PolynomialMatrix<PMType::matfirst,storage,Field> & M,size_t beg,size_t end)216 void copy(const PolynomialMatrix<PMType::matfirst,storage,Field>& M, size_t beg, size_t end){ 217 //std::cout<<"copying.....matfirst to polfirst.....same field"<<std::endl; 218 const size_t ls = COPY_BLOCKSIZE; 219 for (size_t i = beg; i <= end; i+=ls) 220 for (size_t j = 0; j < _col * _row; j+=ls) 221 // Rk: the two loop must be interchanged in some cases 222 for (size_t _i = i; _i < std::min(end+1, i + ls); _i++) 223 for (size_t _j = j; _j < std::min(_col * _row, j + ls);++_j) 224 ref(_j,_i-beg)= M.get(_j,_i); 225 } 226 227 // copy elt from M[beg..end], _size must be >= end-beg+1 228 // M is stored as a Polynomial of Matrices with a different field 229 template<size_t storage, typename OtherField> copy(const PolynomialMatrix<PMType::matfirst,storage,OtherField> & M,size_t beg,size_t end)230 void copy(const PolynomialMatrix<PMType::matfirst,storage,OtherField> & M, size_t beg, size_t end){ 231 //std::cout<<"copying.....matfirst to polfirst.....other field"<<std::endl; 232 const size_t ls = COPY_BLOCKSIZE; 233 Hom<OtherField,Field> hom(M.field(),field()) ; 234 for (size_t i = beg; i <= end; i+=ls) 235 for (size_t j = 0; j < _col * _row; j+=ls) 236 for (size_t _i = i; _i < std::min(end+1, i + ls); _i++) { 237 for (size_t _j = j; _j < std::min(_col * _row, j + ls);++_j) 238 hom.image(ref(_j,_i-beg),M.get(_j,_i) ); 239 } 240 } 241 242 template<typename Mat> copy(const Mat & M)243 void copy(const Mat& M){ 244 linbox_check (M.size()==_size && M.rowdim()== _row && M.coldim()== _col); 245 copy(M,0,M.size()-1); 246 } 247 248 // rebind functor to change base field (e.g. apply modulo reduction) 249 template<typename _Tp1> 250 struct rebind { 251 typedef PolynomialMatrix<PMType::polfirst,PMStorage::plain, Field> Self_t; 252 typedef PolynomialMatrix<PMType::polfirst,PMStorage::plain, _Tp1> Other_t; 253 operatorrebind254 void operator() (Other_t& Ap, 255 const Self_t& A){ 256 Hom<Field, _Tp1> hom(A.field(), Ap.field()) ; 257 for (size_t i = 0; i < A._row * A._col; i++) 258 for (size_t j = 0; j < A._size; j++) 259 hom.image (Ap.ref(i,j), A.get(i,j)); 260 } 261 }; 262 263 // get access to the the k-th coeff of the ith matrix entry ref(size_t i,size_t k)264 inline Element& ref(size_t i, size_t k){ 265 return _rep[i*_store+k]; 266 } ref(size_t i,size_t j,size_t k)267 inline Element& ref(size_t i, size_t j, size_t k){ 268 return ref(i*_col+j,k); 269 } get(size_t i,size_t k)270 inline Element get(size_t i, size_t k)const { 271 return _rep[i*_store+k]; 272 } get(size_t i,size_t j,size_t k)273 inline Element get(size_t i, size_t j, size_t k) const{ 274 return get(i*_col+j,k); 275 } 276 rowdim()277 size_t rowdim() const {return _row;} coldim()278 size_t coldim() const {return _col;} degree()279 size_t degree() const {return _size-1;} size()280 size_t size() const {return _size;} storage()281 size_t storage() const {return _store;} field()282 const Field& field() const {return *_fld;} 283 dump(std::ostream & os)284 void dump(std::ostream& os) const { 285 os<<_row<<" "<<_col<<" "<<_size<<std::endl; 286 for(size_t i=0;i<_row*_col;i++) 287 for(size_t k=0;k<_size;k++) 288 os<<get(i,k)<<" "; 289 290 os<<std::endl; 291 } 292 write(std::ostream & os)293 std::ostream& write(std::ostream& os) const { return write(os,0,_size-1);} 294 write(std::ostream & os,size_t deg_min,size_t deg_max)295 std::ostream& write(std::ostream& os, size_t deg_min, size_t deg_max) const { 296 integer c; 297 int wid=-1,b; 298 field().cardinality (c); 299 300 if (c >0){ 301 wid = (int) ceil (log ((double) c) / M_LN10); 302 b= ((int) ceil (log ((double) _size) / M_LN10)); 303 wid*=10*(b*(b-1)/2.); 304 } 305 std::cout<<"Matrix(["; 306 for (size_t i = 0; i< _row;++i) { 307 os << " [ "; 308 for (size_t j = 0;j<_col;++j){ 309 os.width (wid); 310 field().write(os,get(i,j,deg_min)); 311 for (size_t k=deg_min+1;k<deg_max;++k){ 312 os<<"+"; 313 field().write(os,get(i,j,k))<<"*x^"<<k-deg_min; 314 } 315 if (deg_max-deg_min>0){ 316 os<<"+"; 317 field().write(os,get(i,j,deg_max))<<"*x^"<<deg_max-deg_min; 318 } 319 os<<(j<_col-1?",":"" ); 320 } 321 os << (i<_row-1?"],":"]" )<< std::endl; 322 } 323 std::cout<<"]);"; 324 325 return os; 326 } 327 getWritePointer()328 Element* getWritePointer(){return _rep.data();} getPointer()329 const Element* getPointer() const {return _rep.data();} 330 realmeminfo()331 size_t realmeminfo()const { 332 return _row*_col*(_store*element_storage(field())+sizeof(Polynomial));} 333 // return _rep.capacity()*sizeof(Element)+_repview.capacity()*sizeof(Polynomial);} 334 meminfo()335 size_t meminfo()const { return _rep.size()*sizeof(Element);} 336 changeField(const Field & F)337 void changeField(const Field& F){_fld=&F;} 338 339 private: 340 size_t _store; 341 std::vector<Polynomial> _repview; 342 //std::vector<Element> _rep; 343 VECT _rep; 344 size_t _row; 345 size_t _col; 346 size_t _size; 347 const Field* _fld; 348 }; 349 350 // Class for Polynomial Matrix stored as a Polynomial of Matrices 351 template<class _Field> 352 class PolynomialMatrix<PMType::matfirst,PMStorage::plain,_Field> { 353 public: 354 typedef _Field Field; 355 typedef typename Field::Element Element; 356 typedef BlasMatrix<Field> Matrix; 357 typedef std::vector<Element> Polynomial; 358 typedef PolynomialMatrix<PMType::matfirst,PMStorage::plain,_Field> Self_t; 359 typedef PolynomialMatrix<PMType::polfirst,PMStorage::plain,_Field> Other_t; 360 361 typedef Matrix value_type; 362 PolynomialMatrix()363 PolynomialMatrix() {} 364 365 PolynomialMatrix(const Self_t&) = delete; 366 PolynomialMatrix(const Field & f,size_t r,size_t c,size_t s)367 PolynomialMatrix(const Field& f, size_t r, size_t c, size_t s) : 368 _rep(s,Matrix(f)), _row(r), _col(c), _size(s), _fld(&f) { 369 //_row(r), _col(c), _size(s), _fld(&f) { 370 for(size_t i=0;i<s;i++) 371 _rep[i].init(f,r,c); 372 //integer p; 373 //std::cout<<"(ALLOC) matfirst at "<<this<<" : "<<r<<"x"<<c<<" - size= "<<s<<" ==> "<<MB(realmeminfo())<<" Mo"<<std::endl; 374 ADD_MEM(realmeminfo()); 375 } 376 ~PolynomialMatrix()377 ~PolynomialMatrix(){ 378 DEL_MEM(realmeminfo()); 379 //std::cout<<"(FREE) matfirst at "<<this<<" : "<<_row<<"x"<<_col<<" - size= "<<_size<<" ==> "<<MB(realmeminfo())<<" Mo"<<std::endl; 380 //integer p; 381 //std::cout<<"PMatrix Desallocating : "<<_row*_col*_size*length(_fld->characteristic(p))/1000000.<<"Mo"<<std::endl; 382 } 383 clear()384 void clear(){ 385 _rep.resize(0,Matrix(*_fld)); 386 } 387 388 389 // retrieve the matrix of degree k in the polynomial matrix 390 inline Matrix& operator[](size_t k) {return _rep[k];} 391 inline const Matrix& operator[](size_t k)const {return _rep[k];} 392 393 // retrieve the polynomial at entry (i,j) in the matrix operator()394 inline Polynomial operator()(size_t i, size_t j) const { 395 Polynomial A(_size); 396 for(size_t k=0;k<_size;k++) 397 A[k]=_rep[k].getEntry(i,j); 398 return A; 399 } 400 // retrieve the polynomial at psotion (i) in the storage of the matrix operator()401 inline Polynomial operator()(size_t i) const { 402 Polynomial A(_size); 403 for(size_t k=0;k<_size;k++) 404 A[k]=*(_rep[k].Begin()+i); 405 return A; 406 } 407 408 // resize the polynomial length of the polynomial matrix resize(size_t s)409 void resize(size_t s){ 410 _rep.resize(s,Matrix(*_fld)); 411 if (s>_size){ 412 for(size_t i=_size;i<s;i++) 413 _rep[i].init(field(),_row,_col); 414 } 415 _size=s; 416 } 417 setsize(size_t s)418 void setsize(size_t s){resize(s);} 419 420 // copy elt from M[beg..end], _size must be >= j-i 421 template <size_t storage> 422 void copy(const PolynomialMatrix<PMType::matfirst,storage,Field>& M, size_t beg, size_t end, size_t start=0){ 423 //cout<<"copying.....matfirst to matfirst.....same field"<<endl; 424 for (size_t i=beg;i<=end;i++) 425 _rep[start+i-beg]=M[i]; 426 427 } 428 429 template<size_t storage, typename OtherField> 430 void copy(const PolynomialMatrix<PMType::matfirst,storage,OtherField> & M, size_t beg, size_t end, size_t start=0){ 431 //cout<<"copying.....matfirst to matfirst.....other field"<<endl; 432 Hom<OtherField,Field> hom(M.field(),field()) ; 433 for (size_t i=beg;i<=end;i++) 434 for (size_t j=0;j<_row*_col;j++) 435 hom.image(ref(j,start+i-beg) , M.get(j,i)); 436 } 437 438 // copy elt from M[beg..end], _size must be >= end-beg+1 439 // M is stored as a Matrix of Polynomials 440 void copy(const Other_t& M, size_t beg, size_t end, size_t start=0){ 441 //cout<<"copying.....polfirst to matfirst.....same field"<<endl; 442 const size_t ls = COPY_BLOCKSIZE; // Loop unrooling to be more cache-friendly (cache block transposition) 443 for(size_t i = beg; i < end+1; i+=ls) 444 for (size_t j = 0; j < _col * _row; j+=ls) 445 for (size_t _i = i; _i < std::min(end+1, i + ls); _i++) { 446 for (size_t _j = j; _j < std::min(_col * _row, j + ls); _j++) 447 ref(_j,start+_i-beg) = M.get(_j,_i); 448 } 449 } 450 451 // copy elt from M[beg..end], _size must be >= end-beg+1 452 // M is stored as a Matrix of Polynomials with a different field 453 template<typename OtherField> 454 void copy(const PolynomialMatrix<PMType::polfirst,PMStorage::plain,OtherField> & M, size_t beg, size_t end, size_t start=0){ 455 //cout<<"copying.....polfirst to matfirst.....other field"<<endl; 456 const size_t ls = COPY_BLOCKSIZE; 457 Hom<OtherField,Field> hom(M.field(),field()) ; 458 for(size_t i = beg; i < end+1; i+=ls) 459 for (size_t j = 0; j < _col * _row; j+=ls) 460 for (size_t _i = i; _i < std::min(end+1, i + ls); _i++) { 461 for (size_t _j = j; _j < std::min(_col * _row, j + ls); _j++) 462 hom.image(ref(_j,start+_i-beg) , M.get(_j,_i)); 463 } 464 } 465 466 template<typename Mat> copy(const Mat & M)467 void copy(const Mat& M){ 468 linbox_check (M.size()==_size && M.rowdim()== _row && M.coldim()== _col); 469 copy(M,0,M.size()-1); 470 } 471 472 // rebind functor to change base field (e.g. apply modulo reduction) 473 template<typename _Tp1> 474 struct rebind { 475 typedef PolynomialMatrix<PMType::matfirst, PMStorage::plain,Field> Self_t; 476 typedef PolynomialMatrix<PMType::matfirst, PMStorage::plain,_Tp1> Other_t; 477 478 template<size_t storage> operatorrebind479 void operator() (PolynomialMatrix<PMType::matfirst, storage,_Tp1>& Ap, 480 const PolynomialMatrix<PMType::matfirst, storage,Field>& A){ 481 Hom<Field, _Tp1> hom(A.field(), Ap.field()) ; 482 for (size_t j = 0; j < A._size; j++) 483 for (size_t i = 0; i < A._row * A._col; i++) 484 hom.image (Ap.ref(i,j), A.get(i,j)); 485 } 486 }; 487 488 typedef PolynomialMatrix<PMType::matfirst,PMStorage::plain,Field> plain; 489 typedef PolynomialMatrix<PMType::matfirst,PMStorage::view,Field> view; 490 typedef PolynomialMatrix<PMType::matfirst,PMStorage::const_view,Field> const_view; 491 at(size_t i,size_t j)492 inline view at(size_t i, size_t j) {return view(*this,i,j);} at(size_t i,size_t j)493 inline const_view at(size_t i, size_t j)const {return const_view(*this,i,j);} 494 495 496 // get access to the the k-th coeff of the ith matrix entry ref(size_t i,size_t k)497 inline Element& ref(size_t i, size_t k){ 498 return _rep[k].refRep()[i]; 499 //return *(_rep[k].Begin()+i); 500 } ref(size_t i,size_t j,size_t k)501 inline Element& ref(size_t i, size_t j, size_t k){ 502 return ref(i*_col+j,k); 503 } get(size_t i,size_t k)504 inline Element get(size_t i, size_t k)const { 505 return *(_rep[k].Begin()+i); 506 } get(size_t i,size_t j,size_t k)507 inline Element get(size_t i, size_t j, size_t k) const{ 508 return get(i*_col+j,k); 509 } 510 rowdim()511 inline size_t rowdim() const {return _row;} coldim()512 inline size_t coldim() const {return _col;} degree()513 inline size_t degree() const {return _size-1;} real_degree()514 inline size_t real_degree() const { 515 MatrixDomain<Field> MD(field()); 516 size_t d= _size-1; 517 while(d>0 && MD.isZero(_rep[d])) d--; 518 return d; 519 } size()520 inline size_t size() const {return _size;} storage()521 inline size_t storage()const {return _size;} 522 field()523 inline const Field& field() const {return *_fld;} 524 dump(std::ostream & os)525 void dump(std::ostream& os) const { 526 os<<_row<<" "<<_col<<" "<<_size<<std::endl; 527 for(size_t k=0;k<_size;k++) 528 for(size_t i=0;i<_row*_col;i++) 529 os<<get(i,k)<<" "; 530 os<<std::endl; 531 } 532 533 write(std::ostream & os)534 std::ostream& write(std::ostream& os) const { return write(os,0,real_degree());} 535 write(std::ostream & os,size_t deg_min,size_t deg_max)536 std::ostream& write(std::ostream& os, size_t deg_min, size_t deg_max) const { 537 integer c; 538 int wid,b; 539 field().cardinality (c); 540 541 if (c >0) 542 wid = (int) ceil (log ((double) c) / M_LN10); 543 else 544 wid=(int) ceil (log ((double) _rep[0].getEntry(0,0)) / M_LN10); 545 546 b= ((int) ceil (log ((double) (deg_max-deg_min+1)) / M_LN10)); 547 wid*=10*(b*(b-1)/2.); 548 std::cout<<"Matrix(["; 549 for (size_t i = 0; i< _row;++i) { 550 os << " [ "; 551 for (size_t j = 0;j<_col;++j){ 552 os.width (wid); 553 field().write(os,_rep[deg_min].getEntry(i,j)); 554 for (size_t k=deg_min+1;k<deg_max;++k){ 555 os<<"+"; 556 field().write(os,_rep[k].getEntry(i,j))<<"*x^"<<k-deg_min; 557 } 558 if (deg_max-deg_min>0){ 559 os<<"+"; 560 field().write(os,_rep[deg_max].getEntry(i,j))<<"*x^"<<deg_max-deg_min; 561 } 562 os<<(j<_col-1?",":"" ); 563 } 564 os << (i<_row-1?"],":"]" )<< std::endl; 565 } 566 std::cout<<"]);"; 567 return os; 568 } 569 realmeminfo()570 size_t realmeminfo()const { 571 572 return _size*(_row*_col*element_storage(field())+sizeof(Matrix)); 573 } 574 575 // NEED FOR YUHASZ 576 typedef typename std::vector<Matrix>::const_iterator const_iterator; begin()577 const_iterator begin() const {return _rep.begin();} 578 579 private: 580 std::vector<Matrix> _rep; 581 size_t _row; 582 size_t _col; 583 size_t _size; 584 const Field* _fld; 585 }; 586 587 template<typename _Field, size_t T, size_t S> 588 std::ostream& operator<<(std::ostream& os, const PolynomialMatrix<T,S,_Field>& P) { 589 return P.write(os); 590 } 591 592 // Class to handle the view of a Polynomial Matrix according to some degree range 593 // the view is read/write 594 template<class _Field> 595 class PolynomialMatrix<PMType::matfirst,PMStorage::view,_Field> { 596 public: 597 typedef _Field Field; 598 typedef typename Field::Element Element; 599 typedef BlasMatrix<Field> Matrix; 600 typedef std::vector<Element> Polynomial; 601 PolynomialMatrix()602 PolynomialMatrix() {} 603 604 // constructor of a view between i and j from a plain Polynomial Matrix PolynomialMatrix(PolynomialMatrix<PMType::matfirst,PMStorage::plain,Field> & M,size_t i,size_t j)605 PolynomialMatrix(PolynomialMatrix<PMType::matfirst,PMStorage::plain,Field>& M, size_t i,size_t j) 606 : _ptr(&M), _size(j-i+1), _shift(i) 607 {//cout<<i<<"-"<<j<<" -> "<<M.size()<<endl; 608 linbox_check(i<M.size() && i<=j && j< M.size());} 609 610 // constructor of a view between i and j from a view Polynomial Matrix PolynomialMatrix(PolynomialMatrix<PMType::matfirst,PMStorage::view,Field> & M,size_t i,size_t j)611 PolynomialMatrix(PolynomialMatrix<PMType::matfirst,PMStorage::view,Field>& M, size_t i,size_t j) 612 : _ptr(M._ptr), _size(j-i+1), _shift(i+M._shift) 613 {//cout<<i<<"-"<<j<<" -> "<<M.size()<<endl; 614 linbox_check(i<M.size() && i<=j && j< M.size());} 615 616 617 // retrieve the matrix of degree k in the polynomial matrix 618 inline Matrix& operator[](size_t k) {return (*_ptr)[k+_shift];} 619 inline const Matrix& operator[](size_t k)const {return (*_ptr)[k+_shift];} 620 621 // retrieve the polynomial at entry (i,j) in the matrix operator()622 Polynomial operator()(size_t i, size_t j){ 623 Polynomial A(_size, Matrix(_ptr->field())); 624 for(size_t k=0;k<_size;k++) 625 A[k]=(*_ptr)[k+_shift].refEntry(i,j); 626 return A; 627 } 628 629 // get access to the the k-th coeff of the ith matrix entry ref(size_t i,size_t k)630 inline Element& ref(size_t i, size_t k){ 631 return _ptr->ref(i,k+_shift); 632 } ref(size_t i,size_t j,size_t k)633 inline Element& ref(size_t i, size_t j, size_t k){ 634 return ref(i*coldim()+j,k); 635 } get(size_t i,size_t k)636 inline Element get(size_t i, size_t k)const { 637 return _ptr->get(i,k+_shift); 638 } get(size_t i,size_t j,size_t k)639 inline Element get(size_t i, size_t j, size_t k) const{ 640 return get(i*coldim()+j,k); 641 } 642 rowdim()643 inline size_t rowdim() const {return _ptr->rowdim();} coldim()644 inline size_t coldim() const {return _ptr->coldim();} degree()645 inline size_t degree() const {return _size-1;} size()646 inline size_t size() const {return _size;} field()647 inline const Field& field() const {return _ptr->field();} 648 649 typedef PolynomialMatrix<PMType::matfirst,PMStorage::plain,Field> plain; 650 typedef PolynomialMatrix<PMType::matfirst,PMStorage::view,Field> view; 651 typedef PolynomialMatrix<PMType::matfirst,PMStorage::const_view,Field> const_view; 652 653 template<typename Mat> copy(const Mat & M)654 void copy(const Mat& M){ 655 this->copy(M,0,M.size()-1); 656 } 657 658 template<typename Mat> copy(const Mat & M,size_t beg,size_t end)659 void copy(const Mat& M, size_t beg, size_t end){_ptr->copy(M,beg,end,_shift);} 660 661 at(size_t i,size_t j)662 inline view at(size_t i, size_t j) {return view(*this,i,j);} at(size_t i,size_t j)663 inline const_view at(size_t i, size_t j)const {return const_view(*this,i,j);} 664 write(std::ostream & os)665 std::ostream& write(std::ostream& os) const { return _ptr->write(os,_shift,_shift+_size-1);} 666 667 friend class PolynomialMatrix<PMType::matfirst,PMStorage::const_view,Field>; 668 669 private: 670 PolynomialMatrix<PMType::matfirst,PMStorage::plain,Field>* _ptr; 671 size_t _size; 672 size_t _shift; 673 }; 674 675 // Class to handle the view of a Polynomial Matrix according to some degree range 676 // the view is read only 677 template<class _Field> 678 class PolynomialMatrix<PMType::matfirst,PMStorage::const_view,_Field> { 679 public: 680 typedef _Field Field; 681 typedef typename Field::Element Element; 682 typedef BlasMatrix<Field> Matrix; 683 typedef std::vector<Element> Polynomial; 684 PolynomialMatrix()685 PolynomialMatrix() {} 686 687 // constructor of a view between i and j from a plain Polynomial Matrix PolynomialMatrix(const PolynomialMatrix<PMType::matfirst,PMStorage::plain,Field> & M,size_t i,size_t j)688 PolynomialMatrix(const PolynomialMatrix<PMType::matfirst,PMStorage::plain,Field>& M, size_t i,size_t j) 689 : _ptr(&M), _size(j-i+1), _shift(i) 690 {linbox_check(i<M.size() && i<=j && j< M.size());} 691 692 // constructor of a view between i and j from a view Polynomial Matrix PolynomialMatrix(const PolynomialMatrix<PMType::matfirst,PMStorage::view,Field> & M,size_t i,size_t j)693 PolynomialMatrix(const PolynomialMatrix<PMType::matfirst,PMStorage::view,Field>& M, size_t i,size_t j) 694 : _ptr(M._ptr), _size(j-i+1), _shift(i+M._shift) 695 {linbox_check(i<M.size() && i<=j && j< M.size());} 696 697 698 // constructor of a view between i and j from a view Polynomial Matrix PolynomialMatrix(const PolynomialMatrix<PMType::matfirst,PMStorage::const_view,Field> & M,size_t i,size_t j)699 PolynomialMatrix(const PolynomialMatrix<PMType::matfirst,PMStorage::const_view,Field>& M, size_t i,size_t j) 700 : _ptr(M._ptr), _size(j-i+1), _shift(i+M._shift) 701 {//cout<<"("<<i<<","<<j<<") -> "<<M.size()<<endl; 702 linbox_check(i<M.size()&& i<=j && j<= M.size());} 703 704 705 // retrieve the matrix of degree k in the polynomial matrix 706 inline const Matrix& operator[](size_t k)const {return (*_ptr)[k+_shift];} 707 708 // retrieve the polynomial at entry (i,j) in the matrix operator()709 inline Polynomial operator()(size_t i, size_t j){ 710 Polynomial A(_size, Matrix(_ptr->field())); 711 for(size_t k=0;k<_size;k++) 712 A[k]=(*_ptr)[k+_shift].refEntry(i,j); 713 return A; 714 } 715 get(size_t i,size_t k)716 inline Element get(size_t i, size_t k) const { 717 return _ptr->get(i,k+_shift); 718 } get(size_t i,size_t j,size_t k)719 inline Element get(size_t i, size_t j, size_t k) const{ 720 return get(i*coldim()+j,k); 721 } 722 723 rowdim()724 inline size_t rowdim() const {return _ptr->rowdim();} coldim()725 inline size_t coldim() const {return _ptr->coldim();} degree()726 inline size_t degree() const {return _size-1;} size()727 inline size_t size() const {return _size;} field()728 inline const Field& field() const {return _ptr->field();} 729 730 typedef PolynomialMatrix<PMType::matfirst,PMStorage::plain,Field> plain; 731 typedef PolynomialMatrix<PMType::matfirst,PMStorage::const_view,Field> const_view; 732 at(size_t i,size_t j)733 inline const_view at(size_t i, size_t j)const {return const_view(*this,i,j);} write(std::ostream & os)734 std::ostream& write(std::ostream& os) const { return _ptr->write(os,_shift,_shift+_size-1);} 735 736 private: 737 const PolynomialMatrix<PMType::matfirst,PMStorage::plain,Field>* _ptr; 738 size_t _size; 739 size_t _shift; 740 }; 741 742 743 744 745 746 } //end of namespace LinBox 747 748 #endif // __LINBOX_polynomial_matrix_H 749 750 // Local Variables: 751 // mode: C++ 752 // tab-width: 4 753 // indent-tabs-mode: nil 754 // c-basic-offset: 4 755 // End: 756 // vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 757