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

Revision 4086, 9.4 KB checked in by bboyer, 3 years 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
43namespace 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.