| 32 | if (p) sq = sq * sq; |
| 33 | } return res; } |
| 34 | matrix<T> rref(T &det, int &rank) { |
| 35 | matrix<T> mat(*this); det = T(1), rank = 0; |
| 36 | for (int r = 0, c = 0; c < cols; c++) { |
| 37 | int k = r; |
| 38 | rep(i,k+1,rows) if (abs(mat(i,c)) > abs(mat(k,c))) k = i; |
| 39 | if (k >= rows || eq<T>(mat(k, c), T(0))) continue; |
| 40 | if (k != r) { |
| 41 | det *= T(-1); |
| 42 | rep(i,0,cols) swap(mat.at(k, i), mat.at(r, i)); |
| 43 | } det *= mat(r, r); rank++; |
| 44 | T d = mat(r,c); |
| 45 | rep(i,0,cols) mat(r, i) /= d; |
| 46 | rep(i,0,rows) { |
| 47 | T m = mat(i, c); |
| 48 | if (i != r && !eq<T>(m, T(0))) |
| 49 | rep(j,0,cols) mat(i, j) -= m * mat(r, j); |
| 50 | } r++; |
| 51 | } return mat; } |
| 52 | matrix<T> transpose() { |
| 53 | matrix<T> res(cols, rows); |
| 54 | rep(i,0,rows) rep(j,0,cols) res(j, i) = at(i, j); |