1 | //===- Matrix.cpp - MLIR Matrix Class -------------------------------------===// |
2 | // |
3 | // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. |
4 | // See https://llvm.org/LICENSE.txt for license information. |
5 | // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception |
6 | // |
7 | //===----------------------------------------------------------------------===// |
8 | |
9 | #include "mlir/Analysis/Presburger/Matrix.h" |
10 | #include "mlir/Analysis/Presburger/Fraction.h" |
11 | #include "mlir/Analysis/Presburger/MPInt.h" |
12 | #include "mlir/Analysis/Presburger/Utils.h" |
13 | #include "mlir/Support/LLVM.h" |
14 | #include "llvm/Support/MathExtras.h" |
15 | #include "llvm/Support/raw_ostream.h" |
16 | #include <algorithm> |
17 | #include <cassert> |
18 | #include <utility> |
19 | |
20 | using namespace mlir; |
21 | using namespace presburger; |
22 | |
23 | template <typename T> |
24 | Matrix<T>::Matrix(unsigned rows, unsigned columns, unsigned reservedRows, |
25 | unsigned reservedColumns) |
26 | : nRows(rows), nColumns(columns), |
27 | nReservedColumns(std::max(a: nColumns, b: reservedColumns)), |
28 | data(nRows * nReservedColumns) { |
29 | data.reserve(std::max(a: nRows, b: reservedRows) * nReservedColumns); |
30 | } |
31 | |
32 | /// We cannot use the default implementation of operator== as it compares |
33 | /// fields like `reservedColumns` etc., which are not part of the data. |
34 | template <typename T> |
35 | bool Matrix<T>::operator==(const Matrix<T> &m) const { |
36 | if (nRows != m.getNumRows()) |
37 | return false; |
38 | if (nColumns != m.getNumColumns()) |
39 | return false; |
40 | |
41 | for (unsigned i = 0; i < nRows; i++) |
42 | if (getRow(i) != m.getRow(i)) |
43 | return false; |
44 | |
45 | return true; |
46 | } |
47 | |
48 | template <typename T> |
49 | Matrix<T> Matrix<T>::identity(unsigned dimension) { |
50 | Matrix matrix(dimension, dimension); |
51 | for (unsigned i = 0; i < dimension; ++i) |
52 | matrix(i, i) = 1; |
53 | return matrix; |
54 | } |
55 | |
56 | template <typename T> |
57 | unsigned Matrix<T>::getNumReservedRows() const { |
58 | return data.capacity() / nReservedColumns; |
59 | } |
60 | |
61 | template <typename T> |
62 | void Matrix<T>::reserveRows(unsigned rows) { |
63 | data.reserve(rows * nReservedColumns); |
64 | } |
65 | |
66 | template <typename T> |
67 | unsigned Matrix<T>::() { |
68 | resizeVertically(newNRows: nRows + 1); |
69 | return nRows - 1; |
70 | } |
71 | |
72 | template <typename T> |
73 | unsigned Matrix<T>::(ArrayRef<T> elems) { |
74 | assert(elems.size() == nColumns && "elems must match row length!" ); |
75 | unsigned row = appendExtraRow(); |
76 | for (unsigned col = 0; col < nColumns; ++col) |
77 | at(row, col) = elems[col]; |
78 | return row; |
79 | } |
80 | |
81 | template <typename T> |
82 | Matrix<T> Matrix<T>::transpose() const { |
83 | Matrix<T> transp(nColumns, nRows); |
84 | for (unsigned row = 0; row < nRows; ++row) |
85 | for (unsigned col = 0; col < nColumns; ++col) |
86 | transp(col, row) = at(row, col); |
87 | |
88 | return transp; |
89 | } |
90 | |
91 | template <typename T> |
92 | void Matrix<T>::resizeHorizontally(unsigned newNColumns) { |
93 | if (newNColumns < nColumns) |
94 | removeColumns(pos: newNColumns, count: nColumns - newNColumns); |
95 | if (newNColumns > nColumns) |
96 | insertColumns(pos: nColumns, count: newNColumns - nColumns); |
97 | } |
98 | |
99 | template <typename T> |
100 | void Matrix<T>::resize(unsigned newNRows, unsigned newNColumns) { |
101 | resizeHorizontally(newNColumns); |
102 | resizeVertically(newNRows); |
103 | } |
104 | |
105 | template <typename T> |
106 | void Matrix<T>::resizeVertically(unsigned newNRows) { |
107 | nRows = newNRows; |
108 | data.resize(nRows * nReservedColumns); |
109 | } |
110 | |
111 | template <typename T> |
112 | void Matrix<T>::swapRows(unsigned row, unsigned otherRow) { |
113 | assert((row < getNumRows() && otherRow < getNumRows()) && |
114 | "Given row out of bounds" ); |
115 | if (row == otherRow) |
116 | return; |
117 | for (unsigned col = 0; col < nColumns; col++) |
118 | std::swap(at(row, col), at(otherRow, col)); |
119 | } |
120 | |
121 | template <typename T> |
122 | void Matrix<T>::swapColumns(unsigned column, unsigned otherColumn) { |
123 | assert((column < getNumColumns() && otherColumn < getNumColumns()) && |
124 | "Given column out of bounds" ); |
125 | if (column == otherColumn) |
126 | return; |
127 | for (unsigned row = 0; row < nRows; row++) |
128 | std::swap(at(row, column), at(row, otherColumn)); |
129 | } |
130 | |
131 | template <typename T> |
132 | MutableArrayRef<T> Matrix<T>::getRow(unsigned row) { |
133 | return {&data[row * nReservedColumns], nColumns}; |
134 | } |
135 | |
136 | template <typename T> |
137 | ArrayRef<T> Matrix<T>::getRow(unsigned row) const { |
138 | return {&data[row * nReservedColumns], nColumns}; |
139 | } |
140 | |
141 | template <typename T> |
142 | void Matrix<T>::setRow(unsigned row, ArrayRef<T> elems) { |
143 | assert(elems.size() == getNumColumns() && |
144 | "elems size must match row length!" ); |
145 | for (unsigned i = 0, e = getNumColumns(); i < e; ++i) |
146 | at(row, i) = elems[i]; |
147 | } |
148 | |
149 | template <typename T> |
150 | void Matrix<T>::insertColumn(unsigned pos) { |
151 | insertColumns(pos, count: 1); |
152 | } |
153 | template <typename T> |
154 | void Matrix<T>::insertColumns(unsigned pos, unsigned count) { |
155 | if (count == 0) |
156 | return; |
157 | assert(pos <= nColumns); |
158 | unsigned oldNReservedColumns = nReservedColumns; |
159 | if (nColumns + count > nReservedColumns) { |
160 | nReservedColumns = llvm::NextPowerOf2(A: nColumns + count); |
161 | data.resize(nRows * nReservedColumns); |
162 | } |
163 | nColumns += count; |
164 | |
165 | for (int ri = nRows - 1; ri >= 0; --ri) { |
166 | for (int ci = nReservedColumns - 1; ci >= 0; --ci) { |
167 | unsigned r = ri; |
168 | unsigned c = ci; |
169 | T &dest = data[r * nReservedColumns + c]; |
170 | if (c >= nColumns) { // NOLINT |
171 | // Out of bounds columns are zero-initialized. NOLINT because clang-tidy |
172 | // complains about this branch being the same as the c >= pos one. |
173 | // |
174 | // TODO: this case can be skipped if the number of reserved columns |
175 | // didn't change. |
176 | dest = 0; |
177 | } else if (c >= pos + count) { |
178 | // Shift the data occuring after the inserted columns. |
179 | dest = data[r * oldNReservedColumns + c - count]; |
180 | } else if (c >= pos) { |
181 | // The inserted columns are also zero-initialized. |
182 | dest = 0; |
183 | } else { |
184 | // The columns before the inserted columns stay at the same (row, col) |
185 | // but this corresponds to a different location in the linearized array |
186 | // if the number of reserved columns changed. |
187 | if (nReservedColumns == oldNReservedColumns) |
188 | break; |
189 | dest = data[r * oldNReservedColumns + c]; |
190 | } |
191 | } |
192 | } |
193 | } |
194 | |
195 | template <typename T> |
196 | void Matrix<T>::removeColumn(unsigned pos) { |
197 | removeColumns(pos, count: 1); |
198 | } |
199 | template <typename T> |
200 | void Matrix<T>::removeColumns(unsigned pos, unsigned count) { |
201 | if (count == 0) |
202 | return; |
203 | assert(pos + count - 1 < nColumns); |
204 | for (unsigned r = 0; r < nRows; ++r) { |
205 | for (unsigned c = pos; c < nColumns - count; ++c) |
206 | at(r, c) = at(r, c + count); |
207 | for (unsigned c = nColumns - count; c < nColumns; ++c) |
208 | at(r, c) = 0; |
209 | } |
210 | nColumns -= count; |
211 | } |
212 | |
213 | template <typename T> |
214 | void Matrix<T>::insertRow(unsigned pos) { |
215 | insertRows(pos, count: 1); |
216 | } |
217 | template <typename T> |
218 | void Matrix<T>::insertRows(unsigned pos, unsigned count) { |
219 | if (count == 0) |
220 | return; |
221 | |
222 | assert(pos <= nRows); |
223 | resizeVertically(newNRows: nRows + count); |
224 | for (int r = nRows - 1; r >= int(pos + count); --r) |
225 | copyRow(sourceRow: r - count, targetRow: r); |
226 | for (int r = pos + count - 1; r >= int(pos); --r) |
227 | for (unsigned c = 0; c < nColumns; ++c) |
228 | at(r, c) = 0; |
229 | } |
230 | |
231 | template <typename T> |
232 | void Matrix<T>::removeRow(unsigned pos) { |
233 | removeRows(pos, count: 1); |
234 | } |
235 | template <typename T> |
236 | void Matrix<T>::removeRows(unsigned pos, unsigned count) { |
237 | if (count == 0) |
238 | return; |
239 | assert(pos + count - 1 <= nRows); |
240 | for (unsigned r = pos; r + count < nRows; ++r) |
241 | copyRow(sourceRow: r + count, targetRow: r); |
242 | resizeVertically(newNRows: nRows - count); |
243 | } |
244 | |
245 | template <typename T> |
246 | void Matrix<T>::copyRow(unsigned sourceRow, unsigned targetRow) { |
247 | if (sourceRow == targetRow) |
248 | return; |
249 | for (unsigned c = 0; c < nColumns; ++c) |
250 | at(targetRow, c) = at(sourceRow, c); |
251 | } |
252 | |
253 | template <typename T> |
254 | void Matrix<T>::fillRow(unsigned row, const T &value) { |
255 | for (unsigned col = 0; col < nColumns; ++col) |
256 | at(row, col) = value; |
257 | } |
258 | |
259 | // moveColumns is implemented by moving the columns adjacent to the source range |
260 | // to their final position. When moving right (i.e. dstPos > srcPos), the range |
261 | // of the adjacent columns is [srcPos + num, dstPos + num). When moving left |
262 | // (i.e. dstPos < srcPos) the range of the adjacent columns is [dstPos, srcPos). |
263 | // First, zeroed out columns are inserted in the final positions of the adjacent |
264 | // columns. Then, the adjacent columns are moved to their final positions by |
265 | // swapping them with the zeroed columns. Finally, the now zeroed adjacent |
266 | // columns are deleted. |
267 | template <typename T> |
268 | void Matrix<T>::moveColumns(unsigned srcPos, unsigned num, unsigned dstPos) { |
269 | if (num == 0) |
270 | return; |
271 | |
272 | int offset = dstPos - srcPos; |
273 | if (offset == 0) |
274 | return; |
275 | |
276 | assert(srcPos + num <= getNumColumns() && |
277 | "move source range exceeds matrix columns" ); |
278 | assert(dstPos + num <= getNumColumns() && |
279 | "move destination range exceeds matrix columns" ); |
280 | |
281 | unsigned insertCount = offset > 0 ? offset : -offset; |
282 | unsigned finalAdjStart = offset > 0 ? srcPos : srcPos + num; |
283 | unsigned curAdjStart = offset > 0 ? srcPos + num : dstPos; |
284 | // TODO: This can be done using std::rotate. |
285 | // Insert new zero columns in the positions where the adjacent columns are to |
286 | // be moved. |
287 | insertColumns(pos: finalAdjStart, count: insertCount); |
288 | // Update curAdjStart if insertion of new columns invalidates it. |
289 | if (finalAdjStart < curAdjStart) |
290 | curAdjStart += insertCount; |
291 | |
292 | // Swap the adjacent columns with inserted zero columns. |
293 | for (unsigned i = 0; i < insertCount; ++i) |
294 | swapColumns(column: finalAdjStart + i, otherColumn: curAdjStart + i); |
295 | |
296 | // Delete the now redundant zero columns. |
297 | removeColumns(pos: curAdjStart, count: insertCount); |
298 | } |
299 | |
300 | template <typename T> |
301 | void Matrix<T>::addToRow(unsigned sourceRow, unsigned targetRow, |
302 | const T &scale) { |
303 | addToRow(targetRow, getRow(sourceRow), scale); |
304 | } |
305 | |
306 | template <typename T> |
307 | void Matrix<T>::addToRow(unsigned row, ArrayRef<T> rowVec, const T &scale) { |
308 | if (scale == 0) |
309 | return; |
310 | for (unsigned col = 0; col < nColumns; ++col) |
311 | at(row, col) += scale * rowVec[col]; |
312 | } |
313 | |
314 | template <typename T> |
315 | void Matrix<T>::scaleRow(unsigned row, const T &scale) { |
316 | for (unsigned col = 0; col < nColumns; ++col) |
317 | at(row, col) *= scale; |
318 | } |
319 | |
320 | template <typename T> |
321 | void Matrix<T>::addToColumn(unsigned sourceColumn, unsigned targetColumn, |
322 | const T &scale) { |
323 | if (scale == 0) |
324 | return; |
325 | for (unsigned row = 0, e = getNumRows(); row < e; ++row) |
326 | at(row, targetColumn) += scale * at(row, sourceColumn); |
327 | } |
328 | |
329 | template <typename T> |
330 | void Matrix<T>::negateColumn(unsigned column) { |
331 | for (unsigned row = 0, e = getNumRows(); row < e; ++row) |
332 | at(row, column) = -at(row, column); |
333 | } |
334 | |
335 | template <typename T> |
336 | void Matrix<T>::negateRow(unsigned row) { |
337 | for (unsigned column = 0, e = getNumColumns(); column < e; ++column) |
338 | at(row, column) = -at(row, column); |
339 | } |
340 | |
341 | template <typename T> |
342 | void Matrix<T>::negateMatrix() { |
343 | for (unsigned row = 0; row < nRows; ++row) |
344 | negateRow(row); |
345 | } |
346 | |
347 | template <typename T> |
348 | SmallVector<T, 8> Matrix<T>::preMultiplyWithRow(ArrayRef<T> rowVec) const { |
349 | assert(rowVec.size() == getNumRows() && "Invalid row vector dimension!" ); |
350 | |
351 | SmallVector<T, 8> result(getNumColumns(), T(0)); |
352 | for (unsigned col = 0, e = getNumColumns(); col < e; ++col) |
353 | for (unsigned i = 0, e = getNumRows(); i < e; ++i) |
354 | result[col] += rowVec[i] * at(i, col); |
355 | return result; |
356 | } |
357 | |
358 | template <typename T> |
359 | SmallVector<T, 8> Matrix<T>::postMultiplyWithColumn(ArrayRef<T> colVec) const { |
360 | assert(getNumColumns() == colVec.size() && |
361 | "Invalid column vector dimension!" ); |
362 | |
363 | SmallVector<T, 8> result(getNumRows(), T(0)); |
364 | for (unsigned row = 0, e = getNumRows(); row < e; row++) |
365 | for (unsigned i = 0, e = getNumColumns(); i < e; i++) |
366 | result[row] += at(row, i) * colVec[i]; |
367 | return result; |
368 | } |
369 | |
370 | /// Set M(row, targetCol) to its remainder on division by M(row, sourceCol) |
371 | /// by subtracting from column targetCol an appropriate integer multiple of |
372 | /// sourceCol. This brings M(row, targetCol) to the range [0, M(row, |
373 | /// sourceCol)). Apply the same column operation to otherMatrix, with the same |
374 | /// integer multiple. |
375 | static void modEntryColumnOperation(Matrix<MPInt> &m, unsigned row, |
376 | unsigned sourceCol, unsigned targetCol, |
377 | Matrix<MPInt> &otherMatrix) { |
378 | assert(m(row, sourceCol) != 0 && "Cannot divide by zero!" ); |
379 | assert(m(row, sourceCol) > 0 && "Source must be positive!" ); |
380 | MPInt ratio = -floorDiv(lhs: m(row, targetCol), rhs: m(row, sourceCol)); |
381 | m.addToColumn(sourceColumn: sourceCol, targetColumn: targetCol, scale: ratio); |
382 | otherMatrix.addToColumn(sourceColumn: sourceCol, targetColumn: targetCol, scale: ratio); |
383 | } |
384 | |
385 | template <typename T> |
386 | Matrix<T> Matrix<T>::getSubMatrix(unsigned fromRow, unsigned toRow, |
387 | unsigned fromColumn, |
388 | unsigned toColumn) const { |
389 | assert(fromRow <= toRow && "end of row range must be after beginning!" ); |
390 | assert(toRow < nRows && "end of row range out of bounds!" ); |
391 | assert(fromColumn <= toColumn && |
392 | "end of column range must be after beginning!" ); |
393 | assert(toColumn < nColumns && "end of column range out of bounds!" ); |
394 | Matrix<T> subMatrix(toRow - fromRow + 1, toColumn - fromColumn + 1); |
395 | for (unsigned i = fromRow; i <= toRow; ++i) |
396 | for (unsigned j = fromColumn; j <= toColumn; ++j) |
397 | subMatrix(i - fromRow, j - fromColumn) = at(i, j); |
398 | return subMatrix; |
399 | } |
400 | |
401 | template <typename T> |
402 | void Matrix<T>::print(raw_ostream &os) const { |
403 | for (unsigned row = 0; row < nRows; ++row) { |
404 | for (unsigned column = 0; column < nColumns; ++column) |
405 | os << at(row, column) << ' '; |
406 | os << '\n'; |
407 | } |
408 | } |
409 | |
410 | /// We iterate over the `indicator` bitset, checking each bit. If a bit is 1, |
411 | /// we append it to one matrix, and if it is zero, we append it to the other. |
412 | template <typename T> |
413 | std::pair<Matrix<T>, Matrix<T>> |
414 | Matrix<T>::splitByBitset(ArrayRef<int> indicator) { |
415 | Matrix<T> rowsForOne(0, nColumns), rowsForZero(0, nColumns); |
416 | for (unsigned i = 0; i < nRows; i++) { |
417 | if (indicator[i] == 1) |
418 | rowsForOne.appendExtraRow(getRow(i)); |
419 | else |
420 | rowsForZero.appendExtraRow(getRow(i)); |
421 | } |
422 | return {rowsForOne, rowsForZero}; |
423 | } |
424 | |
425 | template <typename T> |
426 | void Matrix<T>::dump() const { |
427 | print(os&: llvm::errs()); |
428 | } |
429 | |
430 | template <typename T> |
431 | bool Matrix<T>::hasConsistentState() const { |
432 | if (data.size() != nRows * nReservedColumns) |
433 | return false; |
434 | if (nColumns > nReservedColumns) |
435 | return false; |
436 | #ifdef EXPENSIVE_CHECKS |
437 | for (unsigned r = 0; r < nRows; ++r) |
438 | for (unsigned c = nColumns; c < nReservedColumns; ++c) |
439 | if (data[r * nReservedColumns + c] != 0) |
440 | return false; |
441 | #endif |
442 | return true; |
443 | } |
444 | |
445 | namespace mlir { |
446 | namespace presburger { |
447 | template class Matrix<MPInt>; |
448 | template class Matrix<Fraction>; |
449 | } // namespace presburger |
450 | } // namespace mlir |
451 | |
452 | IntMatrix IntMatrix::identity(unsigned dimension) { |
453 | IntMatrix matrix(dimension, dimension); |
454 | for (unsigned i = 0; i < dimension; ++i) |
455 | matrix(i, i) = 1; |
456 | return matrix; |
457 | } |
458 | |
459 | std::pair<IntMatrix, IntMatrix> IntMatrix::computeHermiteNormalForm() const { |
460 | // We start with u as an identity matrix and perform operations on h until h |
461 | // is in hermite normal form. We apply the same sequence of operations on u to |
462 | // obtain a transform that takes h to hermite normal form. |
463 | IntMatrix h = *this; |
464 | IntMatrix u = IntMatrix::identity(dimension: h.getNumColumns()); |
465 | |
466 | unsigned echelonCol = 0; |
467 | // Invariant: in all rows above row, all columns from echelonCol onwards |
468 | // are all zero elements. In an iteration, if the curent row has any non-zero |
469 | // elements echelonCol onwards, we bring one to echelonCol and use it to |
470 | // make all elements echelonCol + 1 onwards zero. |
471 | for (unsigned row = 0; row < h.getNumRows(); ++row) { |
472 | // Search row for a non-empty entry, starting at echelonCol. |
473 | unsigned nonZeroCol = echelonCol; |
474 | for (unsigned e = h.getNumColumns(); nonZeroCol < e; ++nonZeroCol) { |
475 | if (h(row, nonZeroCol) == 0) |
476 | continue; |
477 | break; |
478 | } |
479 | |
480 | // Continue to the next row with the same echelonCol if this row is all |
481 | // zeros from echelonCol onwards. |
482 | if (nonZeroCol == h.getNumColumns()) |
483 | continue; |
484 | |
485 | // Bring the non-zero column to echelonCol. This doesn't affect rows |
486 | // above since they are all zero at these columns. |
487 | if (nonZeroCol != echelonCol) { |
488 | h.swapColumns(column: nonZeroCol, otherColumn: echelonCol); |
489 | u.swapColumns(column: nonZeroCol, otherColumn: echelonCol); |
490 | } |
491 | |
492 | // Make h(row, echelonCol) non-negative. |
493 | if (h(row, echelonCol) < 0) { |
494 | h.negateColumn(column: echelonCol); |
495 | u.negateColumn(column: echelonCol); |
496 | } |
497 | |
498 | // Make all the entries in row after echelonCol zero. |
499 | for (unsigned i = echelonCol + 1, e = h.getNumColumns(); i < e; ++i) { |
500 | // We make h(row, i) non-negative, and then apply the Euclidean GCD |
501 | // algorithm to (row, i) and (row, echelonCol). At the end, one of them |
502 | // has value equal to the gcd of the two entries, and the other is zero. |
503 | |
504 | if (h(row, i) < 0) { |
505 | h.negateColumn(column: i); |
506 | u.negateColumn(column: i); |
507 | } |
508 | |
509 | unsigned targetCol = i, sourceCol = echelonCol; |
510 | // At every step, we set h(row, targetCol) %= h(row, sourceCol), and |
511 | // swap the indices sourceCol and targetCol. (not the columns themselves) |
512 | // This modulo is implemented as a subtraction |
513 | // h(row, targetCol) -= quotient * h(row, sourceCol), |
514 | // where quotient = floor(h(row, targetCol) / h(row, sourceCol)), |
515 | // which brings h(row, targetCol) to the range [0, h(row, sourceCol)). |
516 | // |
517 | // We are only allowed column operations; we perform the above |
518 | // for every row, i.e., the above subtraction is done as a column |
519 | // operation. This does not affect any rows above us since they are |
520 | // guaranteed to be zero at these columns. |
521 | while (h(row, targetCol) != 0 && h(row, sourceCol) != 0) { |
522 | modEntryColumnOperation(m&: h, row, sourceCol, targetCol, otherMatrix&: u); |
523 | std::swap(a&: targetCol, b&: sourceCol); |
524 | } |
525 | |
526 | // One of (row, echelonCol) and (row, i) is zero and the other is the gcd. |
527 | // Make it so that (row, echelonCol) holds the non-zero value. |
528 | if (h(row, echelonCol) == 0) { |
529 | h.swapColumns(column: i, otherColumn: echelonCol); |
530 | u.swapColumns(column: i, otherColumn: echelonCol); |
531 | } |
532 | } |
533 | |
534 | // Make all entries before echelonCol non-negative and strictly smaller |
535 | // than the pivot entry. |
536 | for (unsigned i = 0; i < echelonCol; ++i) |
537 | modEntryColumnOperation(m&: h, row, sourceCol: echelonCol, targetCol: i, otherMatrix&: u); |
538 | |
539 | ++echelonCol; |
540 | } |
541 | |
542 | return {h, u}; |
543 | } |
544 | |
545 | MPInt IntMatrix::normalizeRow(unsigned row, unsigned cols) { |
546 | return normalizeRange(range: getRow(row).slice(N: 0, M: cols)); |
547 | } |
548 | |
549 | MPInt IntMatrix::normalizeRow(unsigned row) { |
550 | return normalizeRow(row, cols: getNumColumns()); |
551 | } |
552 | |
553 | MPInt IntMatrix::determinant(IntMatrix *inverse) const { |
554 | assert(nRows == nColumns && |
555 | "determinant can only be calculated for square matrices!" ); |
556 | |
557 | FracMatrix m(*this); |
558 | |
559 | FracMatrix fracInverse(nRows, nColumns); |
560 | MPInt detM = m.determinant(inverse: &fracInverse).getAsInteger(); |
561 | |
562 | if (detM == 0) |
563 | return MPInt(0); |
564 | |
565 | if (!inverse) |
566 | return detM; |
567 | |
568 | *inverse = IntMatrix(nRows, nColumns); |
569 | for (unsigned i = 0; i < nRows; i++) |
570 | for (unsigned j = 0; j < nColumns; j++) |
571 | inverse->at(row: i, column: j) = (fracInverse.at(row: i, column: j) * detM).getAsInteger(); |
572 | |
573 | return detM; |
574 | } |
575 | |
576 | FracMatrix FracMatrix::identity(unsigned dimension) { |
577 | return Matrix::identity(dimension); |
578 | } |
579 | |
580 | FracMatrix::FracMatrix(IntMatrix m) |
581 | : FracMatrix(m.getNumRows(), m.getNumColumns()) { |
582 | for (unsigned i = 0, r = m.getNumRows(); i < r; i++) |
583 | for (unsigned j = 0, c = m.getNumColumns(); j < c; j++) |
584 | this->at(row: i, column: j) = m.at(row: i, column: j); |
585 | } |
586 | |
587 | Fraction FracMatrix::determinant(FracMatrix *inverse) const { |
588 | assert(nRows == nColumns && |
589 | "determinant can only be calculated for square matrices!" ); |
590 | |
591 | FracMatrix m(*this); |
592 | FracMatrix tempInv(nRows, nColumns); |
593 | if (inverse) |
594 | tempInv = FracMatrix::identity(dimension: nRows); |
595 | |
596 | Fraction a, b; |
597 | // Make the matrix into upper triangular form using |
598 | // gaussian elimination with row operations. |
599 | // If inverse is required, we apply more operations |
600 | // to turn the matrix into diagonal form. We apply |
601 | // the same operations to the inverse matrix, |
602 | // which is initially identity. |
603 | // Either way, the product of the diagonal elements |
604 | // is then the determinant. |
605 | for (unsigned i = 0; i < nRows; i++) { |
606 | if (m(i, i) == 0) |
607 | // First ensure that the diagonal |
608 | // element is nonzero, by swapping |
609 | // it with a nonzero row. |
610 | for (unsigned j = i + 1; j < nRows; j++) { |
611 | if (m(j, i) != 0) { |
612 | m.swapRows(row: j, otherRow: i); |
613 | if (inverse) |
614 | tempInv.swapRows(row: j, otherRow: i); |
615 | break; |
616 | } |
617 | } |
618 | |
619 | b = m.at(row: i, column: i); |
620 | if (b == 0) |
621 | return 0; |
622 | |
623 | // Set all elements above the |
624 | // diagonal to zero. |
625 | if (inverse) { |
626 | for (unsigned j = 0; j < i; j++) { |
627 | if (m.at(row: j, column: i) == 0) |
628 | continue; |
629 | a = m.at(row: j, column: i); |
630 | // Set element (j, i) to zero |
631 | // by subtracting the ith row, |
632 | // appropriately scaled. |
633 | m.addToRow(sourceRow: i, targetRow: j, scale: -a / b); |
634 | tempInv.addToRow(sourceRow: i, targetRow: j, scale: -a / b); |
635 | } |
636 | } |
637 | |
638 | // Set all elements below the |
639 | // diagonal to zero. |
640 | for (unsigned j = i + 1; j < nRows; j++) { |
641 | if (m.at(row: j, column: i) == 0) |
642 | continue; |
643 | a = m.at(row: j, column: i); |
644 | // Set element (j, i) to zero |
645 | // by subtracting the ith row, |
646 | // appropriately scaled. |
647 | m.addToRow(sourceRow: i, targetRow: j, scale: -a / b); |
648 | if (inverse) |
649 | tempInv.addToRow(sourceRow: i, targetRow: j, scale: -a / b); |
650 | } |
651 | } |
652 | |
653 | // Now only diagonal elements of m are nonzero, but they are |
654 | // not necessarily 1. To get the true inverse, we should |
655 | // normalize them and apply the same scale to the inverse matrix. |
656 | // For efficiency we skip scaling m and just scale tempInv appropriately. |
657 | if (inverse) { |
658 | for (unsigned i = 0; i < nRows; i++) |
659 | for (unsigned j = 0; j < nRows; j++) |
660 | tempInv.at(row: i, column: j) = tempInv.at(row: i, column: j) / m(i, i); |
661 | |
662 | *inverse = std::move(tempInv); |
663 | } |
664 | |
665 | Fraction determinant = 1; |
666 | for (unsigned i = 0; i < nRows; i++) |
667 | determinant *= m.at(row: i, column: i); |
668 | |
669 | return determinant; |
670 | } |
671 | |
672 | FracMatrix FracMatrix::gramSchmidt() const { |
673 | // Create a copy of the argument to store |
674 | // the orthogonalised version. |
675 | FracMatrix orth(*this); |
676 | |
677 | // For each vector (row) in the matrix, subtract its unit |
678 | // projection along each of the previous vectors. |
679 | // This ensures that it has no component in the direction |
680 | // of any of the previous vectors. |
681 | for (unsigned i = 1, e = getNumRows(); i < e; i++) { |
682 | for (unsigned j = 0; j < i; j++) { |
683 | Fraction jNormSquared = dotProduct(a: orth.getRow(row: j), b: orth.getRow(row: j)); |
684 | assert(jNormSquared != 0 && "some row became zero! Inputs to this " |
685 | "function must be linearly independent." ); |
686 | Fraction projectionScale = |
687 | dotProduct(a: orth.getRow(row: i), b: orth.getRow(row: j)) / jNormSquared; |
688 | orth.addToRow(sourceRow: j, targetRow: i, scale: -projectionScale); |
689 | } |
690 | } |
691 | return orth; |
692 | } |
693 | |
694 | // Convert the matrix, interpreted (row-wise) as a basis |
695 | // to an LLL-reduced basis. |
696 | // |
697 | // This is an implementation of the algorithm described in |
698 | // "Factoring polynomials with rational coefficients" by |
699 | // A. K. Lenstra, H. W. Lenstra Jr., L. Lovasz. |
700 | // |
701 | // Let {b_1, ..., b_n} be the current basis and |
702 | // {b_1*, ..., b_n*} be the Gram-Schmidt orthogonalised |
703 | // basis (unnormalized). |
704 | // Define the Gram-Schmidt coefficients μ_ij as |
705 | // (b_i • b_j*) / (b_j* • b_j*), where (•) represents the inner product. |
706 | // |
707 | // We iterate starting from the second row to the last row. |
708 | // |
709 | // For the kth row, we first check μ_kj for all rows j < k. |
710 | // We subtract b_j (scaled by the integer nearest to μ_kj) |
711 | // from b_k. |
712 | // |
713 | // Now, we update k. |
714 | // If b_k and b_{k-1} satisfy the Lovasz condition |
715 | // |b_k|^2 ≥ (δ - μ_k{k-1}^2) |b_{k-1}|^2, |
716 | // we are done and we increment k. |
717 | // Otherwise, we swap b_k and b_{k-1} and decrement k. |
718 | // |
719 | // We repeat this until k = n and return. |
720 | void FracMatrix::LLL(Fraction delta) { |
721 | MPInt nearest; |
722 | Fraction mu; |
723 | |
724 | // `gsOrth` holds the Gram-Schmidt orthogonalisation |
725 | // of the matrix at all times. It is recomputed every |
726 | // time the matrix is modified during the algorithm. |
727 | // This is naive and can be optimised. |
728 | FracMatrix gsOrth = gramSchmidt(); |
729 | |
730 | // We start from the second row. |
731 | unsigned k = 1; |
732 | while (k < getNumRows()) { |
733 | for (unsigned j = k - 1; j < k; j--) { |
734 | // Compute the Gram-Schmidt coefficient μ_jk. |
735 | mu = dotProduct(a: getRow(row: k), b: gsOrth.getRow(row: j)) / |
736 | dotProduct(a: gsOrth.getRow(row: j), b: gsOrth.getRow(row: j)); |
737 | nearest = round(f: mu); |
738 | // Subtract b_j scaled by the integer nearest to μ_jk from b_k. |
739 | addToRow(row: k, rowVec: getRow(row: j), scale: -Fraction(nearest, 1)); |
740 | gsOrth = gramSchmidt(); // Update orthogonalization. |
741 | } |
742 | mu = dotProduct(a: getRow(row: k), b: gsOrth.getRow(row: k - 1)) / |
743 | dotProduct(a: gsOrth.getRow(row: k - 1), b: gsOrth.getRow(row: k - 1)); |
744 | // Check the Lovasz condition for b_k and b_{k-1}. |
745 | if (dotProduct(a: gsOrth.getRow(row: k), b: gsOrth.getRow(row: k)) > |
746 | (delta - mu * mu) * |
747 | dotProduct(a: gsOrth.getRow(row: k - 1), b: gsOrth.getRow(row: k - 1))) { |
748 | // If it is satisfied, proceed to the next k. |
749 | k += 1; |
750 | } else { |
751 | // If it is not satisfied, decrement k (without |
752 | // going beyond the second row). |
753 | swapRows(row: k, otherRow: k - 1); |
754 | gsOrth = gramSchmidt(); // Update orthogonalization. |
755 | k = k > 1 ? k - 1 : 1; |
756 | } |
757 | } |
758 | } |
759 | |
760 | IntMatrix FracMatrix::normalizeRows() const { |
761 | unsigned numRows = getNumRows(); |
762 | unsigned numColumns = getNumColumns(); |
763 | IntMatrix normalized(numRows, numColumns); |
764 | |
765 | MPInt lcmDenoms = MPInt(1); |
766 | for (unsigned i = 0; i < numRows; i++) { |
767 | // For a row, first compute the LCM of the denominators. |
768 | for (unsigned j = 0; j < numColumns; j++) |
769 | lcmDenoms = lcm(a: lcmDenoms, b: at(row: i, column: j).den); |
770 | // Then, multiply by it throughout and convert to integers. |
771 | for (unsigned j = 0; j < numColumns; j++) |
772 | normalized(i, j) = (at(row: i, column: j) * lcmDenoms).getAsInteger(); |
773 | } |
774 | return normalized; |
775 | } |
776 | |