modify sequence container for wiedemann minpoly to ensure that generatedsequence is never the zero sequence, add proper test for charpoly

Author: pascalgiorgi

linbox/algorithms/bbcharpoly.h

@@ -285,8 +285,8 @@ namespace LinBox
 
 			/* Computation of the minimal polynomial */
 			Polynomial minPoly(F);
-			minpoly (minPoly, A, M);
-//             PD.write(std::cerr<<"Minpoly = ",minPoly) << std::endl;
+			minpoly (minPoly, A, M);            
+            //PD.write(std::cerr<<"Minpoly = ",minPoly) << std::endl;
 
 			if (minPoly.size() == n+1){
 				commentator().stop ("done", NULL, "MbbCharpoly");

linbox/algorithms/blackbox-container-base.h

@@ -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

@@ -87,13 +87,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,7 @@ namespace LinBox
 
 			commentator().stop ("done", NULL, "masseyd");
 			//		commentator().stop ("Done", "Done", "LinBox::MasseyDomain::massey");
-
+            std::cerr<<"SEQ: "<<S<<std::endl;
 			return L;
 		}
 
@@ -398,7 +398,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

@@ -600,6 +600,37 @@ bool testZeroMatrixCharPoly() {
 
     return success;
 }
+bool testZeroSymmetricMatrixCharPoly() {
+    bool success;
+	using Ring = Givaro::Modular<double>;
+	using Matrix = SparseMatrix<Ring>;
+    Ring R(3);
+
+    Matrix A(R, 1, 1);
+    A.setEntry(0, 0, R.zero);
+
+    PolynomialRing<Ring>::Element c_A, Ex;
+
+
+    charpoly(c_A, A, RingCategories::ModularTag(), Method::Blackbox(Shape::Symmetric));
+
+    PolynomialRing<Ring> PZ(R,'X'); PZ.assign(Ex, Givaro::Degree(1), R.one);
+
+    success = PZ.areEqual(c_A, Ex);
+
+    if (!success) {
+        if (writing) {
+            std::clog<<"**** ERROR **** Fail ZSMCP " <<std::endl;
+
+            PZ.write(std::clog << "Ex: ", Ex) << std::endl;
+            PZ.write(std::clog << "cA: ", c_A) << std::endl;
+        }
+        return false;
+    } else
+        if (writing) std::cout << "ZSMCP: PASSED" << std::endl;
+
+    return success;
+}
 
 bool testFourFourMatrix() {
     bool success;
@@ -676,6 +707,7 @@ int main (int argc, char **argv)
     pass &= testZeroDimensionalCharPoly ();
     pass &= testZeroDimensionalMinPoly ();
     pass &= testZeroMatrixCharPoly();
+    pass &= testZeroSymmetricMatrixCharPoly();  
     pass &= testBigScalarCharPoly ();
     pass &= testLocalSmith ();
     pass &= testInconsistent<DenseMatrix<ZRingInts>> (Method::DenseElimination());