1 2 /* tests/__LINBOX_smith_form_kannan_bachem.h 3 * Copyright (C) 2017 Gavin Harrison, 4 * 5 * Written by Gavin Harrison <gmh33@drexel.edu>, 6 * 7 * ========LICENCE======== 8 * This file is part of the library LinBox. 9 * 10 * LinBox is free software: you can redistribute it and/or modify 11 * it under the terms of the GNU Lesser General Public 12 * License as published by the Free Software Foundation; either 13 * version 2.1 of the License, or (at your option) any later version. 14 * 15 * This library is distributed in the hope that it will be useful, 16 * but WITHOUT ANY WARRANTY; without even the implied warranty of 17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 18 * Lesser General Public License for more details. 19 * 20 * You should have received a copy of the GNU Lesser General Public 21 * License along with this library; if not, write to the Free Software 22 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 23 * ========LICENCE======== 24 *. 25 */ 26 27 #include <iostream> 28 #include <stdlib.h> 29 #include <algorithm> 30 31 #include "linbox/matrix/densematrix/blas-matrix.h" 32 #include "linbox/matrix/matrixdomain/matrix-domain.h" 33 #include "linbox/ring/polynomial-local-x.h" 34 35 #ifndef __LINBOX_poly_smith_form_local_x_domain_H 36 #define __LINBOX_poly_smith_form_local_x_domain_H 37 38 namespace LinBox 39 { 40 // Specialized Polynomial Smith Form over Polynomials mod x^n 41 // BaseField should be a Polynomial Ring 42 template<class BaseField> 43 class PolySmithFormLocalXDomain 44 { 45 public: 46 typedef typename BaseField::Element Polynomial; 47 typedef typename BaseField::Coeff Coeff; 48 49 typedef PolynomialLocalX<BaseField> Field; 50 typedef typename Field::Element Element; 51 52 typedef MatrixDomain<Field> MatrixDom; 53 54 typedef typename MatrixDom::OwnMatrix OwnMatrix; 55 typedef typename MatrixDom::Matrix SubMatrix; 56 57 private: 58 BaseField _BF; 59 Field _F; 60 MatrixDom _MD; 61 62 public: PolySmithFormLocalXDomain(const BaseField & F,size_t e)63 PolySmithFormLocalXDomain(const BaseField &F, size_t e) : _BF(F), _F(F, e), _MD(_F) {} PolySmithFormLocalXDomain(const PolySmithFormLocalXDomain & D)64 PolySmithFormLocalXDomain(const PolySmithFormLocalXDomain &D) : _BF(D._BF), _F(D._F), _MD(_F) {} 65 66 private: 67 68 template<class Matrix1> printMatrix(const Matrix1 & A)69 void printMatrix(const Matrix1 &A) const { 70 std::cout << "[" << std::endl; 71 for (size_t i = 0; i < A.rowdim(); i++) { 72 std::cout << "\t["; 73 for (size_t j = 0; j < A.coldim(); j++) { 74 _F.write(std::cout, A.getEntry(i, j)); 75 if (j < A.coldim() - 1) { 76 std::cout << ", "; 77 } 78 } 79 std::cout << "]"; 80 if (i < A.rowdim() - 1) { 81 std::cout << ","; 82 } 83 std::cout << std::endl; 84 } 85 std::cout << "]" << std::endl; 86 } 87 88 template<typename Matrix> swapRows(Matrix & M,size_t r1,size_t r2)89 void swapRows(Matrix &M, size_t r1, size_t r2) const { 90 SubMatrix row1(M, r1, 0, 1, M.coldim()); 91 SubMatrix row2(M, r2, 0, 1, M.coldim()); 92 93 row1.swap(row2); 94 } 95 96 template<typename Matrix> swapCols(Matrix & M,size_t c1,size_t c2)97 void swapCols(Matrix &M, size_t c1, size_t c2) const { 98 SubMatrix col1(M, 0, c1, M.rowdim(), 1); 99 SubMatrix col2(M, 0, c2, M.rowdim(), 1); 100 101 col1.swap(col2); 102 } 103 104 template<class Matrix> findPivot(Matrix & A)105 bool findPivot(Matrix &A) const { 106 int d_row = -1, d_col = -1; // location of divisor 107 size_t d_firstNonZero = -1; 108 Element divisor; 109 bool d_isUnit = false; 110 111 for (size_t i = 0; i < A.rowdim() && !d_isUnit; i++) { 112 for (size_t j = 0; j < A.coldim() && !d_isUnit; j++) { 113 Element tmp; 114 _F.assign(tmp, A.getEntry(i, j)); 115 116 if (_F.isZero(tmp)) { 117 continue; 118 } 119 120 size_t firstNonZero = _F.firstNonZeroCoeff(tmp); 121 if (firstNonZero == 0) { 122 d_row = i; 123 d_col = j; 124 d_isUnit = true; 125 continue; 126 } 127 128 if (d_row > -1 && d_col > -1 && d_firstNonZero <= firstNonZero) { 129 continue; 130 } 131 132 _F.assign(divisor, tmp); 133 d_row = i; 134 d_col = j; 135 d_firstNonZero = firstNonZero; 136 } 137 } 138 139 if (d_row == -1 && d_col == -1) { 140 return false; 141 } 142 143 swapRows(A, 0, d_row); 144 swapCols(A, 0, d_col); 145 return true; 146 } 147 148 template<class Matrix> reduceMatrix(Matrix & A)149 size_t reduceMatrix(Matrix &A) { 150 Element pivot; 151 A.getEntry(pivot, 0, 0); 152 153 if (_F.isUnit(pivot)) { 154 return 0; 155 } 156 157 size_t e = _F.firstNonZeroCoeff(pivot); 158 159 for (size_t i = 0; i < A.rowdim(); i++) { 160 for (size_t j = 0; j < A.coldim(); j++) { 161 Element tmp; 162 A.getEntry(tmp, i, j); 163 _F.rightShiftIn(tmp, e); 164 A.setEntry(i, j, tmp); 165 } 166 } 167 168 _F.setExponent(_F.getExponent() - e); 169 170 return e; 171 } 172 173 template<class Matrix> zeroRow(Matrix & A)174 void zeroRow(Matrix &A) const { 175 for (size_t i = 1; i < A.coldim(); i++) { 176 A.setEntry(0, i, _F.zero); 177 } 178 } 179 180 template<class Matrix> eliminateCol(Element & pivot,Matrix & A)181 void eliminateCol(Element &pivot, Matrix &A) const { 182 Element pivot_inv; 183 A.getEntry(pivot, 0, 0); 184 185 _F.inv(pivot_inv, pivot); 186 _F.negin(pivot_inv); 187 188 SubMatrix pivotRow(A, 0, 0, 1, A.coldim()); 189 _MD.mulin(pivotRow, pivot_inv); 190 191 for (size_t r = 1; r < A.rowdim(); r++) { 192 SubMatrix otherRow(A, r, 0, 1, A.coldim()); 193 194 Element other; 195 otherRow.getEntry(other, 0, 0); 196 197 _MD.saxpyin(otherRow, other, pivotRow); 198 } 199 } 200 201 template<class Matrix> solveHelper(std::vector<Element> & L,std::vector<size_t> & es,size_t e,Matrix & A)202 void solveHelper(std::vector<Element> &L, std::vector<size_t> &es, size_t e, Matrix &A) { 203 if (A.rowdim() == 0 || A.coldim() == 0) { 204 return; 205 } 206 207 if (!findPivot(A)) { 208 size_t dim = std::min(A.rowdim(), A.coldim()); 209 for (size_t i = 0; i < dim; i++) { 210 L.push_back(_F.zero); 211 } 212 return; 213 } 214 215 e += reduceMatrix(A); 216 217 Element pivot; 218 eliminateCol(pivot, A); 219 220 zeroRow(A); 221 222 L.push_back(pivot); 223 es.push_back(e); 224 SubMatrix B(A, 1, 1, A.rowdim() - 1, A.coldim() - 1); 225 solveHelper(L, es, e, B); 226 } 227 228 public: 229 template<class Matrix> solve(std::vector<Polynomial> & result,Matrix & A)230 void solve(std::vector<Polynomial> &result, Matrix &A) { 231 size_t initial_exp = _F.getExponent(); 232 233 std::vector<Element> L; 234 std::vector<size_t> es; 235 236 solveHelper(L, es, 0, A); 237 238 _F.setExponent(initial_exp); 239 240 for (size_t i = 0; i < L.size(); i++) { 241 Element tmp; 242 _F.leftShift(tmp, L[i], es[i]); 243 244 Polynomial r; 245 _F.denormalize(r, tmp); 246 result.push_back(r); 247 } 248 } 249 250 template<class Matrix> solveDet(Polynomial & det,Matrix & A)251 void solveDet(Polynomial &det, Matrix &A) { 252 size_t initial_exp = _F.getExponent(); 253 254 std::vector<Element> L; 255 std::vector<size_t> es; 256 solveHelper(L, es, 0, A); 257 258 _F.setExponent(initial_exp); 259 260 size_t e = 0; 261 for (size_t i = 0; i < es.size(); i++) { 262 e += es[i]; 263 } 264 265 Element tmp; 266 _F.leftShift(tmp, L[0], e); 267 for (size_t i = 1; i < L.size(); i++) { 268 _F.mulin(tmp, L[i]); 269 } 270 _F.denormalize(det, tmp); 271 } 272 }; // end of class PolySmithFormLocalXDomain 273 } 274 275 #endif // __LINBOX_poly_smith_form_local_x_domain_H 276