Skip to content

Commit 4c70455

Browse files
ClementPernetClementPernet
authored
Merge pull request #214 from linbox-team/ciadebug
Ciadebug
2 parents 411b132 + 75f88ac commit 4c70455

File tree

5 files changed

+53
-29
lines changed

5 files changed

+53
-29
lines changed

‎linbox/algorithms/bbcharpoly.h‎

Lines changed: 13 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -182,7 +182,7 @@ namespace LinBox
182182
/* Factorization over the integers */
183183
std::vector<IntPoly> intFactors;
184184
std::vector<uint64_t> exp;
185-
IPD.factor (intFactors, exp, IntPoly(intMinPoly.getRep().begin(),intMinPoly.getRep().end()));
185+
IPD.factor (intFactors, exp, intMinPoly);
186186
size_t factnum = intFactors.size();
187187

188188
/* Choose a modular prime field */
@@ -195,14 +195,17 @@ namespace LinBox
195195

196196
for (size_t i = 0; i < intFactors.size(); ++i) {
197197
uint64_t deg = (intFactors[i].size()-1);
198+
199+
198200
FactorMult<FieldPoly,IntPoly>* FFM=NULL;
199201
if (exp[i] > 1) {
200202
IntPoly *tmp = new IntPoly(intFactors[i]);
201203
FM* depend = NULL;
202204
for (size_t j = 1; j <= exp[i]; ++j){
203205
IntPoly * tmp2 = new IntPoly(*tmp);
204-
FieldPoly *tmp2p = new FieldPoly(tmp2->size());
205-
typename IntPoly::template rebind<Field>() (*tmp2p, *tmp2, F);
206+
// FieldPoly *tmp2p = new FieldPoly(tmp2->size());
207+
// typename IntPoly::template rebind<Field>() (*tmp2p, *tmp2, F);
208+
FieldPoly *tmp2p = new FieldPoly(*tmp2, F);
206209

207210
FFM = new FM (tmp2p, tmp2, 0, depend);
208211
factCharPoly.insert (std::pair<size_t, FM*> (deg, FFM));
@@ -219,8 +222,10 @@ namespace LinBox
219222
leadingBlocks.insert (std::pair<FM*,bool>(FFM,false));
220223
}
221224
else {
222-
FieldPoly* fp=new FieldPoly(intFactors[i].size());
223-
typename IntPoly::template rebind<Field>() (*fp, (intFactors[i]), F);
225+
// FieldPoly* fp=new FieldPoly(intFactors[i].size());
226+
// typename IntPoly::template rebind<Field>() (*fp, (intFactors[i]), F);
227+
FieldPoly* fp=new FieldPoly(intFactors[i], F);
228+
224229
FFM = new FM (fp,&intFactors[i],1,NULL);
225230
factCharPoly.insert (std::pair<size_t, FM* > (intFactors[i].size()-1, FFM));
226231
leadingBlocks.insert (std::pair<FM*,bool>(FFM,false));
@@ -237,6 +242,8 @@ namespace LinBox
237242
IntPoly tmpP;
238243
intRing.assign(intCharPoly[0], intRing.one);
239244
for (FactPolyIterator it_f = factCharPoly.begin(); it_f != factCharPoly.end(); ++it_f){
245+
// it_f->second->write(std::clog) << std::endl;
246+
240247
IPD.pow (tmpP, *it_f->second->intP, (long) it_f->second->multiplicity);
241248
IPD.mulin (intCharPoly, tmpP);
242249
delete it_f->second->intP;
@@ -245,7 +252,7 @@ namespace LinBox
245252
}
246253
commentator().stop ("done", NULL, "IbbCharpoly");
247254

248-
return P = Polynomial(A.field(), typename Polynomial::Rep(intCharPoly.begin(),intCharPoly.end()));
255+
return P = Polynomial(A.field(), intCharPoly);
249256
}
250257

251258
/** Algorithm computing the characteristic polynomial

‎linbox/algorithms/cia.h‎

Lines changed: 16 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -51,21 +51,27 @@ namespace LinBox
5151
{
5252
commentator().start ("Integer Givaro::Dense Charpoly ", "CIA");
5353

54-
typename Blackbox::Field intRing = A.field();
54+
using Ring = typename Blackbox::Field;
55+
Ring ZZ = A.field();
5556
typedef Givaro::Modular<double> Field;
5657
typedef typename Blackbox::template rebind<Field>::other FBlackbox;
57-
typedef Givaro::Poly1FactorDom<typename Blackbox::Field, Givaro::Dense> IntPolyDom;
58-
typedef Givaro::Poly1FactorDom<Field, Givaro::Dense> FieldPolyDom;
58+
typedef PolynomialRing<Ring> IntPolyDom;
59+
typedef PolynomialRing<Field> FieldPolyDom;
5960
typedef typename IntPolyDom::Element IntPoly;
6061
typedef typename FieldPolyDom::Element FieldPoly;
6162

62-
IntPolyDom IPD(intRing);
63+
IntPolyDom IPD(ZZ);
6364

64-
FieldPoly fieldCharPoly(A.coldim());
6565
/* Computation of the integer minimal polynomial */
66-
IntPoly intMinPoly;
66+
IntPoly intMinPoly(ZZ);
6767
minpoly (intMinPoly, A, RingCategories::IntegerTag(), M);
6868

69+
if (intMinPoly.size() == A.coldim()+1){
70+
commentator().stop ("done", NULL, "CIA");
71+
return P = intMinPoly;
72+
}
73+
// IPD.write(std::cerr<<"Minpoly = ", intMinPoly) << std::endl;
74+
6975
/* Factorization over the integers */
7076
std::vector<IntPoly> intFactors;
7177
std::vector<uint64_t> mult;
@@ -77,6 +83,7 @@ namespace LinBox
7783
++primeg;
7884
Field F(*primeg);
7985
FBlackbox fbb(F, A.rowdim(), A.coldim());
86+
FieldPoly fieldCharPoly(F);
8087
MatrixHom::map(fbb, A);
8188
charpoly (fieldCharPoly, fbb, M);
8289
/* Determination of the multiplicities */
@@ -88,7 +95,7 @@ namespace LinBox
8895
fieldFactors[i].resize(d);
8996
for (size_t j = 0; j < d; ++j)
9097
//F.init ((fieldFactors[i])[j], (*intFactors[i])[j]);
91-
F.init ((fieldFactors[i])[j], intRing.convert(tmp_convert,(intFactors[i])[j]));// PG 2005-08-04
98+
F.init ((fieldFactors[i])[j], ZZ.convert(tmp_convert,(intFactors[i])[j]));// PG 2005-08-04
9299
}
93100

94101
FieldPoly currPol = fieldCharPoly;
@@ -107,8 +114,8 @@ namespace LinBox
107114
multip[i] = m-1;
108115
}
109116

