From 37c418c602485be9b023378a1b60ebdb4a592f2a Mon Sep 17 00:00:00 2001
From: Jens Vollinga <jensv@nikhef.nl>
Date: Wed, 20 Aug 2008 22:16:13 +0200
Subject: [PATCH 1/1] Fixed bug in the Q matrix calculation and the univariate
 factorization function.

---
 check/exam_factor.cpp |  9 +++++++-
 ginac/factor.cpp      | 51 +++++++++++++++++++++++++++++--------------
 2 files changed, 43 insertions(+), 17 deletions(-)

diff --git a/check/exam_factor.cpp b/check/exam_factor.cpp
index ed064588..bf1751a5 100644
--- a/check/exam_factor.cpp
+++ b/check/exam_factor.cpp
@@ -58,6 +58,9 @@ static unsigned exam_factor1()
 	e = ex("(1+x)^3", syms);
 	result += check_factor(e);
 
+	e = ex("(x+1)*(x+4)", syms);
+	result += check_factor(e);
+
 	e = ex("x^6-3*x^5+x^4-3*x^3-x^2-3*x+1", syms);
 	result += check_factor(e);
 
@@ -79,7 +82,11 @@ static unsigned exam_factor1()
 	e = ex("x^8-40*x^6+352*x^4-960*x^2+576", syms);
 	result += check_factor(e);
 
-	e = ex("x*(2+x^2)*(1+x+x^3+x^2+x^6+x^5+x^4)*(1+x^3)^2*(-1+x)", syms);
+	e = ex("x*(2+x^2)*(1+x+x^3+x^2+x^6+x^5+x^4)*(1+x)^2*(1-x+x^2)^2*(-1+x)", syms);
+
+	result += check_factor(e);
+
+	e = ex("(x+4+x^2-x^3+43*x^4)*(x+1-x^2-3*x^3+4*x^4)", syms);
 	result += check_factor(e);
 
 	return result;
diff --git a/ginac/factor.cpp b/ginac/factor.cpp
index 18e7c04c..5489b4ce 100644
--- a/ginac/factor.cpp
+++ b/ginac/factor.cpp
@@ -700,19 +700,32 @@ static void q_matrix(const UniPoly& a, Matrix& Q)
 {
 	unsigned int n = a.degree();
 	unsigned int q = cl_I_to_uint(a.R->modulus);
-	vector<cl_MI> r(n, a.R->zero());
-	r[0] = a.R->one();
-	Q.set_row(0, r);
-	unsigned int max = (n-1) * q;
-	for ( size_t m=1; m<=max; ++m ) {
-		cl_MI rn_1 = r.back();
-		for ( size_t i=n-1; i>0; --i ) {
-			r[i] = r[i-1] - rn_1 * a[i];
-		}
-		r[0] = -rn_1 * a[0];
-		if ( (m % q) == 0 ) {
-			Q.set_row(m/q, r);
+// fast and buggy
+//	vector<cl_MI> r(n, a.R->zero());
+//	r[0] = a.R->one();
+//	Q.set_row(0, r);
+//	unsigned int max = (n-1) * q;
+//	for ( size_t m=1; m<=max; ++m ) {
+//		cl_MI rn_1 = r.back();
+//		for ( size_t i=n-1; i>0; --i ) {
+//			r[i] = r[i-1] - rn_1 * a[i];
+//		}
+//		r[0] = -rn_1 * a[0];
+//		if ( (m % q) == 0 ) {
+//			Q.set_row(m/q, r);
+//		}
+//	}
+// slow and (hopefully) correct
+	for ( size_t i=0; i<n; ++i ) {
+		UniPoly qk(a.R);
+		qk.set(i*q, a.R->one());
+		UniPoly r(a.R);
+		rem(qk, a, r);
+		Vec rvec;
+		for ( size_t j=0; j<n; ++j ) {
+			rvec.push_back(r[j]);
 		}
+		Q.set_row(i, rvec);
 	}
 }
 
@@ -795,7 +808,12 @@ static void berlekamp(const UniPoly& a, UniPolyVec& upv)
 					*u = uo;
 				}
 				factors.push_back(g);
-				++size;
+				size = 0;
+				list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
+				while ( i != end ) {
+					if ( i->degree() ) ++size; 
+					++i;
+				}
 				if ( size == k ) {
 					list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
 					while ( i != end ) {
@@ -803,9 +821,9 @@ static void berlekamp(const UniPoly& a, UniPolyVec& upv)
 					}
 					return;
 				}
-				if ( u->degree() < nur.degree() ) {
-					break;
-				}
+//				if ( u->degree() < nur.degree() ) {
+//					break;
+//				}
 			}
 		}
 		if ( ++r == k ) {
@@ -1107,6 +1125,7 @@ static ex factor_univariate(const ex& poly, const ex& x)
 					mf.factors = newfactors2;
 					mf.poly = answer.op(1);
 					tocheck.push(mf);
+					break;
 				}
 			}
 			else {
-- 
2.44.0