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