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 | std::vector<typ> data; |
---|
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()),data(a.data){ |
---|
28 | } |
---|
29 | Matrix(int height_, int width_):width(width_),height(height_),data(width_*height_){ |
---|
30 | assert(height>=0); |
---|
31 | assert(width>=0); |
---|
32 | }; |
---|
33 | ~Matrix(){ |
---|
34 | }; |
---|
35 | Matrix():width(0),height(0){ |
---|
36 | }; |
---|
37 | static Matrix rowVectorMatrix(Vector<typ> const &v) |
---|
38 | { |
---|
39 | Matrix ret(1,v.size()); |
---|
40 | for(int i=0;i<v.size();i++)ret[0][i]=v[i]; |
---|
41 | return ret; |
---|
42 | } |
---|
43 | Vector<typ> column(int i)const |
---|
44 | { |
---|
45 | assert(i>=0); |
---|
46 | assert(i<getWidth()); |
---|
47 | Vector<typ> ret(getHeight()); |
---|
48 | for(int j=0;j<getHeight();j++)ret[j]=(*this)[j][i]; |
---|
49 | return ret; |
---|
50 | } |
---|
51 | Matrix transposed()const |
---|
52 | { |
---|
53 | Matrix ret(getWidth(),getHeight()); |
---|
54 | for(int i=0;i<getWidth();i++) |
---|
55 | for(int j=0;j<getHeight();j++) |
---|
56 | ret[i][j]=(*this)[j][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[i][i]=typ(1); |
---|
63 | return m; |
---|
64 | } |
---|
65 | void append(Matrix const &m) |
---|
66 | { |
---|
67 | assert(m.getWidth()==width); |
---|
68 | data.resize((height+m.height)*width); |
---|
69 | int oldHeight=height; |
---|
70 | height+=m.height; |
---|
71 | for(int i=0;i<m.height;i++) |
---|
72 | { |
---|
73 | for(int j=0;j<m.width;j++) |
---|
74 | (*this)[i+oldHeight][j]=m[i][j]; |
---|
75 | } |
---|
76 | } |
---|
77 | void appendRow(Vector<typ> const &v) |
---|
78 | { |
---|
79 | assert(v.size()==width); |
---|
80 | data.resize((height+1)*width); |
---|
81 | height++; |
---|
82 | for(int j=0;j<width;j++) |
---|
83 | (*this)[height-1][j]=v[j]; |
---|
84 | } |
---|
85 | void eraseLastRow() |
---|
86 | { |
---|
87 | assert(height>0); |
---|
88 | data.resize((height-1)*width); |
---|
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++) |
---|
114 | for(int j=0;j<q.width;j++) |
---|
115 | p[i][j]=s*(q[i][j]); |
---|
116 | return p; |
---|
117 | } |
---|
118 | friend Matrix operator*(const Matrix& a, const Matrix& b) |
---|
119 | { |
---|
120 | assert(a.width==b.height); |
---|
121 | Matrix ret(b.width,a.height); |
---|
122 | for(int i=0;i<b.width;i++) |
---|
123 | ret[i]=a.vectormultiply(b.column(i)); |
---|
124 | return ret.transposed(); |
---|
125 | } |
---|
126 | /* template<class T> |
---|
127 | Matrix<T>(const Matrix<T>& c):v(c.size()){ |
---|
128 | for(int i=0;i<size();i++)v[i]=typ(c[i]);} |
---|
129 | */ |
---|
130 | friend Matrix operator-(const Matrix &b) |
---|
131 | { |
---|
132 | Matrix ret(b.height,b.width); |
---|
133 | for(int i=0;i<b.height;i++)ret[i]=-b[i]; |
---|
134 | return ret; |
---|
135 | } |
---|
136 | |
---|
137 | /** |
---|
138 | Returns the specified submatrix. The endRow and endColumn are not included. |
---|
139 | */ |
---|
140 | Matrix submatrix(int startRow, int startColumn, int endRow, int endColumn)const |
---|
141 | { |
---|
142 | assert(startRow>=0); |
---|
143 | assert(startColumn>=0); |
---|
144 | assert(endRow>=startRow); |
---|
145 | assert(endColumn>=startColumn); |
---|
146 | assert(endRow<=height); |
---|
147 | assert(endColumn<=width); |
---|
148 | Matrix ret(endRow-startRow,endColumn-startColumn); |
---|
149 | for(int i=startRow;i<endRow;i++) |
---|
150 | for(int j=startColumn;j<endColumn;j++) |
---|
151 | ret[i-startRow][j-startColumn]=(*this)[i][j]; |
---|
152 | return ret; |
---|
153 | } |
---|
154 | |
---|
155 | class RowRef; |
---|
156 | class const_RowRef{ |
---|
157 | int rowNumM; |
---|
158 | Matrix const &matrix; |
---|
159 | friend class RowRef; |
---|
160 | public: |
---|
161 | inline const_RowRef(const Matrix &matrix_, int rowNum_)__attribute__((always_inline)): |
---|
162 | rowNumM(rowNum_*matrix_.width), |
---|
163 | matrix(matrix_) |
---|
164 | { |
---|
165 | } |
---|
166 | inline typ const &operator[](int j)const __attribute__((always_inline)) |
---|
167 | { |
---|
168 | assert(j>=0); |
---|
169 | assert(j<matrix.width); |
---|
170 | return matrix.data[rowNumM+j]; |
---|
171 | } |
---|
172 | inline typ const &UNCHECKEDACCESS(int j)const __attribute__((always_inline)) |
---|
173 | { |
---|
174 | return matrix.data[rowNumM+j]; |
---|
175 | } |
---|
176 | const Vector<typ> toVector()const |
---|
177 | { |
---|
178 | Vector<typ> ret(matrix.width); |
---|
179 | for(int j=0;j<matrix.width;j++) |
---|
180 | ret[j]=matrix.data[rowNumM+j]; |
---|
181 | return ret; |
---|
182 | } |
---|
183 | operator Vector<typ>()const |
---|
184 | { |
---|
185 | return toVector(); |
---|
186 | } |
---|
187 | bool operator==(Vector<typ> const &b)const |
---|
188 | { |
---|
189 | return toVector()==b; |
---|
190 | } |
---|
191 | /* typ dot(Vector<typ> const &b)const |
---|
192 | { |
---|
193 | return dot(toVector(),b); |
---|
194 | }*/ |
---|
195 | Vector<typ> operator-()const |
---|
196 | { |
---|
197 | return -toVector(); |
---|
198 | } |
---|
199 | }; |
---|
200 | class RowRef{ |
---|
201 | int rowNumM; |
---|
202 | Matrix &matrix; |
---|
203 | public: |
---|
204 | inline RowRef(Matrix &matrix_, int rowNum_): |
---|
205 | rowNumM(rowNum_*matrix_.width), |
---|
206 | matrix(matrix_) |
---|
207 | { |
---|
208 | } |
---|
209 | inline typ &operator[](int j) __attribute__((always_inline)) |
---|
210 | { |
---|
211 | assert(j>=0); |
---|
212 | assert(j<matrix.width); |
---|
213 | return matrix.data[rowNumM+j]; |
---|
214 | } |
---|
215 | inline typ &UNCHECKEDACCESS(int j) |
---|
216 | { |
---|
217 | return matrix.data[rowNumM+j]; |
---|
218 | } |
---|
219 | RowRef &operator=(Vector<typ> const &v) |
---|
220 | { |
---|
221 | assert(v.size()==matrix.width); |
---|
222 | for(int j=0;j<matrix.width;j++) |
---|
223 | matrix.data[rowNumM+j]=v[j]; |
---|
224 | |
---|
225 | return *this; |
---|
226 | } |
---|
227 | RowRef &operator=(RowRef const &v) |
---|
228 | { |
---|
229 | assert(v.matrix.width==matrix.width); |
---|
230 | for(int j=0;j<matrix.width;j++) |
---|
231 | matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j]; |
---|
232 | |
---|
233 | return *this; |
---|
234 | } |
---|
235 | /* RowRef &operator+=(Vector<typ> const &v) |
---|
236 | { |
---|
237 | assert(v.size()==matrix.width); |
---|
238 | for(int j=0;j<matrix.width;j++) |
---|
239 | matrix.data[rowNumM+j]+=v.v[j]; |
---|
240 | |
---|
241 | return *this; |
---|
242 | }*/ |
---|
243 | RowRef &operator+=(RowRef const &v) |
---|
244 | { |
---|
245 | assert(v.matrix.width==matrix.width); |
---|
246 | for(int j=0;j<matrix.width;j++) |
---|
247 | matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j]; |
---|
248 | |
---|
249 | return *this; |
---|
250 | } |
---|
251 | RowRef &operator+=(const_RowRef const &v) |
---|
252 | { |
---|
253 | assert(v.matrix.width==matrix.width); |
---|
254 | for(int j=0;j<matrix.width;j++) |
---|
255 | matrix.data[rowNumM+j]+=v.matrix.data[v.rowNumM+j]; |
---|
256 | |
---|
257 | return *this; |
---|
258 | } |
---|
259 | RowRef &operator=(const_RowRef const &v) |
---|
260 | { |
---|
261 | assert(v.matrix.width==matrix.width); |
---|
262 | for(int j=0;j<matrix.width;j++) |
---|
263 | matrix.data[rowNumM+j]=v.matrix.data[v.rowNumM+j]; |
---|
264 | |
---|
265 | return *this; |
---|
266 | } |
---|
267 | const Vector<typ> toVector()const |
---|
268 | { |
---|
269 | Vector<typ> ret(matrix.width); |
---|
270 | for(int j=0;j<matrix.width;j++) |
---|
271 | ret[j]=matrix.data[rowNumM+j]; |
---|
272 | return ret; |
---|
273 | } |
---|
274 | operator Vector<typ>()const |
---|
275 | { |
---|
276 | return toVector(); |
---|
277 | } |
---|
278 | /* typ dot(Vector<typ> const &b)const |
---|
279 | { |
---|
280 | return dot(toVector(),b); |
---|
281 | }*/ |
---|
282 | bool isZero()const |
---|
283 | { |
---|
284 | for(int j=0;j<matrix.width;j++)if(!(matrix.data[rowNumM+j].isZero()))return false; |
---|
285 | return true; |
---|
286 | } |
---|
287 | }; |
---|
288 | |
---|
289 | |
---|
290 | inline RowRef operator[](int i) __attribute__((always_inline)) |
---|
291 | { |
---|
292 | assert(i>=0); |
---|
293 | assert(i<height); |
---|
294 | return RowRef(*this,i); |
---|
295 | } |
---|
296 | inline const_RowRef operator[](int i)const __attribute__((always_inline)) |
---|
297 | { |
---|
298 | assert(i>=0); |
---|
299 | assert(i<height); |
---|
300 | return const_RowRef(*this,i); |
---|
301 | } |
---|
302 | |
---|
303 | |
---|
304 | const typ& UNCHECKEDACCESS(int i,int j)const __attribute__((always_inline)){ |
---|
305 | /* assert(i>=0); |
---|
306 | assert(i<height); |
---|
307 | assert(j>=0); |
---|
308 | assert(j<width);*/ |
---|
309 | return data[i*width+j];} |
---|
310 | typ& UNCHECKEDACCESS(int i,int j) __attribute__((always_inline)){ |
---|
311 | /* assert(i>=0); |
---|
312 | assert(i<height); |
---|
313 | assert(j>=0); |
---|
314 | assert(j<width);*/ |
---|
315 | return data[i*width+j];} |
---|
316 | |
---|
317 | |
---|
318 | |
---|
319 | bool operator<(const Matrix & b)const |
---|
320 | { |
---|
321 | if(getWidth()<b.getWidth())return true; |
---|
322 | if(b.getWidth()<getWidth())return false; |
---|
323 | if(getHeight()<b.getHeight())return true; |
---|
324 | if(b.getHeight()<getHeight())return false; |
---|
325 | |
---|
326 | for(int i=0;i<getHeight();i++) |
---|
327 | { |
---|
328 | if((*this)[i].toVector()<b[i].toVector())return true; |
---|
329 | if(b[i].toVector()<(*this)[i].toVector())return false; |
---|
330 | } |
---|
331 | return false; |
---|
332 | } |
---|
333 | /** |
---|
334 | Adds a times the i th row to the j th row. |
---|
335 | */ |
---|
336 | void madd(int i, typ a, int j) |
---|
337 | { |
---|
338 | assert(i!=j); |
---|
339 | assert(i>=0 && i<height); |
---|
340 | assert(j>=0 && j<height); |
---|
341 | |
---|
342 | if(!a.isZero()) |
---|
343 | for(int k=0;k<width;k++) |
---|
344 | if(!(*this)[i][k].isZero()) |
---|
345 | (*this)[j][k].madd((*this)[i][k],a); |
---|
346 | } |
---|
347 | |
---|
348 | friend std::ostream &operator<<(std::ostream &f, Matrix const &a){ |
---|
349 | f<<"{"; |
---|
350 | for(int i=0;i<a.getHeight();i++) |
---|
351 | { |
---|
352 | if(i)f<<","<<std::endl; |
---|
353 | f<<a[i].toVector(); |
---|
354 | } |
---|
355 | f<<"}"<<std::endl; |
---|
356 | return f; |
---|
357 | } |
---|
358 | |
---|
359 | std::string toString()const |
---|
360 | { |
---|
361 | std::stringstream f; |
---|
362 | f<<*this; |
---|
363 | return f.str(); |
---|
364 | } |
---|
365 | |
---|
366 | /** |
---|
367 | Swaps the i th and the j th row. |
---|
368 | */ |
---|
369 | void swapRows(int i, int j) |
---|
370 | { |
---|
371 | for(int a=0;a<width;a++)std::swap((*this)[i][a],(*this)[j][a]); |
---|
372 | } |
---|
373 | /** |
---|
374 | This method is used for iterating through the pivots in a matrix |
---|
375 | in row echelon form. To find the first pivot put i=-1 and |
---|
376 | j=-1 and call this routine. The (i,j) th entry of the matrix is a |
---|
377 | pivot. Call the routine again to get the next pivot. When no more |
---|
378 | pivots are found the routine returns false. |
---|
379 | */ |
---|
380 | bool nextPivot(int &i, int &j)const |
---|
381 | { |
---|
382 | i++; |
---|
383 | if(i>=height)return false; |
---|
384 | while(++j<width) |
---|
385 | { |
---|
386 | if(!(*this)[i][j].isZero()) return true; |
---|
387 | } |
---|
388 | return false; |
---|
389 | } |
---|
390 | /** |
---|
391 | Returns the indices of the columns containing a pivot. |
---|
392 | The returned list is sorted. |
---|
393 | The matrix must be in row echelon form. |
---|
394 | */ |
---|
395 | std::vector<int> pivotColumns()const |
---|
396 | { |
---|
397 | std::vector<int> ret; |
---|
398 | int pivotI=-1; |
---|
399 | int pivotJ=-1; |
---|
400 | while(nextPivot(pivotI,pivotJ))ret.push_back(pivotJ); |
---|
401 | return ret; |
---|
402 | } |
---|
403 | /** |
---|
404 | Returns the indices of the columns not containing a pivot. |
---|
405 | The returned list is sorted. |
---|
406 | The matrix must be in row echelon form. |
---|
407 | */ |
---|
408 | std::vector<int> nonPivotColumns()const |
---|
409 | { |
---|
410 | std::vector<int> ret; |
---|
411 | int pivotI=-1; |
---|
412 | int pivotJ=-1; |
---|
413 | int firstPossiblePivot=0; |
---|
414 | while(nextPivot(pivotI,pivotJ)) |
---|
415 | { |
---|
416 | for(int j=firstPossiblePivot;j<pivotJ;j++) |
---|
417 | ret.push_back(j); |
---|
418 | firstPossiblePivot=pivotJ+1; |
---|
419 | } |
---|
420 | for(int j=firstPossiblePivot;j<getWidth();j++) |
---|
421 | ret.push_back(j); |
---|
422 | |
---|
423 | return ret; |
---|
424 | } |
---|
425 | /** |
---|
426 | This routine removes the zero rows of the matrix. |
---|
427 | */ |
---|
428 | void removeZeroRows() |
---|
429 | { |
---|
430 | int nonZeros=0; |
---|
431 | for(int i=0;i<height;i++)if(!(*this)[i].isZero())nonZeros++; |
---|
432 | if(nonZeros==height)return; |
---|
433 | |
---|
434 | Matrix b(nonZeros,width); |
---|
435 | |
---|
436 | int j=0; |
---|
437 | for(int i=0;i<height;i++) |
---|
438 | { |
---|
439 | if(!(*this)[i].isZero()) |
---|
440 | { |
---|
441 | b[j]=(*this)[i]; |
---|
442 | j++; |
---|
443 | } |
---|
444 | } |
---|
445 | *this=b; |
---|
446 | } |
---|
447 | /** |
---|
448 | Returns the index of a row whose index is at least currentRow and |
---|
449 | which has a non-zero element on the column th entry. If no such |
---|
450 | row exists then -1 is returned. This routine is used in the Gauss |
---|
451 | reduction. To make the reduction more efficient the routine |
---|
452 | chooses its row with as few non-zero entries as possible. |
---|
453 | */ |
---|
454 | int findRowIndex(int column, int currentRow)const |
---|
455 | { |
---|
456 | int best=-1; |
---|
457 | int bestNumberOfNonZero=0; |
---|
458 | for(int i=currentRow;i<height;i++) |
---|
459 | if(!(*this)[i][column].isZero()) |
---|
460 | { |
---|
461 | int nz=0; |
---|
462 | for(int k=column+1;k<width;k++) |
---|
463 | if(!(*this)[i][k].isZero())nz++; |
---|
464 | if(best==-1) |
---|
465 | { |
---|
466 | best=i; |
---|
467 | bestNumberOfNonZero=nz; |
---|
468 | } |
---|
469 | else if(nz<bestNumberOfNonZero) |
---|
470 | { |
---|
471 | best=i; |
---|
472 | bestNumberOfNonZero=nz; |
---|
473 | } |
---|
474 | } |
---|
475 | return best; |
---|
476 | } |
---|
477 | /** |
---|
478 | Performs a Gauss reduction and returns the number of row swaps (and negative scalings) |
---|
479 | done. The result is a matrix in row echelon form. The pivots may |
---|
480 | not be all 1. In terms of Groebner bases, what is computed is a |
---|
481 | minimal (not necessarily reduced) Groebner basis of the linear |
---|
482 | ideal generated by the rows. The number of swaps is need if one |
---|
483 | wants to compute the determinant afterwards. In this case it is |
---|
484 | also a good idea to set the flag returnIfZeroDeterminant which |
---|
485 | make the routine terminate before completion if it discovers that |
---|
486 | the determinant is zero. |
---|
487 | */ |
---|
488 | int reduce(bool returnIfZeroDeterminant=false, bool integral=false, bool makePivotsOne=false) |
---|
489 | { |
---|
490 | assert(integral || typ::isField()); |
---|
491 | assert(!makePivotsOne || !integral); |
---|
492 | |
---|
493 | int retSwaps=0; |
---|
494 | int currentRow=0; |
---|
495 | |
---|
496 | for(int i=0;i<width;i++) |
---|
497 | { |
---|
498 | int s=findRowIndex(i,currentRow); |
---|
499 | |
---|
500 | if(s!=-1) |
---|
501 | { |
---|
502 | if(s!=currentRow) |
---|
503 | { |
---|
504 | swapRows(currentRow,s); |
---|
505 | retSwaps++; |
---|
506 | } |
---|
507 | if(makePivotsOne) |
---|
508 | {//THE PIVOT SHOULD BE SET TO ONE IF INTEGRAL IS FALSE |
---|
509 | if((*this)[currentRow][i].sign()>=0)retSwaps++; |
---|
510 | typ inverse=typ(1)/(*this)[currentRow][i]; |
---|
511 | // if(!rows[currentRow][i].isOne()) |
---|
512 | { |
---|
513 | for(int k=0;k<width;k++) |
---|
514 | if(!(*this)[currentRow][k].isZero()) |
---|
515 | (*this)[currentRow][k]*=inverse; |
---|
516 | } |
---|
517 | } |
---|
518 | for(int j=currentRow+1;j<height;j++) |
---|
519 | if(integral) |
---|
520 | { |
---|
521 | if(!(*this)[j][i].isZero()) |
---|
522 | { |
---|
523 | typ s; |
---|
524 | typ t; |
---|
525 | |
---|
526 | typ g=typ::gcd((*this)[currentRow][i],(*this)[j][i],s,t); |
---|
527 | typ u=-(*this)[j][i]/g; |
---|
528 | typ v=(*this)[currentRow][i]/g; |
---|
529 | /* We want the (s,t) vector to be as small as possible. |
---|
530 | * We are allowed to adjust by multiples of (u,v). |
---|
531 | * The following computes the correct multiplier (in most cases). |
---|
532 | */ |
---|
533 | /* { |
---|
534 | FieldElement multiplier=(s*u+t*v)*((u*u+v*v).inverse()); |
---|
535 | double d=mpq_get_d(*(multiplier.getGmpRationalTemporaryPointer())); |
---|
536 | multiplier=multiplier.getField()->zHomomorphism(-(((int)(d+10000.5))-10000)); |
---|
537 | s.madd(multiplier,u); |
---|
538 | t.madd(multiplier,v); |
---|
539 | }*/ |
---|
540 | for(int k=0;k<width;k++) |
---|
541 | { |
---|
542 | typ A=(*this)[currentRow][k]; |
---|
543 | typ B=(*this)[j][k]; |
---|
544 | |
---|
545 | (*this)[currentRow][k]=s*A+t*B; |
---|
546 | (*this)[j][k]=u*A+v*B; |
---|
547 | } |
---|
548 | } |
---|
549 | } |
---|
550 | else |
---|
551 | { |
---|
552 | if(!(*this)[j][i].isZero()) |
---|
553 | madd(currentRow,-(*this)[j][i]/(*this)[currentRow][i],j); |
---|
554 | } |
---|
555 | currentRow++; |
---|
556 | } |
---|
557 | else |
---|
558 | if(returnIfZeroDeterminant)return -1; |
---|
559 | } |
---|
560 | |
---|
561 | return retSwaps; |
---|
562 | } |
---|
563 | /** |
---|
564 | Computes a reduced row echelon form from a row echelon form. In |
---|
565 | Groebner basis terms this is the same as tranforming a minimal |
---|
566 | Groebner basis to a reduced one except that we do not force |
---|
567 | pivots to be 1 unless the scalePivotsToOne parameter is set. |
---|
568 | */ |
---|
569 | int REformToRREform(bool scalePivotsToOne=false) |
---|
570 | { |
---|
571 | int ret=0; |
---|
572 | int pivotI=-1; |
---|
573 | int pivotJ=-1; |
---|
574 | while(nextPivot(pivotI,pivotJ)) |
---|
575 | { |
---|
576 | if(scalePivotsToOne) |
---|
577 | (*this)[pivotI]=(*this)[pivotI].toVector()/(*this)[pivotI][pivotJ]; |
---|
578 | for(int i=0;i<pivotI;i++) |
---|
579 | if(!(*this)[i][pivotJ].isZero()) |
---|
580 | madd(pivotI,-(*this)[i][pivotJ]/(*this)[pivotI][pivotJ],i); |
---|
581 | } |
---|
582 | return ret; |
---|
583 | } |
---|
584 | /** |
---|
585 | This function may be called if the FieldMatrix is in Row Echelon |
---|
586 | Form. The input is a FieldVector which is rewritten modulo the |
---|
587 | rows of the matrix. The result is unique and is the same as the |
---|
588 | normal form of the input vector modulo the Groebner basis of the |
---|
589 | linear ideal generated by the rows of the matrix. |
---|
590 | */ |
---|
591 | Vector<typ> canonicalize(Vector<typ> v)const |
---|
592 | { |
---|
593 | assert(typ::isField()); |
---|
594 | assert((int)v.size()==getWidth()); |
---|
595 | |
---|
596 | int pivotI=-1; |
---|
597 | int pivotJ=-1; |
---|
598 | |
---|
599 | while(nextPivot(pivotI,pivotJ)) |
---|
600 | if(!v[pivotJ].isZero()) |
---|
601 | { |
---|
602 | typ s=-v[pivotJ]/(*this)[pivotI][pivotJ]; |
---|
603 | |
---|
604 | for(int k=0;k<width;k++) |
---|
605 | if(!(*this)[pivotI][k].isZero()) |
---|
606 | v[k].madd((*this)[pivotI][k],s); |
---|
607 | } |
---|
608 | return v; |
---|
609 | } |
---|
610 | /** |
---|
611 | Calls reduce() and constructs matrix whose rows forms a basis of |
---|
612 | the kernel of the linear map defined by the original matrix. The |
---|
613 | return value is the new matrix. |
---|
614 | */ |
---|
615 | Matrix reduceAndComputeKernel() |
---|
616 | { |
---|
617 | Matrix ret(width-reduceAndComputeRank(),width); |
---|
618 | |
---|
619 | REformToRREform(); |
---|
620 | |
---|
621 | int k=0; |
---|
622 | int pivotI=-1; |
---|
623 | int pivotJ=-1; |
---|
624 | bool pivotExists=nextPivot(pivotI,pivotJ); |
---|
625 | for(int j=0;j<width;j++) |
---|
626 | { |
---|
627 | if(pivotExists && (pivotJ==j)) |
---|
628 | { |
---|
629 | pivotExists=nextPivot(pivotI,pivotJ); |
---|
630 | continue; |
---|
631 | } |
---|
632 | int pivot2I=-1; |
---|
633 | int pivot2J=-1; |
---|
634 | while(nextPivot(pivot2I,pivot2J)) |
---|
635 | { |
---|
636 | ret[k][pivot2J]=(*this)[pivot2I][j]/(*this)[pivot2I][pivot2J]; |
---|
637 | } |
---|
638 | ret[k][j]=typ(-1); |
---|
639 | k++; |
---|
640 | } |
---|
641 | return ret; |
---|
642 | } |
---|
643 | /** |
---|
644 | Assumes that the matrix has a kernel of dimension 1. |
---|
645 | Calls reduce() and returns a non-zero vector in the kernel. |
---|
646 | If the matrix is an (n-1)x(n) matrix then the returned vector has |
---|
647 | the property that if it was appended as a row to the original matrix |
---|
648 | then the determinant of that matrix would be positive. Of course |
---|
649 | this property, as described here, only makes sense for ordered fields. |
---|
650 | Only allowed for fields at the moment. |
---|
651 | */ |
---|
652 | Vector<typ> reduceAndComputeVectorInKernel() |
---|
653 | { |
---|
654 | assert(typ::isField()); |
---|
655 | // TODO: (optimization) if the field is ordered, then it is better to just keep track of signs rather than |
---|
656 | // multiplying by sign*diagonalProduct*lastEntry at the end. |
---|
657 | int nswaps=this->reduce(); |
---|
658 | typ sign=typ(1-2*(nswaps&1)); |
---|
659 | int rank=reduceAndComputeRank(); |
---|
660 | assert(rank+1==width); |
---|
661 | |
---|
662 | REformToRREform(); |
---|
663 | |
---|
664 | Vector<typ> ret(width); |
---|
665 | |
---|
666 | typ diagonalProduct(1); |
---|
667 | { |
---|
668 | int pivot2I=-1; |
---|
669 | int pivot2J=-1; |
---|
670 | while(nextPivot(pivot2I,pivot2J)) |
---|
671 | { |
---|
672 | diagonalProduct*=(*this)[pivot2I][pivot2J]; |
---|
673 | } |
---|
674 | } |
---|
675 | { |
---|
676 | int j=nonPivotColumns().front(); |
---|
677 | int pivot2I=-1; |
---|
678 | int pivot2J=-1; |
---|
679 | ret[j]=typ(-1); |
---|
680 | // Pretend that we are appending ret to the matrix, and reducing this |
---|
681 | // new row by the previous ones. The last entry of the resulting matrix |
---|
682 | // is lastEntry. |
---|
683 | typ lastEntry=ret[j]; |
---|
684 | while(nextPivot(pivot2I,pivot2J)) |
---|
685 | { |
---|
686 | ret[pivot2J]=(*this)[pivot2I][j]/(*this)[pivot2I][pivot2J]; |
---|
687 | lastEntry-=ret[pivot2J]*ret[pivot2J]; |
---|
688 | } |
---|
689 | ret=(sign*(diagonalProduct*lastEntry))*ret; |
---|
690 | } |
---|
691 | |
---|
692 | return ret; |
---|
693 | } |
---|
694 | |
---|
695 | /** |
---|
696 | Calls reduce() and returns the rank of the matrix. |
---|
697 | */ |
---|
698 | int reduceAndComputeRank() |
---|
699 | { |
---|
700 | reduce(false,!typ::isField(),false); |
---|
701 | int ret=0; |
---|
702 | int pivotI=-1; |
---|
703 | int pivotJ=-1; |
---|
704 | while(nextPivot(pivotI,pivotJ))ret++; |
---|
705 | return ret; |
---|
706 | } |
---|
707 | /** |
---|
708 | * Sort the rows of the matrix. |
---|
709 | */ |
---|
710 | struct rowComparer{ |
---|
711 | bool operator()(std::pair<Matrix*,int> i, std::pair<Matrix*,int> j) {return ((*i.first)[i.second].toVector()<(*j.first)[j.second].toVector());} |
---|
712 | } theRowComparer; |
---|
713 | void sortRows() |
---|
714 | { |
---|
715 | std::vector<std::pair<Matrix*,int> > v; |
---|
716 | for(int i=0;i<height;i++)v.push_back(std::pair<Matrix*,int>(this,i)); |
---|
717 | std::sort(v.begin(),v.end(),theRowComparer); |
---|
718 | Matrix result(height,width); |
---|
719 | for(int i=0;i<height;i++) |
---|
720 | result[i]=(*this)[v[i].second].toVector(); |
---|
721 | data=result.data; |
---|
722 | } |
---|
723 | /** |
---|
724 | * Sort the rows of the matrix and remove duplicate rows. |
---|
725 | */ |
---|
726 | void sortAndRemoveDuplicateRows() |
---|
727 | { |
---|
728 | sortRows(); |
---|
729 | if(getHeight()==0)return; |
---|
730 | Matrix B(0,getWidth()); |
---|
731 | B.appendRow((*this)[0]); |
---|
732 | for(int i=1;i<getHeight();i++) |
---|
733 | if((*this)[i].toVector()!=(*this)[i-1].toVector())B.appendRow((*this)[i].toVector()); |
---|
734 | *this=B; |
---|
735 | } |
---|
736 | /** |
---|
737 | Takes two matrices with the same number of columns and construct |
---|
738 | a new matrix which has the rows of the matrix top on the top and |
---|
739 | the rows of the matrix bottom at the bottom. The return value is |
---|
740 | the constructed matrix. |
---|
741 | */ |
---|
742 | friend Matrix combineOnTop(Matrix const &top, Matrix const &bottom) |
---|
743 | { |
---|
744 | assert(bottom.getWidth()==top.getWidth()); |
---|
745 | Matrix ret(top.getHeight()+bottom.getHeight(),top.getWidth()); |
---|
746 | for(int i=0;i<top.getHeight();i++)ret[i]=top[i]; |
---|
747 | for(int i=0;i<bottom.getHeight();i++)ret[i+top.getHeight()]=bottom[i]; |
---|
748 | |
---|
749 | return ret; |
---|
750 | } |
---|
751 | /** |
---|
752 | Takes two matrices with the same number of rows and construct |
---|
753 | a new matrix which has the columns of the matrix left on the left and |
---|
754 | the columns of the matrix right on the right. The return value is |
---|
755 | the constructed matrix. |
---|
756 | */ |
---|
757 | friend Matrix combineLeftRight(Matrix const &left, Matrix const &right) |
---|
758 | { |
---|
759 | assert(left.getHeight()==right.getHeight()); |
---|
760 | Matrix ret(left.getHeight(),left.getWidth()+right.getWidth()); |
---|
761 | for(int i=0;i<left.getHeight();i++) |
---|
762 | { |
---|
763 | for(int j=0;j<left.getWidth();j++)ret[i][j]=left[i][j]; |
---|
764 | for(int j=0;j<right.getWidth();j++)ret[i][j+left.getWidth()]=right[i][j]; |
---|
765 | } |
---|
766 | return ret; |
---|
767 | } |
---|
768 | }; |
---|
769 | |
---|
770 | typedef Matrix<Integer> ZMatrix; |
---|
771 | typedef Matrix<Rational> QMatrix; |
---|
772 | typedef Matrix<int> IntMatrix; |
---|
773 | |
---|
774 | inline QMatrix ZToQMatrix(ZMatrix const &m) |
---|
775 | { |
---|
776 | QMatrix ret(m.getHeight(),m.getWidth()); |
---|
777 | for(int i=0;i<m.getHeight();i++)ret[i]=ZToQVector(m[i]); |
---|
778 | return ret; |
---|
779 | } |
---|
780 | |
---|
781 | inline ZMatrix QToZMatrixPrimitive(QMatrix const &m) |
---|
782 | { |
---|
783 | ZMatrix ret(m.getHeight(),m.getWidth()); |
---|
784 | for(int i=0;i<m.getHeight();i++)ret[i]=QToZVectorPrimitive(m[i]); |
---|
785 | return ret; |
---|
786 | } |
---|
787 | |
---|
788 | |
---|
789 | inline IntMatrix ZToIntMatrix(ZMatrix const &m) |
---|
790 | { |
---|
791 | IntMatrix ret(m.getHeight(),m.getWidth()); |
---|
792 | for(int i=0;i<m.getHeight();i++)ret[i]=ZToIntVector(m[i]); |
---|
793 | return ret; |
---|
794 | } |
---|
795 | |
---|
796 | |
---|
797 | inline ZMatrix IntToZMatrix(IntMatrix const &m) |
---|
798 | { |
---|
799 | ZMatrix ret(m.getHeight(),m.getWidth()); |
---|
800 | for(int i=0;i<m.getHeight();i++)ret[i]=IntToZVector(m[i]); |
---|
801 | return ret; |
---|
802 | } |
---|
803 | |
---|
804 | inline QMatrix canonicalizeSubspace(QMatrix const &m) |
---|
805 | { |
---|
806 | QMatrix temp=m; |
---|
807 | temp.reduce(); |
---|
808 | temp.REformToRREform(); |
---|
809 | temp.removeZeroRows(); |
---|
810 | return temp; |
---|
811 | } |
---|
812 | |
---|
813 | inline ZMatrix canonicalizeSubspace(ZMatrix const &m) |
---|
814 | { |
---|
815 | return QToZMatrixPrimitive(canonicalizeSubspace(ZToQMatrix(m))); |
---|
816 | } |
---|
817 | |
---|
818 | |
---|
819 | inline QMatrix kernel(QMatrix const &m) |
---|
820 | { |
---|
821 | QMatrix temp=m; |
---|
822 | return temp.reduceAndComputeKernel(); |
---|
823 | } |
---|
824 | |
---|
825 | inline ZMatrix kernel(ZMatrix const &m) |
---|
826 | { |
---|
827 | return QToZMatrixPrimitive(kernel(ZToQMatrix(m))); |
---|
828 | } |
---|
829 | |
---|
830 | } |
---|
831 | |
---|
832 | |
---|
833 | #endif /* LIB_ZMATRIX_H_ */ |
---|