1#ifndef __LINBOX_matrix_SlicedPolynomialMatrix_SlicedPolynomialMatrix_INL 2#define __LINBOX_matrix_SlicedPolynomialMatrix_SlicedPolynomialMatrix_INL 3 4namespace Linbox 5{ 6 ////////////////////////// 7 //irreducible polynomial// 8 ////////////////////////// 9 10 template < class _Field, class _Rep, class _MatrixElement > 11 polynomial& SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::modulo(polynomial& g, polynomial&h, polynomial& f) 12 { 13 polynomial w1; 14 Poly1Dom<Domain,Dense>::div(w1, h, f); 15 polynomial w2; 16 Poly1Dom<Domain,Dense>::mul(w2, w1, f); 17 Poly1Dom<Domain,Dense>::sub(g, h, w2); 18 return g; 19 } 20 21 /* 22 algorithm description 23 http://en.wikipedia.org/wiki/Factorization_of_polynomials_over_finite_fields#Rabin.27s_test_of_irreducibility 24 Algorithm Rabin Irreducibility Test 25 Input: A monic polynomial f in Fq[x] of degree n, 26 p1, ..., pk all distinct prime divisors of n. 27 Output: Either "f is irreducible" or "f is reducible". 28 Begin 29 for j = 1 to k do 30 n_j=n/p_j; 31 for i = 1 to k do 32 h:=x^{q^{n_i}}-x \bmod f; 33 g := gcd(f, h); 34 if g ? 1, then return 'f is reducible' and STOP; 35 end for; 36 h:= x^{q^{n}}-x \bmod f; 37 if h = 0, then return "f is irreducible", 38 else return "f is reducible" 39 end. 40 */ 41 /* 42 function description 43 Fq[x], f = x^n + bx + a 44 returns true if f is irreducible and false if f is reducible 45 */ 46 template < class _Field, class _Rep, class _MatrixElement > 47 bool SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::rabin(int n, int q, int a, int b) 48 { 49 polynomial f(n + 1); 50 f[0] = a; 51 f[1] = b; 52 for (int i = 2; i < n; i++) 53 { 54 f[i] = 0; 55 } 56 f[n] = 1; 57 std::vector<int> pd; 58 factorize(n, pd); //factorizes n into primes 59 //n = pd[0]^alpha[0] * ... * pd[k - 1]^alpha[k - 1] 60 int k = pd.size(); 61 std::vector<int> nd(k); 62 for (int j = 0; j < k; j++) 63 { 64 nd[j] = n / pd[j]; 65 } 66 polynomial vector_zero(); 67 assign(vector_zero, 0); 68 polynomial vector_one(); 69 assign(vector_one, 1); 70 for (int j = 0; j < k; j++) 71 { 72 int h_size = pow(q, nd[j]) + 1; 73 polynomial h(h_size); 74 h[0] = 0; 75 h[1] = -1; 76 for (int i = 2; i < h_size - 1; i++) 77 { 78 h[i] = 0; 79 } 80 h[h_size - 1] = 1; 81 polynomial g; 82 Poly1Dom<IntField, Dense>::gcd (g, f, h); 83 bool g_equals_1 = g.areEqual(vector_one); 84 if (! g_equals_1) 85 { 86 return true; //f is irreducible 87 } 88 } 89 int h_size = pow(q, n) + 1; 90 polynomial h(h_size); 91 h[0] = 0; 92 h[1] = -1; 93 for (int i = 2; i < h_size - 1; i++) 94 { 95 h[i] = 0; 96 } 97 h[h_size - 1] = 1; 98 polynomial g; 99 modulo(g, h, f);//!!! 100 bool g_equals_0 = g.areEqual(vector_zero); 101 return g_equals_0; //if g = 0, then return "f is irreducible", else return "f is reducible" 102 } 103 104 template < class _Field, class _Rep, class _MatrixElement > 105 void SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::setIrreduciblePlynomial(int max_steps) 106 { 107 bool found = false; 108 int w = F.characteristic(); 109 for (int step = 0; (step < max_steps) && (!found); step++) 110 { 111 int a = rand() % w; 112 int b = rand() % w; 113 if (a + b > 0) 114 { 115 found = rabin(n, w, a, b); 116 } 117 } 118 irreducible.push_back(a); 119 irreducible.push_back(b); 120 for (int i = 1; i < n; i++) 121 { 122 irreducible.push_back(0); 123 } 124 irreducible.push_back(1); 125 return; 126 } 127 128 //////////////// 129 //Constructors// 130 //////////////// 131 132 template < class _Field, class _Rep, class _MatrixElement > 133 SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::SlicedPolynomialMatrix (const _Field &BF) 134 { 135 GF = BF; 136 F_temp = IntField(GF.characteristic()); //public function to set characteristic? 137 F = F_temp; 138 int n = GF.exponent(); //GF = GF(p^n) 139 for (size_t m = 0; m < n; m++) 140 { 141 V.emplace_back(BlasMatrix<IntField>(F)); 142 } 143 setIrreduciblePolynomial(); 144 } 145 146 template < class _Field, class _Rep, class _MatrixElement > 147 SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::SlicedPolynomialMatrix (const _Field &BF, const size_t & m1, const size_t &m2) 148 { 149 GF = BF; 150 F_temp = IntField(GF.characteristic()); //public function to set characteristic? 151 F = F_temp; 152 int n = GF.exponent(); //GF = GF(p^n) 153 for (size_t m = 0; m < n; m++) 154 { 155 V.emplace_back(BlasMatrix<IntField>(F, m1, m2)); 156 } 157 setIrreduciblePolynomial(); 158 } 159 160 template < class _Field, class _Rep, class _MatrixElement > 161 SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::SlicedPolynomialMatrix (const _Field &BF, polynomial& pp) 162 { 163 GF = BF; 164 F_temp = IntField(GF.characteristic()); //public function to set characteristic? 165 F = F_temp; 166 int n = GF.exponent(); //GF = GF(p^n) 167 for (size_t m = 0; m < n; m++) 168 { 169 V.emplace_back(BlasMatrix<IntField>(F)); 170 } 171 irreducible = pp; 172 } 173 174 template < class _Field, class _Rep, class _MatrixElement > 175 SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::SlicedPolynomialMatrix (const _Field &BF, const size_t & m1, const size_t &m2, polynomial& pp) 176 { 177 GF = BF; 178 F_temp = IntField(GF.characteristic()); //public function to set characteristic? 179 F = F_temp; 180 int n = GF.exponent(); //GF = GF(p^n) 181 for (size_t m = 0; m < n; m++) 182 { 183 V.emplace_back(BlasMatrix<IntField>(F, m1, m2)); 184 } 185 irreducible = pp; 186 } 187 188 /////////////// 189 // Destructor// 190 /////////////// 191 192 template < class _Field, class _Rep, class _MatrixElement > 193 SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::~SlicedPolynomialMatrix() 194 { 195 //LidiaGfq has a destructor, GivaroGfq doesn't, so currently field type members GF and F are not destroyed 196 V.~vector(); 197 //if some members are added, delete them 198 } 199 //////////////////////// 200 //dimensions of vector// 201 //////////////////////// 202 203 template < class _Field, class _Rep, class _MatrixElement > 204 size_t SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::length() const 205 { 206 return V.size(); 207 } 208 209 template < class _Field, class _Rep, class _MatrixElement > 210 size_t SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::rowdim() const 211 { 212 return V[0].rowdim(); 213 } 214 215 template < class _Field, class _Rep, class _MatrixElement > 216 size_t SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::coldim() const 217 { 218 return V[0].coldim(); 219 } 220 221 ///////////////// 222 //return fields// 223 ///////////////// 224 225 template < class _Field, class _Rep, class _MatrixElement > 226 const Field& SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::fieldGF() const 227 { 228 return GF; 229 } 230 231 template < class _Field, class _Rep, class _MatrixElement > 232 const IntField& SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::fieldF() const 233 { 234 return F; 235 } 236 237 ///////////////////////// 238 //functions for entries// 239 ///////////////////////// 240 241 template < class _Field, class _Rep, class _MatrixElement > 242 const MatrixElement& SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::setEntry (size_t m, size_t i, size_t j, const MatrixElement &a_mij) 243 { 244 V[m].setEntry(i, j, a_mij); 245 return a_mij; 246 } 247 248 template < class _Field, class _Rep, class _MatrixElement > 249 _MatrixElement & SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::refEntry (size_t m, size_t i, size_t j) 250 { 251 return V[m].refEntry(i, j); 252 253 } 254 255 template < class _Field, class _Rep, class _MatrixElement > 256 _MatrixElement & SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::getEntry (size_t m, size_t i, size_t j) 257 { 258 return V[m].getEntry(i, j); 259 260 } 261 262 ///////////////////////////////////// 263 //functions for matrix-coefficients// 264 ///////////////////////////////////// 265 266 template < class _Field, class _Rep, class _MatrixElement > 267 void SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::setMatrixCoefficient (size_t m, const BlasMatrix<IntField> &V_m) 268 { 269 V[m] = V_m; 270 } 271 272 template < class _Field, class _Rep, class _MatrixElement > 273 BlasMatrix<IntField> &SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::refMatrixCoefficient (size_t m) 274 { 275 return V[m]; 276 } 277 278 template < class _Field, class _Rep, class _MatrixElement > 279 const BlasMatrix<IntField> &SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::getMatrixCoefficient (size_t m) const 280 { 281 return V[m]; 282 } 283 284 ///////// 285 //swaps// 286 ///////// 287 288 template < class _Field, class _Rep, class _MatrixElement > 289 void SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::swapRows(size_t i1, size_t i2) 290 { 291 for (size_t m = 0; m < this->length(); m++) 292 { 293 for (size_t j = 0; j < this->coldim(); j++) 294 { 295 MatrixElement c = this->getEntry(m, i1, j); 296 this->setEntry(m, i1, j, this->getEntry(m, i2, j)); 297 this->setEntry(m, i2, j, c); 298 } 299 } 300 } 301 302 template < class _Field, class _Rep, class _MatrixElement > 303 void SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::swapCols(size_t j1, size_t j2) 304 { 305 for (size_t m = 0; m < this->length(); m++) 306 { 307 for (size_t i = 0; i < this->colrow(); i++) 308 { 309 MatrixElement c = this->getEntry(m, i, j1); 310 this->setEntry(m, i, j1, this->getEntry(m, i, j2)); 311 this->setEntry(m, i, j2, c); 312 } 313 } 314 } 315 316 ///////////// 317 //transpose// 318 ///////////// 319 320 template < class _Field, class _Rep, class _MatrixElement > 321 SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement > SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::transpose(SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement > & tV) const 322 { 323 //check dimensions 324 for (size_t m = 0; m < this->length(); m++) 325 { 326 this->getMatrixCoefficent(m).transpose(tV.refMatrixCoefficent(m)); 327 } 328 return tV; 329 } 330 331 ////////////////// 332 //input / output// 333 ////////////////// 334 template < class _Field, class _Rep, class _MatrixElement > 335 std::istream& SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::read (std::istream &file) 336 { 337 int K = this->length(); 338 int I = this->rowdim(); 339 int J = this->coldim(); 340 MatrixElement c; 341 for (int k = 0; k < K; k++) 342 { 343 for (int i = 0; i < I; i++) 344 { 345 for (int j = 0; j < J; j++) 346 { 347 file >> c; 348 this->setEntry(k, i, j, c); 349 } 350 } 351 } 352 return file; 353 } 354 355 template < class _Field, class _Rep, class _MatrixElement > 356 std::ostream& SlicedPolynomialMatrix< _Field, _Rep, _MatrixElement >::write (std::ostream &file) 357 { 358 int K = this->length(); 359 int I = this->rowdim(); 360 int J = this->coldim(); 361 for (int k = 0; k < K; k++) 362 { 363 for (int i = 0; i < I; i++) 364 { 365 for (int j = 0; j < J; j++) 366 { 367 file << this->getEntry(k, i, j) << " "; 368 } 369 file << std::endl; 370 } 371 file << std::endl; 372 } 373 return file; 374 } 375} 376 377#endif 378 379// Local Variables: 380// mode: C++ 381// tab-width: 4 382// indent-tabs-mode: nil 383// c-basic-offset: 4 384// End: 385// vim:sts=4:sw=4:ts=4:et:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s 386