Context Navigation

source: trunk/linbox/linbox/algorithms/charpoly-rational.h @ 4086

Revision 4086, 9.4 KB checked in by bboyer, 19 months ago (diff)
clang++ emits no error (almost) and compiles linbox \!

Line
1	/* -- mode: C++; tab-width: 8; indent-tabs-mode: t; c-basic-offset: 8 -- */
2	// vim:sts=8:sw=8:ts=8:noet:sr:cino=>s,f0,{0,g0,(0,\:0,t0,+0,=s
3	/* linbox/blackbox/rational-reconstruction-base.h
4	* Copyright (C) 2009 Anna Marszalek
5	*
6	* Written by Anna Marszalek <aniau@astronet.pl>
7	*
8	* This library is free software; you can redistribute it and/or
9	* modify it under the terms of the GNU Lesser General Public
10	* License as published by the Free Software Foundation; either
11	* version 2 of the License, or (at your option) any later version.
12	*
13	* This library is distributed in the hope that it will be useful,
14	* but WITHOUT ANY WARRANTY; without even the implied warranty of
15	* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16	* Lesser General Public License for more details.
17	*
18	* You should have received a copy of the GNU Lesser General Public
19	* License along with this library; if not, write to the
20	* Free Software Foundation, Inc., 51 Franklin Street - Fifth Floor,
21	* Boston, MA 02110-1301, USA.
22	*/
23
24	#ifndef __LINBOX_charpoly_rational_H
25	#define __LINBOX_charpoly_rational_H
26
27	#include "linbox/util/commentator.h"
28	#include "linbox/util/timer.h"
29	#include "linbox/field/modular.h"
30
31	//#include "linbox/field/gmp-rational.h"
32	#include "linbox/field/PID-integer.h"
33	#include "linbox/blackbox/rational-matrix-factory.h"
34	#include "linbox/algorithms/cra-early-multip.h"
35	#include "linbox/algorithms/cra-domain.h"
36	//#include "linbox/algorithms/rational-cra.h"
37	#include "linbox/algorithms/rational-reconstruction-base.h"
38	#include "linbox/algorithms/classic-rational-reconstruction.h"
39	#include "linbox/solutions/charpoly.h"
40	#include "linbox/blackbox/compose.h"
41	#include "linbox/blackbox/diagonal.h"
42
43	namespace LinBox
44	{
45	//typedef GMPRationalField Rationals;
46	//typedef Rationals::Element Quotient;
47
48	/*
49	* Computes the characteristic polynomial of a rational dense matrix
50	*/
51
52	template<class T1, class T2>
53	struct MyModularCharpoly{
54	T1* t1;
55	T2* t2;
56
57	int switcher;
58
59	MyModularCharpoly(T1* s1, T2* s2, int s = 1) {t1=s1; t2=s2;switcher = s;}
60
61	int setSwitcher(int s) {return switcher = s;}
62
63	template<typename Polynomial, typename Field>
64	Polynomial& operator()(Polynomial& P, const Field& F) const
65	{
66	if (switcher ==1) {
67	t1->operator()(P,F);
68	}
69	else {
70	t2->operator()(P,F);
71	}
72	return P;
73	}
74	};
75
76	template <class Blackbox, class MyMethod>
77	struct MyRationalModularCharpoly {
78	const Blackbox &A;
79	const MyMethod &M;
80	const std::vector<Integer> &mul;//multiplicative prec;
81
82	MyRationalModularCharpoly(const Blackbox& b, const MyMethod& n,
83	const std::vector<Integer >& p) :
84	A(b), M(n), mul(p)
85	{}
86	MyRationalModularCharpoly(MyRationalModularCharpoly& C) :
87	MyRationalModularCharpoly(C.A,C.M,C.mul)
88	{}
89
90	template<typename Polynomial, typename Field>
91	Polynomial& operator()(Polynomial& P, const Field& F) const
92	{
93	typedef typename Blackbox::template rebind<Field>::other FBlackbox;
94	FBlackbox * Ap;
95	MatrixHom::map(Ap, A, F);
96	charpoly( P, *Ap, typename FieldTraits<Field>::categoryTag(), M);
97	typename std::vector<Integer >::const_iterator it = mul.begin();
98	typename Polynomial::iterator it_p = P.begin();
99	for (;it_p !=P.end(); ++it, ++it_p) {
100	typename Field::Element e;
101	F.init(e, *it);
102	F.mulin(*it_p,e);
103	}
104
105	delete Ap;
106	return P;
107	}
108	};
109
110	template <class Blackbox, class MyMethod>
111	struct MyIntegerModularCharpoly {
112	const Blackbox &A;
113	const MyMethod &M;
114	const std::vector<typename Blackbox::Field::Element> &vD;//diagonal div. prec;
115	const std::vector<typename Blackbox::Field::Element > &mul;//multiplicative prec;
116
117	MyIntegerModularCharpoly(const Blackbox& b, const MyMethod& n,
118	const std::vector<typename Blackbox::Field::Element>& ve,
119	const std::vector<typename Blackbox::Field::Element >& p) :
120	A(b), M(n), vD(ve), mul(p) {}
121
122	MyIntegerModularCharpoly(MyIntegerModularCharpoly& C) :
123	MyIntegerModularCharpoly(C.A,C.M,C.vD,C.mul)
124	{}
125
126	template<typename Polynomial, typename Field>
127	Polynomial& operator()(Polynomial& P, const Field& F) const
128	{
129	typedef typename Blackbox::template rebind<Field>::other FBlackbox;
130	FBlackbox * Ap;
131	MatrixHom::map(Ap, A, F);
132
133	typename std::vector<typename Blackbox::Field::Element>::const_iterator it;
134
135	int i=0;
136	for (it = vD.begin(); it != vD.end(); ++it,++i) {
137	typename Field::Element t,tt;
138	F.init(t,*it);
139	F.invin(t);
140	for (int j=0; j < A.coldim(); ++j) {
141	F.mulin(Ap->refEntry(i,j),t);
142	}
143	}
144
145	charpoly( P, *Ap, typename FieldTraits<Field>::categoryTag(), M);
146	typename std::vector<typename Blackbox::Field::Element >::const_iterator it2 = mul.begin();
147	typename Polynomial::iterator it_p = P.begin();
148	for (;it_p !=P.end(); ++it2, ++it_p) {
149	typename Field::Element e;
150	F.init(e, *it2);
151	F.mulin(*it_p,e);
152	}
153
154	delete Ap;
155	return P;
156	}
157	};
158
159	template <class Rationals, template <class> class Vector, class MyMethod >
160	Vector<typename Rationals::Element>& rational_charpoly (Vector<typename Rationals::Element> &p,
161	const BlasMatrix<Rationals > &A,
162	const MyMethod &Met= Method::Hybrid())
163	{
164
165	typedef Modular<double> myModular;
166	typedef typename Rationals::Element Quotient;
167
168	commentator.start ("Rational Charpoly", "Rminpoly");
169
170	RandomPrimeIterator genprime( 26-(int)ceil(log((double)A.rowdim())*0.7213475205));
171
172	std::vector<Integer> F(A.rowdim()+1,1);
173	std::vector<Integer> M(A.rowdim()+1,1);
174	std::vector<Integer> Di(A.rowdim());
175
176	RationalMatrixFactory<PID_integer,Rationals,BlasMatrix<Rationals > > FA(&A);
177	Integer da=1, di=1; Integer D=1;
178	FA.denominator(da);
179
180	for (int i=M.size()-2; i >= 0 ; --i) {
181	//c[m]=1, c[0]=det(A);
182	FA.denominator(di,i);
183	D *=di;
184	Di[i]=di;
185	M[i] = M[i+1]*da;
186	}
187	for (int i=0; i < M.size() ; ++i ) {
188	gcd(M[i],M[i],D);
189	}
190
191	PID_integer Z;
192	BlasMatrix<PID_integer> Atilde(Z,A.rowdim(), A.coldim());
193	FA.makeAtilde(Atilde);
194
195	ChineseRemainder< EarlyMultipCRA<Modular<double> > > cra(4UL);
196	MyRationalModularCharpoly<BlasMatrix<Rationals > , MyMethod> iteration1(A, Met, M);
197	MyIntegerModularCharpoly<BlasMatrix<PID_integer>, MyMethod> iteration2(Atilde, Met, Di, M);
198	MyModularCharpoly<MyRationalModularCharpoly<BlasMatrix<Rationals > , MyMethod>,
199	MyIntegerModularCharpoly<BlasMatrix<PID_integer>, MyMethod> > iteration(&iteration1,&iteration2);
200
201	RReconstruction<PID_integer, ClassicMaxQRationalReconstruction<PID_integer> > RR;
202
203	std::vector<Integer> PP; // use of integer due to non genericity of cra. PG 2005-08-04
204	UserTimer t1,t2;
205	t1.clear();
206	t2.clear();
207	t1.start();
208	cra(2,PP,iteration1,genprime);
209	t1.stop();
210	t2.start();
211	cra(2,PP,iteration2,genprime);
212	t2.stop();
213
214
215
216	if (t1.time() < t2.time()) {
217	//cout << "ratim";
218	iteration.setSwitcher(1);
219	}
220	else {
221	//cout << "intim";
222	iteration.setSwitcher(2);
223	}
224
225	int k=4;
226	while (! cra(k,PP, iteration, genprime)) {
227	k *=2;
228	Integer m; //Integer r; Integer a,b;
229	cra.getModulus(m);
230	cra.result(PP);//need to divide
231	for (int i=0; i < PP.size(); ++i) {
232	Integer D_1;
233	inv(D_1,M[i],m);
234	PP[i] = (PP[i]*D_1) % m;
235	}
236	Integer den,den1;
237	std::vector<Integer> num(A.rowdim()+1);
238	std::vector<Integer> num1(A.rowdim()+1);
239	if (RR.reconstructRational(num,den,PP,m,-1)) {//performs reconstruction strating form c[m], use c[i] as prec for c[i-1]
240	cra(1,PP,iteration,genprime);
241	cra.getModulus(m);
242	for (int i=0; i < PP.size(); ++i) {
243	Integer D_1;
244	inv(D_1,M[i],m);
245	PP[i] = (PP[i]*D_1) % m;
246	}
247	bool terminated = true;
248	if (RR.reconstructRational(num1,den1,PP,m,-1)) {
249	if (den==den1) {
250	for (int i=0; i < num.size(); ++i) {
251	if (num[i] != num1[i]) {
252	terminated =false;
253	break;
254	}
255	}
256	}
257	else {
258	terminated = false;
259	}
260	//set p
261	if (terminated) {
262	size_t i =0;
263	integer t,tt,ttt;
264	integer err;
265	// size_t max_err = 0;
266	Quotient qerr;
267	p.resize(PP.size());
268	typename Vector <typename Rationals::Element>::iterator it;
269	Rationals Q;
270	for (it= p.begin(); it != p.end(); ++it, ++i) {
271	A.field().init(*it, num[i],den);
272	Q.get_den(t,*it);
273	if (it != p.begin()) Q.get_den(tt,*(it-1));
274	else tt = 1;
275	Q.init(qerr,t,tt);
276
277	}
278	return p;
279	break;
280	}
281	}
282	}
283	}
284
285	cra.result(PP);
286
287	size_t i =0;
288	integer t,tt;
289	integer err;
290	size_t max_res=0;int max_i;
291	// double rel;
292	size_t max_resu=0; int max_iu;
293	// size_t max_err = 0;
294	Quotient qerr;
295	p.resize(PP.size());
296
297	typename Vector <typename Rationals::Element>::iterator it;
298
299	Rationals Q;
300	for (it= p.begin(); it != p.end(); ++it, ++i) {
301	A.field().init(*it, PP[i],M[i]);
302	Q.get_den(t, *it);
303	Q.get_num(tt,*it);
304	err = M[i]/t;
305	size_t resi = err.bitsize() + tt.bitsize() -1;
306	size_t resu = t.bitsize() + tt.bitsize() -1;
307	if (resi > max_res) {max_res = resi; max_i=i;}
308	if (resu > max_resu) {max_resu = resu; max_iu =i;}
309	//size_t resu = t.bitsize() + tt.bitsize() -1;
310	//if (err.bitsize() > max_err) max_err = err.bitsize();
311	}
312
313	max_res=0;
314	for (it= p.begin()+1; it != p.end(); ++it) {
315	//A.field().init(*it, PP[i],M[i]);
316	Q.get_den(t, *it);
317	Q.get_den(tt, *(it-1));
318	Q.init(qerr,t,tt);
319	Q.get_num(tt, *it);
320	size_t resi = Q.bitsize(t,qerr) + tt.bitsize() -2;
321	if (resi > max_res) {max_res = resi; max_i=i;}
322	//if (err.bitsize() > max_err) max_err = err.bitsize();
323	}
324
325	commentator.stop ("done", NULL, "Iminpoly");
326
327	return p;
328
329	}
330
331	}
332
333	#endif //__LINBOX_charpoly_rational_H

Note: See TracBrowser for help on using the repository browser.

Download in other formats: