31 void makeColumnIntegralPrimitive(
Matrix& mat,
size_t col) {
37 mpz_class numGcd = mat(0, col).get_num();
38 mpz_class denLcm = mat(0, col).get_den();
39 for (
size_t row = 1; row < mat.
getRowCount(); ++row) {
40 mpz_gcd(numGcd.get_mpz_t(), numGcd.get_mpz_t(),
41 mat(row, col).get_num_mpz_t());
42 mpz_lcm(denLcm.get_mpz_t(), denLcm.get_mpz_t(),
43 mat(row, col).get_den_mpz_t());
48 for (
size_t row = 0; row < mat.
getRowCount(); ++row) {
49 mat(row, col) /= numGcd;
50 mat(row, col) *= denLcm;
52 ASSERT(mat(row, col).get_den() == 1);
58 _rowCount(rowCount), _colCount(colCount), _entries(rowCount * colCount) {
64 Matrix mat(rowCount, colCount);
66 size_t copyRowCount = std::min(rowCount,
getRowCount());
67 size_t copyColCount = std::min(colCount,
getColCount());
68 for (
size_t row = 0; row < copyRowCount; ++row)
69 for (
size_t col = 0; col < copyColCount; ++col)
70 mat(row, col) = (*this)(row, col);
89 if (a(row, col) != b(row, col))
104 fputs(out.str().c_str(), file);
111 for (
size_t col = 0; col < mat.
getColCount(); ++col)
112 for (
size_t row = 0; row < mat.
getRowCount(); ++row)
113 pr[baseCol + col] << mat(row, col) <<
'\n';
124 prod(r, c) += a(r, i) * b(i, c);
130 if (&trans == &mat) {
136 for (
size_t row = 0; row < mat.
getRowCount(); ++row)
137 for (
size_t col = 0; col < mat.
getColCount(); ++col)
138 trans(col, row) = mat(row, col);
148 size_t sourceRow,
const mpq_class& mult) {
152 for (
size_t col = 0; col < mat.
getColCount(); ++col)
153 mat(resultRow, col) += mat(sourceRow, col) * mult;
157 for (
size_t col = 0; col < mat.
getColCount(); ++col)
158 mat(row, col) *= mult;
165 for (
size_t col = 0; col < mat.
getColCount(); ++col)
166 swap(mat(row1, col), mat(row2, col));
170 bool permutationOdd =
false;
173 for (
size_t pivotCol = 0; pivotCol < mat.
getColCount(); ++pivotCol) {
174 size_t pivotRow = rowsDone;
176 if (mat(pivotRow, pivotCol) != 0)
181 if (rowsDone != pivotRow) {
182 permutationOdd = !permutationOdd;
188 for (
size_t row = pivotRow + 1; row < mat.
getRowCount(); ++row) {
189 if (row != pivotRow && mat(row, pivotCol) != 0) {
191 -mat(row, pivotCol) / mat(pivotRow, pivotCol));
192 ASSERT(mat(row, pivotCol) == 0);
197 return permutationOdd;
208 if (mat(pivotRow, pivotCol) == 0) {
211 multiplyRow(mat, pivotRow, 1 / mat(pivotRow, pivotCol));
212 ASSERT(mat(pivotRow, pivotCol) == 1);
213 for (
size_t row = 0; row < pivotRow; ++row) {
214 if (row != pivotRow && mat(row, pivotCol) != 0) {
216 ASSERT(mat(row, pivotCol) == 0);
226 size_t rowBegin,
size_t rowEnd,
227 size_t colBegin,
size_t colEnd) {
228 ASSERT(rowBegin <= rowEnd);
230 ASSERT(colBegin <= colEnd);
235 subMatrix(tmp, mat, rowBegin, rowEnd, colBegin, colEnd);
240 sub.
resize(rowEnd - rowBegin, colEnd - colBegin);
241 for (
size_t row = rowBegin; row < rowEnd; ++row)
242 for (
size_t col = colBegin; col < colEnd; ++col)
243 sub(row - rowBegin, col - colBegin) = mat(row, col);
247 const Matrix& source,
size_t sourceRow) {
253 for (
size_t col = 0; col < colCount; ++col)
254 target(targetRow, col) = source(sourceRow, col);
264 inv.
resize(size, size + size);
265 for (
size_t i = 0; i < size; ++i)
266 inv(i, size + i) = 1;
269 if (inv(size - 1, size - 1) == 0)
272 subMatrix(inv, inv, 0, size, size, 2 * size);
285 if (mat(row, col) == 0) {
307 if (mat(row, col) == 0) {
311 colHasPivot[col] =
true;
320 for (
size_t col = 0; col < mat.
getColCount(); ++col) {
323 if (colHasPivot[col])
329 for (
size_t nullRow = 0; nullRow < basis.
getRowCount(); ++nullRow) {
330 if (colHasPivot[nullRow]) {
331 basis(nullRow, nullCol) = -mat(row, col);
333 }
else if (nullRow == col)
334 basis(nullRow, nullCol) = 1;
336 basis(nullRow, nullCol) = 0;
344 for (
size_t col = 0; col < basis.
getColCount(); ++col)
345 makeColumnIntegralPrimitive(basis, col);
355 for (
size_t col = 0; col < rhs.
getColCount(); ++col)
356 for (
size_t row = 0; row < rhs.
getRowCount(); ++row)
357 system(row, midCol + col) = rhs(row, col);
362 for (
size_t row = 0; row < system.
getRowCount(); ++row) {
363 for (
size_t col = 0; col < midCol; ++col)
364 if (system(row, col) != 0)
366 for (
size_t col = midCol; col < system.
getColCount(); ++col)
367 if (system(row, col) != 0)
376 for (
size_t col = 0; col < midCol; ++col) {
377 if (row == system.
getRowCount() || system(row, col) == 0) {
381 ASSERT(system(row, col) == 1);
383 sol(col, r) = system(row, midCol + r);
407 return solve(dummy, a, b) &&
solve(dummy, b, a);
414 bool permutationOdd =
rowReduce(reduced);
416 mpq_class det = permutationOdd ? -1 : 1;
418 det *= reduced(i, i);
423 size_t getOppositeZeroRow(
const Matrix& mat) {
433 for (
size_t opposite = 1; opposite < 4; ++opposite) {
435 for (
size_t col = 0; col < mat.
getColCount(); ++col) {
436 tmp = mat(0, col) + mat(opposite, col);
437 for (
size_t row = 1; row < 4; ++row)
439 tmp -= mat(row, col);
453 return getOppositeZeroRow(mat) != mat.
getRowCount();
458 size_t opposite = getOppositeZeroRow(mat);
461 for (a = 1; a < 4; ++a)
467 for (b = a + 1; b < 4; ++b)
474 for (
size_t col = 0; col < mat.
getColCount(); ++col) {
475 tmp(0, col) = mat(a, col) - mat(0, col);
476 tmp(1, col) = mat(b, col) - mat(0, col);