110-
IntPoly intCharPoly (A.coldim());
111-
intRing.assign (intCharPoly[0], intRing.one);
117+
IntPoly intCharPoly (ZZ);
118+
IPD.assign (intCharPoly, IPD.one);
112119
for (size_t i = 0; i < nf; ++i){
113120
IPD.pow( P, intFactors[i], multip[i] );
114121
IPD.mulin( intCharPoly, P );

‎linbox/algorithms/wiedemann.h‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -130,6 +130,15 @@ namespace LinBox
130130
#endif // INCLUDE_TIMING
131131
}
132132

133+
// std::cerr << "P: " << P << std::endl;
134+
// std::cerr << "WD deg: " << deg << std::endl;
135+
if (!deg) {
136+
// zero sequence, matrix minpoly is X
137+
P.resize(2);
138+
A.field().assign(P[0],A.field().zero);
139+
A.field().assign(P[1],A.field().one);
140+
}
141+
133142
commentator().stop ("done", NULL, "minpoly");
134143

135144
return P;

‎linbox/solutions/charpoly.h‎

Lines changed: 2 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -203,8 +203,7 @@ namespace LinBox
203203

204204
}
205205

206-
#if 0 // CIA is buggy (try with matrix(ZZ,4,[1..16]) )
207-
//#if defined(__LINBOX_HAVE_NTL)
206+
#if defined(__LINBOX_HAVE_NTL)
208207

209208
#include "linbox/algorithms/cia.h"
210209
namespace LinBox
@@ -227,9 +226,7 @@ namespace LinBox
227226
{
228227
if (A.coldim() != A.rowdim())
229228
throw LinboxError("LinBox ERROR: matrix must be square for characteristic polynomial computation\n");
230-
typename Givaro::Poly1Dom<typename Blackbox::Field, Givaro::Dense>::Element Pg;
231-
std::cerr<<"CIA"<<std::endl;
232-
return P = cia (Pg, A, M);
229+
return cia (P, A, M);
233230
}
234231

235232
/** Compute the characteristic polynomial over {\bf Z}

‎tests/test-regression.C‎

Lines changed: 13 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -35,14 +35,15 @@
3535
#include "linbox/vector/blas-vector.h"
3636
#include "linbox/solutions/solve.h"
3737
#include "linbox/solutions/charpoly.h"
38-
#include "linbox/algorithms/smith-form-sparseelim-poweroftwo.h"
39-
#include "test-smith-form.h"
4038

4139
using namespace LinBox;
4240
typedef Givaro::ZRing<Givaro::Integer> ZRingInts;
4341

4442
bool writing=false;
4543

44+
#include "linbox/algorithms/smith-form-sparseelim-poweroftwo.h"
45+
#include "test-smith-form.h"
46+
4647
bool testSolveSparse(){
4748

4849
typedef Givaro::Modular<unsigned int> Field;
@@ -587,11 +588,12 @@ bool testZeroMatrixCharPoly() {
587588
success = PZ.areEqual(c_A, Ex);
588589

589590
if (!success) {
590-
if (writing) std::clog<<"**** ERROR **** Fail ZMCP " <<std::endl;
591-
592-
PZ.write(std::clog << "Ex: ", Ex) << std::endl;
593-
PZ.write(std::clog << "cA: ", c_A) << std::endl;
591+
if (writing) {
592+
std::clog<<"**** ERROR **** Fail ZMCP " <<std::endl;
594593

594+
PZ.write(std::clog << "Ex: ", Ex) << std::endl;
595+
PZ.write(std::clog << "cA: ", c_A) << std::endl;
596+
}
595597
return false;
596598
} else
597599
if (writing) std::cout << "ZMCP: PASSED" << std::endl;
@@ -621,10 +623,12 @@ bool testFourFourMatrix() {
621623
success = PZ.areEqual(c_A, Res);
622624

623625
if (!success) {
624-
if (writing) std::clog<<"**** ERROR **** Fail tFFM " <<std::endl;
626+
if (writing) {
627+
std::clog<<"**** ERROR **** Fail tFFM " <<std::endl;
625628

626-
PZ.write(std::clog << "Ex: ", Res) << std::endl;
627-
PZ.write(std::clog << "cA: ", c_A) << std::endl;
629+
PZ.write(std::clog << "Ex: ", Res) << std::endl;
630+
PZ.write(std::clog << "cA: ", c_A) << std::endl;
631+
}
628632

629633
return false;
630634
} else

0 commit comments

Comments
 (0)