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