Merge pull request linbox-team#288 from linbox-team/zeromatrixminpoly

problem with zero sequence for minpoly/charpoly

Author: jgdumas

linbox/algorithms/bbcharpoly.h

@@ -204,8 +204,8 @@ namespace LinBox
 					FM* depend = NULL;
 					for (size_t j = 1; j <= exp[i]; ++j){
 						IntPoly * tmp2 = new IntPoly(*tmp);
-// 						FieldPoly *tmp2p = new FieldPoly(tmp2->size());
-// 						typename IntPoly::template rebind<Field>() (*tmp2p, *tmp2, F);
+//						FieldPoly *tmp2p = new FieldPoly(tmp2->size());
+//						typename IntPoly::template rebind<Field>() (*tmp2p, *tmp2, F);
 						FieldPoly *tmp2p = new FieldPoly(*tmp2, F);
 
 						FFM = new FM (tmp2p, tmp2, 0, depend);
@@ -223,8 +223,8 @@ namespace LinBox
 					leadingBlocks.insert (std::pair<FM*,bool>(FFM,false));
 				}
 				else {
-// 					FieldPoly* fp=new FieldPoly(intFactors[i].size());
-// 					typename IntPoly::template rebind<Field>() (*fp, (intFactors[i]), F);
+//					FieldPoly* fp=new FieldPoly(intFactors[i].size());
+//					typename IntPoly::template rebind<Field>() (*fp, (intFactors[i]), F);
                     IntPoly * ip=new IntPoly(intFactors[i]);
 					FieldPoly* fp=new FieldPoly(intFactors[i], F);
 
@@ -286,14 +286,14 @@ namespace LinBox
 			/* Computation of the minimal polynomial */
 			Polynomial minPoly(F);
 			minpoly (minPoly, A, M);
-//             PD.write(std::cerr<<"Minpoly = ",minPoly) << std::endl;
+            //PD.write(std::cerr<<"Minpoly = ",minPoly) << std::endl;
 
 			if (minPoly.size() == n+1){
 				commentator().stop ("done", NULL, "MbbCharpoly");
 				return P = minPoly;
 			}
 
-            
+
 
 			Polynomial charPoly (F);
 

linbox/algorithms/blackbox-container-base.h

@@ -32,7 +32,7 @@
 // - wait   : waits for the end of the current computation
 // ================================================================
 
-
+#include "linbox/solutions/constants.h"
 #include "linbox/vector/vector-domain.h"
 
 namespace LinBox
@@ -157,11 +157,20 @@ namespace LinBox
 		{
 			casenumber = 1;
 			u.resize (_BB->coldim ());
-			for (long i = (long)u.size (); i--;)
-				g.random (u[(size_t)i]);
-			v.resize (_BB->rowdim ());
-			return _VD.dot (_value, u, u);
-		}
+            v.resize (_BB->rowdim ());
+
+            size_t trials=0;
+            do {
+                for (long i = (long)this->u.size (); i--;)
+                    g.random (this->u[(size_t)i]);
+                this->_VD.dot (this->_value, this->u, this->u);
+            } while(_field->isZero(this->_value) && ++trials<= LINBOX_DEFAULT_TRIALS_BEFORE_FAILURE);
+
+            if (trials >= LINBOX_DEFAULT_TRIALS_BEFORE_FAILURE)
+                std::cerr<<"ERROR in "<<__FILE__<<" at line "<<__LINE__<<" -> projection always auto-orthogonal after "<<LINBOX_DEFAULT_TRIALS_BEFORE_FAILURE<<" attempts\n";;
+
+            return _value;
+        }
 
 		/// User Left vectors, Zero Right vector
 		template<class Vector>

linbox/algorithms/blackbox-container.h

@@ -33,6 +33,7 @@
 #include "linbox/randiter/archetype.h"
 #include "linbox/algorithms/blackbox-container-base.h"
 #include "linbox/util/timer.h"
