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