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