+#include "linbox/solutions/constants.h"
 
 namespace LinBox
 {
@@ -87,13 +88,23 @@ namespace LinBox
 		{
 			this->casenumber = 1;
 			this->u.resize (this->_BB->coldim ());
-			for (long i = (long)this->u.size (); i--;)
-				g.random (this->u[(size_t)i]);
-			this->w.resize (this->_BB->coldim ());
-			for (long i = (long)this->w.size (); i--;)
-				g.random (this->w[(size_t)i]);
-			this->v.resize (this->_BB->rowdim ());
-			this->_VD.dot (this->_value, this->u, this->w);
+            this->w.resize (this->_BB->coldim ());
+            this->v.resize (this->_BB->rowdim ());
+
+            size_t trials=0;
+            do {
+                for (long i = (long)this->u.size (); i--;)
+                    g.random (this->u[(size_t)i]);
+                
+                for (long i = (long)this->w.size (); i--;)
+                    g.random (this->w[(size_t)i]);
+			
+                this->_VD.dot (this->_value, this->u, this->w);
+            } while(F.isZero(this->_value) && ++trials<= LINBOX_DEFAULT_TRIALS_BEFORE_FAILURE);
+
+            if (trials >= LINBOX_DEFAULT_TRIALS_BEFORE_FAILURE)
+                std::cerr<<"ERROR in "<<__FILE__<<" at line "<<__LINE__<<" -> projection always orthogonal after "<<LINBOX_DEFAULT_TRIALS_BEFORE_FAILURE<<" attempts\n";;
+                
 #ifdef INCLUDE_TIMING
 			_applyTime = _dotTime = 0.0;
 #endif

linbox/algorithms/massey-domain.h

@@ -348,7 +348,6 @@ namespace LinBox
 
 			commentator().stop ("done", NULL, "masseyd");
 			//		commentator().stop ("Done", "Done", "LinBox::MasseyDomain::massey");
-
 			return L;
 		}
 
@@ -398,7 +397,8 @@ namespace LinBox
 		{
 			long dp = massey (phi, full_poly);
 			rank = (size_t) (v_degree(phi) - v_val (phi));
-
+            //std::cerr<<"dp="<<dp<<" rank="<<rank<<" -> "<<phi.size()<<" : "<<phi<<std::endl;
+            
             if (dp==0){// empty sequence, its minpoly is 1		
                  phi.resize(1);		
                  field().assign(phi[0],field().one);		

linbox/algorithms/wiedemann.h

@@ -129,7 +129,7 @@ namespace LinBox
 			<< "Time required for LSR fix:      " << WD.fixTime () << std::endl;
 #endif // INCLUDE_TIMING
 		}
-               
+
 		commentator().stop ("done", NULL, "minpoly");
 
 		return P;

tests/test-regression.C

@@ -570,7 +570,8 @@ bool testDixonSmallFat() {
     return success = ZZ.areEqual(r[0], b[0]);
 }
 
-bool testZeroMatrixCharPoly() {
+template<typename CMethod=Method::Blackbox>
+bool testZeroMatrixCharPoly(const CMethod& cMeth=CMethod()) {
     bool success;
 	using Ring = Givaro::Modular<double>;
 	using Matrix = SparseMatrix<Ring>;
@@ -581,7 +582,7 @@ bool testZeroMatrixCharPoly() {
 
     PolynomialRing<Ring>::Element c_A, Ex;
 
-    charpoly(c_A, A);
+    charpoly(c_A, A, RingCategories::ModularTag(), cMeth );
 
     PolynomialRing<Ring> PZ(R,'X'); PZ.assign(Ex, Givaro::Degree(1), R.one);
 
@@ -596,7 +597,7 @@ bool testZeroMatrixCharPoly() {
         }
         return false;
     } else
-        if (writing) std::cout << "ZMCP: PASSED" << std::endl;
+        if (writing) std::cout << "ZMCP" << cMeth.name() << cMeth.shapeFlags.flags << " : PASSED" << std::endl;
 
     return success;
 }
@@ -675,7 +676,10 @@ int main (int argc, char **argv)
     pass &= testSparseDiagDet(46);
     pass &= testZeroDimensionalCharPoly ();
     pass &= testZeroDimensionalMinPoly ();
-    pass &= testZeroMatrixCharPoly();
+    pass &= testZeroMatrixCharPoly<>();
+    pass &= testZeroMatrixCharPoly(Method::Blackbox(Shape::Symmetric));
+    pass &= testZeroMatrixCharPoly(Method::DenseElimination());
+//     pass &= testZeroMatrixCharPoly(Method::SparseElimination());
     pass &= testBigScalarCharPoly ();
     pass &= testLocalSmith ();
     pass &= testInconsistent<DenseMatrix<ZRingInts>> (Method::DenseElimination());