[def863] | 1 | /* |
---|
| 2 | * lib_zmatrix.h |
---|
| 3 | * |
---|
| 4 | * Created on: Sep 28, 2010 |
---|
| 5 | * Author: anders |
---|
| 6 | */ |
---|
| 7 | |
---|
| 8 | #ifndef LIB_ZMATRIX_H_ |
---|
| 9 | #define LIB_ZMATRIX_H_ |
---|
| 10 | |
---|
| 11 | #include <vector> |
---|
| 12 | #include <algorithm> |
---|
| 13 | #include "gfanlib_vector.h" |
---|
| 14 | |
---|
| 15 | namespace gfan{ |
---|
| 16 | |
---|
| 17 | template <class typ> class Matrix{ |
---|
| 18 | int width,height; |
---|
| 19 | std::vector<Vector<typ> > rows; |
---|
| 20 | public: |
---|
| 21 | inline int getHeight()const{return height;}; |
---|
| 22 | inline int getWidth()const{return width;}; |
---|
[26b713] | 23 | Matrix(const Matrix &a):width(a.getWidth()),height(a.getHeight()),rows(a.rows){ |
---|
[def863] | 24 | } |
---|
[26b713] | 25 | Matrix(int height_, int width_):width(width_),height(height_),rows(height_){ |
---|
[def863] | 26 | assert(height>=0); |
---|
| 27 | assert(width>=0); |
---|
| 28 | for(int i=0;i<getHeight();i++)rows[i]=Vector<typ>(width); |
---|
| 29 | }; |
---|
| 30 | ~Matrix(){ |
---|
| 31 | }; |
---|
| 32 | Matrix():width(0),height(0){ |
---|
| 33 | }; |
---|
| 34 | static Matrix rowVectorMatrix(Vector<typ> const &v) |
---|
| 35 | { |
---|
| 36 | Matrix ret(1,v.size()); |
---|
| 37 | for(int i=0;i<v.size();i++)ret[0][i]=v[i]; |
---|
| 38 | return ret; |
---|
| 39 | } |
---|
| 40 | Vector<typ> column(int i)const |
---|
| 41 | { |
---|
| 42 | assert(i>=0); |
---|
| 43 | assert(i<getWidth()); |
---|
| 44 | Vector<typ> ret(getHeight()); |
---|
| 45 | for(int j=0;j<getHeight();j++)ret[j]=rows[j][i]; |
---|
| 46 | return ret; |
---|
| 47 | } |
---|
| 48 | Matrix transposed()const |
---|
| 49 | { |
---|
| 50 | Matrix ret(getWidth(),getHeight()); |
---|
| 51 | for(int i=0;i<getWidth();i++) |
---|
| 52 | ret.rows[i]=column(i); |
---|
| 53 | return ret; |
---|
| 54 | } |
---|
| 55 | static Matrix identity(int n) |
---|
| 56 | { |
---|
| 57 | Matrix m(n,n); |
---|
| 58 | for(int i=0;i<n;i++)m.rows[i]=Vector<typ>::standardVector(n,i); |
---|
| 59 | return m; |
---|
| 60 | } |
---|
| 61 | void append(Matrix const &m) |
---|
| 62 | { |
---|
| 63 | for(int i=0;i<m.height;i++) |
---|
| 64 | { |
---|
| 65 | rows.push_back(m[i]); |
---|
| 66 | } |
---|
| 67 | height+=m.height; |
---|
| 68 | } |
---|
| 69 | void appendRow(Vector<typ> const &v) |
---|
| 70 | { |
---|
| 71 | assert(v.size()==width); |
---|
| 72 | rows.push_back(v); |
---|
| 73 | height++; |
---|
| 74 | } |
---|
| 75 | void eraseLastRow() |
---|
| 76 | { |
---|
| 77 | assert(rows.size()>0); |
---|
| 78 | rows.resize(rows.size()-1); |
---|
| 79 | height--; |
---|
| 80 | } |
---|
| 81 | /*IntegerVector vectormultiply(IntegerVector const &v)const |
---|
| 82 | { |
---|
| 83 | assert(v.size()==width); |
---|
| 84 | IntegerVector ret(height); |
---|
| 85 | for(int i=0;i<height;i++) |
---|
| 86 | ret[i]=dot(rows[i],v); |
---|
| 87 | return ret; |
---|
| 88 | }*/ |
---|
| 89 | /** |
---|
| 90 | * Decides if v is in the kernel of the matrix. |
---|
| 91 | */ |
---|
| 92 | /* bool inKernel(IntegerVector const &v)const |
---|
| 93 | { |
---|
| 94 | assert(v.size()==width); |
---|
| 95 | for(int i=0;i<height;i++) |
---|
| 96 | if(dotLong(rows[i],v)!=0)return false; |
---|
| 97 | return true; |
---|
| 98 | } |
---|
| 99 | */ |
---|
| 100 | friend Matrix operator*(const typ &s, const Matrix& q) |
---|
| 101 | { |
---|
| 102 | Matrix p=q; |
---|
| 103 | for(int i=0;i<q.height;i++)p[i]=s*(q[i]); |
---|
| 104 | return p; |
---|
| 105 | } |
---|
[26b713] | 106 | friend Matrix operator*(const Matrix& a, const Matrix& b) |
---|
[def863] | 107 | { |
---|
| 108 | assert(a.width==b.height); |
---|
| 109 | Matrix ret(b.width,a.height); |
---|
| 110 | for(int i=0;i<b.width;i++) |
---|
| 111 | ret[i]=a.vectormultiply(b.column(i)); |
---|
| 112 | return ret.transposed(); |
---|
[26b713] | 113 | } |
---|
[def863] | 114 | /* template<class T> |
---|
| 115 | Matrix<T>(const Matrix<T>& c):v(c.size()){ |
---|
| 116 | for(int i=0;i<size();i++)v[i]=typ(c[i]);} |
---|
| 117 | */ |
---|
| 118 | friend Matrix operator-(const Matrix &b) |
---|
| 119 | { |
---|
| 120 | Matrix ret(b.height,b.width); |
---|
| 121 | for(int i=0;i<b.height;i++)ret[i]=-b[i]; |
---|
| 122 | return ret; |
---|
| 123 | } |
---|
| 124 | |
---|
| 125 | /** |
---|
| 126 | Returns the specified submatrix. The endRow and endColumn are not included. |
---|
| 127 | */ |
---|
| 128 | Matrix submatrix(int startRow, int startColumn, int endRow, int endColumn)const |
---|
| 129 | { |
---|
| 130 | assert(startRow>=0); |
---|
| 131 | assert(startColumn>=0); |
---|
| 132 | assert(endRow>=startRow); |
---|
| 133 | assert(endColumn>=startColumn); |
---|
| 134 | assert(endRow<=height); |
---|
| 135 | assert(endColumn<=width); |
---|
| 136 | Matrix ret(endRow-startRow,endColumn-startColumn); |
---|
| 137 | for(int i=startRow;i<endRow;i++) |
---|
| 138 | for(int j=startColumn;j<endColumn;j++) |
---|
| 139 | ret[i-startRow][j-startColumn]=rows[i][j]; |
---|
| 140 | return ret; |
---|
| 141 | } |
---|
| 142 | const Vector<typ>& operator[](int n)const{assert(n>=0 && n<getHeight());return (rows[n]);} |
---|
| 143 | // Bugfix for gcc4.5 (passing assertion to the above operator): |
---|
| 144 | Vector<typ>& operator[](int n){if(!(n>=0 && n<getHeight())){(*(const Matrix<typ>*)(this))[n];}return (rows[n]);} |
---|
| 145 | |
---|
| 146 | |
---|
| 147 | bool operator<(const Matrix & b)const |
---|
| 148 | { |
---|
| 149 | if(getWidth()<b.getWidth())return true; |
---|
| 150 | if(b.getWidth()<getWidth())return false; |
---|
| 151 | if(getHeight()<b.getHeight())return true; |
---|
| 152 | if(b.getHeight()<getHeight())return false; |
---|
| 153 | |
---|
| 154 | for(int i=0;i<getHeight();i++) |
---|
| 155 | { |
---|
| 156 | if((*this)[i]<b[i])return true; |
---|
| 157 | if(b[i]<(*this)[i])return false; |
---|
| 158 | } |
---|
| 159 | return false; |
---|
| 160 | } |
---|
| 161 | /** |
---|
| 162 | Adds a times the i th row to the j th row. |
---|
| 163 | */ |
---|
| 164 | void madd(int i, typ a, int j) |
---|
| 165 | { |
---|
| 166 | assert(i!=j); |
---|
| 167 | assert(i>=0 && i<height); |
---|
| 168 | assert(j>=0 && j<height); |
---|
| 169 | |
---|
| 170 | if(!a.isZero()) |
---|
| 171 | for(int k=0;k<width;k++) |
---|
| 172 | if(!rows[i][k].isZero()) |
---|
| 173 | rows[j][k].madd(rows[i][k],a); |
---|
| 174 | } |
---|
| 175 | |
---|
| 176 | friend std::ostream &operator<<(std::ostream &f, Matrix const &a){ |
---|
| 177 | f<<"{"; |
---|
| 178 | for(int i=0;i<a.getHeight();i++) |
---|
| 179 | { |
---|
| 180 | if(i)f<<","<<std::endl; |
---|
| 181 | f<<a.rows[i]; |
---|
| 182 | } |
---|
| 183 | f<<"}"<<std::endl; |
---|
| 184 | return f; |
---|
| 185 | } |
---|
| 186 | /** |
---|
| 187 | Swaps the i th and the j th row. |
---|
| 188 | */ |
---|
| 189 | void swapRows(int i, int j) |
---|
| 190 | { |
---|
| 191 | std::swap(rows[i],rows[j]); |
---|
| 192 | } |
---|
| 193 | /** |
---|
| 194 | This method is used for iterating through the pivots in a matrix |
---|
| 195 | in row echelon form. To find the first pivot put i=-1 and |
---|
| 196 | j=-1 and call this routine. The (i,j) th entry of the matrix is a |
---|
| 197 | pivot. Call the routine again to get the next pivot. When no more |
---|
| 198 | pivots are found the routine returns false. |
---|
| 199 | */ |
---|
| 200 | bool nextPivot(int &i, int &j)const |
---|
| 201 | { |
---|
| 202 | i++; |
---|
| 203 | if(i>=height)return false; |
---|
| 204 | while(++j<width) |
---|
| 205 | { |
---|
| 206 | if(!rows[i][j].isZero()) return true; |
---|
| 207 | } |
---|
| 208 | return false; |
---|
| 209 | } |
---|
| 210 | /** |
---|
| 211 | Returns the indices of the columns containing a pivot. |
---|
| 212 | The returned list is sorted. |
---|
| 213 | The matrix must be in row echelon form. |
---|
| 214 | */ |
---|
| 215 | std::vector<int> pivotColumns()const |
---|
| 216 | { |
---|
| 217 | std::vector<int> ret; |
---|
| 218 | int pivotI=-1; |
---|
| 219 | int pivotJ=-1; |
---|
| 220 | while(nextPivot(pivotI,pivotJ))ret.push_back(pivotJ); |
---|
| 221 | return ret; |
---|
| 222 | } |
---|
| 223 | /** |
---|
| 224 | Returns the indices of the columns not containing a pivot. |
---|
| 225 | The returned list is sorted. |
---|
| 226 | The matrix must be in row echelon form. |
---|
| 227 | */ |
---|
| 228 | std::vector<int> nonPivotColumns()const |
---|
| 229 | { |
---|
| 230 | std::vector<int> ret; |
---|
| 231 | int pivotI=-1; |
---|
| 232 | int pivotJ=-1; |
---|
| 233 | int firstPossiblePivot=0; |
---|
| 234 | while(nextPivot(pivotI,pivotJ)) |
---|
| 235 | { |
---|
| 236 | for(int j=firstPossiblePivot;j<pivotJ;j++) |
---|
| 237 | ret.push_back(j); |
---|
| 238 | firstPossiblePivot=pivotJ+1; |
---|
| 239 | } |
---|
| 240 | for(int j=firstPossiblePivot;j<getWidth();j++) |
---|
| 241 | ret.push_back(j); |
---|
| 242 | |
---|
| 243 | return ret; |
---|
| 244 | } |
---|
| 245 | /** |
---|
| 246 | This routine removes the zero rows of the matrix. |
---|
| 247 | */ |
---|
| 248 | void removeZeroRows() |
---|
| 249 | { |
---|
| 250 | int nonZeros=0; |
---|
| 251 | for(int i=0;i<height;i++)if(!(*this)[i].isZero())nonZeros++; |
---|
| 252 | if(nonZeros==height)return; |
---|
| 253 | |
---|
| 254 | Matrix b(nonZeros,width); |
---|
| 255 | |
---|
| 256 | int j=0; |
---|
| 257 | for(int i=0;i<height;i++) |
---|
| 258 | { |
---|
| 259 | if(!(*this)[i].isZero()) |
---|
| 260 | { |
---|
| 261 | b[j]=(*this)[i]; |
---|
| 262 | j++; |
---|
| 263 | } |
---|
| 264 | } |
---|
| 265 | *this=b; |
---|
| 266 | } |
---|
| 267 | /** |
---|
| 268 | Returns the index of a row whose index is at least currentRow and |
---|
| 269 | which has a non-zero element on the column th entry. If no such |
---|
| 270 | row exists then -1 is returned. This routine is used in the Gauss |
---|
| 271 | reduction. To make the reduction more efficient the routine |
---|
| 272 | chooses its row with as few non-zero entries as possible. |
---|
| 273 | */ |
---|
| 274 | int findRowIndex(int column, int currentRow)const |
---|
| 275 | { |
---|
| 276 | int best=-1; |
---|
| 277 | int bestNumberOfNonZero=0; |
---|
| 278 | for(int i=currentRow;i<height;i++) |
---|
| 279 | if(!rows[i][column].isZero()) |
---|
| 280 | { |
---|
| 281 | int nz=0; |
---|
| 282 | for(int k=column+1;k<width;k++) |
---|
| 283 | if(!rows[i][k].isZero())nz++; |
---|
| 284 | if(best==-1) |
---|
| 285 | { |
---|
| 286 | best=i; |
---|
| 287 | bestNumberOfNonZero=nz; |
---|
| 288 | } |
---|
| 289 | else if(nz<bestNumberOfNonZero) |
---|
| 290 | { |
---|
| 291 | best=i; |
---|
| 292 | bestNumberOfNonZero=nz; |
---|
| 293 | } |
---|
| 294 | } |
---|
| 295 | return best; |
---|
| 296 | } |
---|
| 297 | /** |
---|
| 298 | Performs a Gauss reduction and returns the number of row swaps |
---|
| 299 | done. The result is a matrix in row echelon form. The pivots may |
---|
| 300 | not be all 1. In terms of Groebner bases, what is computed is a |
---|
| 301 | minimal (not necessarily reduced) Groebner basis of the linear |
---|
| 302 | ideal generated by the rows. The number of swaps is need if one |
---|
| 303 | wants to compute the determinant afterwards. In this case it is |
---|
| 304 | also a good idea to set the flag returnIfZeroDeterminant which |
---|
| 305 | make the routine terminate before completion if it discovers that |
---|
| 306 | the determinant is zero. |
---|
| 307 | */ |
---|
[26b713] | 308 | int reduce(bool returnIfZeroDeterminant=false, bool integral=false, bool makePivotsOne=false) |
---|
[def863] | 309 | { |
---|
| 310 | assert(integral || typ::isField()); |
---|
| 311 | int retSwaps=0; |
---|
| 312 | int currentRow=0; |
---|
| 313 | |
---|
| 314 | for(int i=0;i<width;i++) |
---|
| 315 | { |
---|
| 316 | int s=findRowIndex(i,currentRow); |
---|
| 317 | |
---|
| 318 | if(s!=-1) |
---|
| 319 | { |
---|
| 320 | if(s!=currentRow) |
---|
| 321 | { |
---|
| 322 | swapRows(currentRow,s); |
---|
| 323 | retSwaps++; |
---|
| 324 | } |
---|
[26b713] | 325 | if(makePivotsOne) |
---|
| 326 | {//THE PIVOT SHOULD BE SET TO ONE IF INTEGRAL IS FALSE |
---|
| 327 | if(!rows[currentRow][i].sign()<0)retSwaps++; |
---|
| 328 | typ inverse=typ(1)/rows[currentRow][i]; |
---|
| 329 | // if(!rows[currentRow][i].isOne()) |
---|
| 330 | { |
---|
| 331 | for(int k=0;k<width;k++) |
---|
| 332 | if(!rows[currentRow][k].isZero()) |
---|
| 333 | rows[currentRow][k]*=inverse; |
---|
| 334 | } |
---|
| 335 | } |
---|
[def863] | 336 | for(int j=currentRow+1;j<height;j++) |
---|
| 337 | if(integral) |
---|
| 338 | { |
---|
| 339 | if(!rows[j][i].isZero()) |
---|
| 340 | { |
---|
| 341 | typ s; |
---|
| 342 | typ t; |
---|
| 343 | |
---|
| 344 | typ g=typ::gcd(rows[currentRow][i],rows[j][i],s,t); |
---|
| 345 | typ u=-rows[j][i]/g; |
---|
| 346 | typ v=rows[currentRow][i]/g; |
---|
| 347 | /* We want the (s,t) vector to be as small as possible. |
---|
| 348 | * We are allowed to adjust by multiples of (u,v). |
---|
| 349 | * The following computes the correct multiplier (in most cases). |
---|
| 350 | */ |
---|
| 351 | /* { |
---|
| 352 | FieldElement multiplier=(s*u+t*v)*((u*u+v*v).inverse()); |
---|
| 353 | double d=mpq_get_d(*(multiplier.getGmpRationalTemporaryPointer())); |
---|
| 354 | multiplier=multiplier.getField()->zHomomorphism(-(((int)(d+10000.5))-10000)); |
---|
| 355 | s.madd(multiplier,u); |
---|
| 356 | t.madd(multiplier,v); |
---|
| 357 | }*/ |
---|
| 358 | for(int k=0;k<width;k++) |
---|
| 359 | { |
---|
| 360 | typ A=rows[currentRow][k]; |
---|
| 361 | typ B=rows[j][k]; |
---|
| 362 | |
---|
| 363 | rows[currentRow][k]=s*A+t*B; |
---|
| 364 | rows[j][k]=u*A+v*B; |
---|
| 365 | } |
---|
| 366 | } |
---|
| 367 | } |
---|
| 368 | else |
---|
| 369 | { |
---|
| 370 | if(!rows[j][i].isZero()) |
---|
| 371 | madd(currentRow,-rows[j][i]/rows[currentRow][i],j); |
---|
| 372 | } |
---|
| 373 | currentRow++; |
---|
| 374 | } |
---|
| 375 | else |
---|
| 376 | if(returnIfZeroDeterminant)return -1; |
---|
| 377 | } |
---|
| 378 | |
---|
| 379 | return retSwaps; |
---|
| 380 | } |
---|
| 381 | /** |
---|
| 382 | Computes a reduced row echelon form from a row echelon form. In |
---|
| 383 | Groebner basis terms this is the same as tranforming a minimal |
---|
| 384 | Groebner basis to a reduced one except that we do not force |
---|
| 385 | pivots to be 1 unless the scalePivotsToOne parameter is set. |
---|
| 386 | */ |
---|
| 387 | int REformToRREform(bool scalePivotsToOne=false) |
---|
| 388 | { |
---|
| 389 | int ret=0; |
---|
| 390 | int pivotI=-1; |
---|
| 391 | int pivotJ=-1; |
---|
| 392 | while(nextPivot(pivotI,pivotJ)) |
---|
| 393 | { |
---|
| 394 | if(scalePivotsToOne) |
---|
| 395 | rows[pivotI]=rows[pivotI]/rows[pivotI][pivotJ]; |
---|
| 396 | for(int i=0;i<pivotI;i++) |
---|
| 397 | if(!rows[i][pivotJ].isZero()) |
---|
| 398 | madd(pivotI,-rows[i][pivotJ]/rows[pivotI][pivotJ],i); |
---|
| 399 | } |
---|
| 400 | return ret; |
---|
| 401 | } |
---|
| 402 | /** |
---|
| 403 | This function may be called if the FieldMatrix is in Row Echelon |
---|
| 404 | Form. The input is a FieldVector which is rewritten modulo the |
---|
| 405 | rows of the matrix. The result is unique and is the same as the |
---|
| 406 | normal form of the input vector modulo the Groebner basis of the |
---|
| 407 | linear ideal generated by the rows of the matrix. |
---|
| 408 | */ |
---|
| 409 | Vector<typ> canonicalize(Vector<typ> v)const |
---|
| 410 | { |
---|
| 411 | assert(typ::isField()); |
---|
| 412 | assert(v.size()==getWidth()); |
---|
| 413 | |
---|
| 414 | int pivotI=-1; |
---|
| 415 | int pivotJ=-1; |
---|
| 416 | |
---|
| 417 | while(nextPivot(pivotI,pivotJ)) |
---|
| 418 | if(!v[pivotJ].isZero()) |
---|
| 419 | { |
---|
| 420 | typ s=-v[pivotJ]/rows[pivotI][pivotJ]; |
---|
| 421 | |
---|
| 422 | for(int k=0;k<width;k++) |
---|
| 423 | if(!rows[pivotI][k].isZero()) |
---|
| 424 | v[k].madd(rows[pivotI][k],s); |
---|
| 425 | } |
---|
| 426 | return v; |
---|
| 427 | } |
---|
| 428 | /** |
---|
| 429 | Calls reduce() and constructs matrix whose rows forms a basis of |
---|
| 430 | the kernel of the linear map defined by the original matrix. The |
---|
| 431 | return value is the new matrix. |
---|
| 432 | */ |
---|
| 433 | Matrix reduceAndComputeKernel() |
---|
| 434 | { |
---|
| 435 | Matrix ret(width-reduceAndComputeRank(),width); |
---|
| 436 | |
---|
| 437 | REformToRREform(); |
---|
| 438 | |
---|
| 439 | int k=0; |
---|
| 440 | int pivotI=-1; |
---|
| 441 | int pivotJ=-1; |
---|
| 442 | bool pivotExists=nextPivot(pivotI,pivotJ); |
---|
| 443 | for(int j=0;j<width;j++) |
---|
| 444 | { |
---|
| 445 | if(pivotExists && (pivotJ==j)) |
---|
| 446 | { |
---|
| 447 | pivotExists=nextPivot(pivotI,pivotJ); |
---|
| 448 | continue; |
---|
| 449 | } |
---|
| 450 | int pivot2I=-1; |
---|
| 451 | int pivot2J=-1; |
---|
| 452 | while(nextPivot(pivot2I,pivot2J)) |
---|
| 453 | { |
---|
| 454 | ret[k][pivot2J]=rows[pivot2I][j]/rows[pivot2I][pivot2J]; |
---|
| 455 | } |
---|
| 456 | ret[k][j]=typ(-1); |
---|
| 457 | k++; |
---|
| 458 | } |
---|
| 459 | return ret; |
---|
| 460 | } |
---|
| 461 | /** |
---|
| 462 | Assumes that the matrix has a kernel of dimension 1. |
---|
| 463 | Calls reduce() and returns a non-zero vector in the kernel. |
---|
| 464 | If the matrix is an (n-1)x(n) matrix then the returned vector has |
---|
| 465 | the property that if it was appended as a row to the original matrix |
---|
| 466 | then the determinant of that matrix would be positive. Of course |
---|
| 467 | this property, as described here, only makes sense for ordered fields. |
---|
| 468 | Only allowed for fields at the moment. |
---|
| 469 | */ |
---|
| 470 | Vector<typ> reduceAndComputeVectorInKernel() |
---|
| 471 | { |
---|
| 472 | assert(typ::isField()); |
---|
| 473 | // TODO: (optimization) if the field is ordered, then it is better to just keep track of signs rather than |
---|
| 474 | // multiplying by sign*diagonalProduct*lastEntry at the end. |
---|
| 475 | int nswaps=this->reduce(); |
---|
| 476 | typ sign=typ(1-2*(nswaps&1)); |
---|
| 477 | int rank=reduceAndComputeRank(); |
---|
| 478 | assert(rank+1==width); |
---|
| 479 | |
---|
| 480 | REformToRREform(); |
---|
| 481 | |
---|
| 482 | Vector<typ> ret(width); |
---|
| 483 | |
---|
| 484 | typ diagonalProduct(1); |
---|
| 485 | { |
---|
| 486 | int pivot2I=-1; |
---|
| 487 | int pivot2J=-1; |
---|
| 488 | while(nextPivot(pivot2I,pivot2J)) |
---|
| 489 | { |
---|
| 490 | diagonalProduct*=rows[pivot2I][pivot2J]; |
---|
| 491 | } |
---|
| 492 | } |
---|
| 493 | { |
---|
| 494 | int j=nonPivotColumns().front(); |
---|
| 495 | int pivot2I=-1; |
---|
| 496 | int pivot2J=-1; |
---|
| 497 | ret[j]=typ(-1); |
---|
| 498 | // Pretend that we are appending ret to the matrix, and reducing this |
---|
| 499 | // new row by the previous ones. The last entry of the resulting matrix |
---|
| 500 | // is lastEntry. |
---|
| 501 | typ lastEntry=ret[j]; |
---|
| 502 | while(nextPivot(pivot2I,pivot2J)) |
---|
| 503 | { |
---|
| 504 | ret[pivot2J]=rows[pivot2I][j]/rows[pivot2I][pivot2J]; |
---|
| 505 | lastEntry-=ret[pivot2J]*ret[pivot2J]; |
---|
| 506 | } |
---|
| 507 | ret=(sign*(diagonalProduct*lastEntry))*ret; |
---|
| 508 | } |
---|
| 509 | |
---|
| 510 | return ret; |
---|
| 511 | } |
---|
| 512 | |
---|
| 513 | /** |
---|
| 514 | Calls reduce() and returns the rank of the matrix. |
---|
| 515 | */ |
---|
| 516 | int reduceAndComputeRank() |
---|
| 517 | { |
---|
[26b713] | 518 | if (typ::isField()) |
---|
| 519 | reduce(); |
---|
| 520 | else |
---|
| 521 | reduce(false,true,false); |
---|
[def863] | 522 | int ret=0; |
---|
| 523 | int pivotI=-1; |
---|
| 524 | int pivotJ=-1; |
---|
| 525 | while(nextPivot(pivotI,pivotJ))ret++; |
---|
| 526 | return ret; |
---|
| 527 | } |
---|
| 528 | /** |
---|
| 529 | * Sort the rows of the matrix. |
---|
| 530 | */ |
---|
| 531 | void sortRows() |
---|
| 532 | { |
---|
| 533 | std::sort(rows.begin(),rows.end()); |
---|
| 534 | } |
---|
| 535 | /** |
---|
| 536 | * Sort the rows of the matrix and remove duplicate rows. |
---|
| 537 | */ |
---|
| 538 | void sortAndRemoveDuplicateRows() |
---|
| 539 | { |
---|
| 540 | sortRows(); |
---|
| 541 | if(getHeight()==0)return; |
---|
| 542 | Matrix B(0,getWidth()); |
---|
| 543 | B.appendRow((*this)[0]); |
---|
| 544 | for(int i=1;i<getHeight();i++) |
---|
| 545 | if(rows[i]!=rows[i-1])B.appendRow((*this)[i]); |
---|
| 546 | *this=B; |
---|
| 547 | } |
---|
| 548 | /** |
---|
| 549 | Takes two matrices with the same number of columns and construct |
---|
| 550 | a new matrix which has the rows of the matrix top on the top and |
---|
| 551 | the rows of the matrix bottom at the bottom. The return value is |
---|
| 552 | the constructed matrix. |
---|
| 553 | */ |
---|
| 554 | friend Matrix combineOnTop(Matrix const &top, Matrix const &bottom) |
---|
| 555 | { |
---|
| 556 | assert(bottom.getWidth()==top.getWidth()); |
---|
| 557 | Matrix ret(top.getHeight()+bottom.getHeight(),top.getWidth()); |
---|
| 558 | for(int i=0;i<top.getHeight();i++)ret.rows[i]=top.rows[i]; |
---|
| 559 | for(int i=0;i<bottom.getHeight();i++)ret.rows[i+top.getHeight()]=bottom.rows[i]; |
---|
| 560 | |
---|
| 561 | return ret; |
---|
| 562 | } |
---|
| 563 | /** |
---|
| 564 | Takes two matrices with the same number of rows and construct |
---|
| 565 | a new matrix which has the columns of the matrix left on the left and |
---|
| 566 | the columns of the matrix right on the right. The return value is |
---|
| 567 | the constructed matrix. |
---|
| 568 | */ |
---|
| 569 | friend Matrix combineLeftRight(Matrix const &left, Matrix const &right) |
---|
| 570 | { |
---|
| 571 | assert(left.getHeight()==right.getHeight()); |
---|
| 572 | Matrix ret(left.getHeight(),left.getWidth()+right.getWidth()); |
---|
| 573 | for(int i=0;i<left.getHeight();i++) |
---|
| 574 | { |
---|
| 575 | for(int j=0;j<left.getWidth();j++)ret.rows[i][j]=left.rows[i][j]; |
---|
| 576 | for(int j=0;j<right.getWidth();j++)ret.rows[i][j+left.getWidth()]=right.rows[i][j]; |
---|
| 577 | } |
---|
| 578 | return ret; |
---|
| 579 | } |
---|
| 580 | }; |
---|
| 581 | |
---|
| 582 | typedef Matrix<Integer> ZMatrix; |
---|
| 583 | typedef Matrix<Rational> QMatrix; |
---|
[c4d065] | 584 | typedef Matrix<int> IntMatrix; |
---|
[def863] | 585 | |
---|
| 586 | inline QMatrix ZToQMatrix(ZMatrix const &m) |
---|
| 587 | { |
---|
| 588 | QMatrix ret(m.getHeight(),m.getWidth()); |
---|
| 589 | for(int i=0;i<m.getHeight();i++)ret[i]=ZToQVector(m[i]); |
---|
| 590 | return ret; |
---|
| 591 | } |
---|
| 592 | |
---|
| 593 | inline ZMatrix QToZMatrixPrimitive(QMatrix const &m) |
---|
| 594 | { |
---|
| 595 | ZMatrix ret(m.getHeight(),m.getWidth()); |
---|
| 596 | for(int i=0;i<m.getHeight();i++)ret[i]=QToZVectorPrimitive(m[i]); |
---|
| 597 | return ret; |
---|
| 598 | } |
---|
[74a91c9] | 599 | |
---|
[c4d065] | 600 | |
---|
| 601 | inline IntMatrix ZToIntMatrix(ZMatrix const &m) |
---|
| 602 | { |
---|
| 603 | IntMatrix ret(m.getHeight(),m.getWidth()); |
---|
| 604 | for(int i=0;i<m.getHeight();i++)ret[i]=ZToIntVector(m[i]); |
---|
| 605 | return ret; |
---|
| 606 | } |
---|
| 607 | |
---|
| 608 | |
---|
| 609 | inline ZMatrix IntToZMatrix(IntMatrix const &m) |
---|
| 610 | { |
---|
| 611 | ZMatrix ret(m.getHeight(),m.getWidth()); |
---|
| 612 | for(int i=0;i<m.getHeight();i++)ret[i]=IntToZVector(m[i]); |
---|
| 613 | return ret; |
---|
| 614 | } |
---|
| 615 | |
---|
[74a91c9] | 616 | inline QMatrix canonicalizeSubspace(QMatrix const &m) |
---|
| 617 | { |
---|
| 618 | QMatrix temp=m; |
---|
| 619 | temp.reduce(); |
---|
| 620 | temp.REformToRREform(); |
---|
| 621 | temp.removeZeroRows(); |
---|
| 622 | return temp; |
---|
| 623 | } |
---|
| 624 | |
---|
| 625 | inline ZMatrix canonicalizeSubspace(ZMatrix const &m) |
---|
| 626 | { |
---|
| 627 | return QToZMatrixPrimitive(canonicalizeSubspace(ZToQMatrix(m))); |
---|
| 628 | } |
---|
| 629 | |
---|
| 630 | |
---|
| 631 | inline QMatrix kernel(QMatrix const &m) |
---|
| 632 | { |
---|
| 633 | QMatrix temp=m; |
---|
| 634 | return temp.reduceAndComputeKernel(); |
---|
| 635 | } |
---|
| 636 | |
---|
| 637 | inline ZMatrix kernel(ZMatrix const &m) |
---|
| 638 | { |
---|
| 639 | return QToZMatrixPrimitive(kernel(ZToQMatrix(m))); |
---|
| 640 | } |
---|
| 641 | |
---|
[def863] | 642 | } |
---|
| 643 | |
---|
| 644 | |
---|
| 645 | #endif /* LIB_ZMATRIX_H_ */ |
---|