From a12c1d88131c5a301d35767f0c4c947064703418 Mon Sep 17 00:00:00 2001
From: Jens Vollinga <jensv@balin.nikhef.nl>
Date: Thu, 18 Dec 2008 11:10:50 +0100
Subject: [PATCH] Fixed bug in multivariate factorization. Fractional numerical
 content was not extracted and led to a not so finite running behavior. Added
 new checks for that bug.

Added static keywords to hide debugging output operators.
---
 check/exam_factor.cpp | 22 ++++++++++++++++++++--
 ginac/factor.cpp      | 15 +++++----------
 2 files changed, 25 insertions(+), 12 deletions(-)

diff --git a/check/exam_factor.cpp b/check/exam_factor.cpp
index 8851b7cc..71a63f24 100644
--- a/check/exam_factor.cpp
+++ b/check/exam_factor.cpp
@@ -41,7 +41,7 @@ static unsigned check_factor(const ex& e)
 static unsigned exam_factor1()
 {
 	unsigned result = 0;
-	ex e, d;
+	ex e;
 	symbol x("x");
 	lst syms;
 	syms.append(x);
@@ -110,7 +110,7 @@ static unsigned exam_factor1()
 static unsigned exam_factor2()
 {
 	unsigned result = 0;
-	ex e, d;
+	ex e;
 	symbol x("x"), y("y"), z("z");
 	lst syms;
 	syms = x, y, z;
@@ -166,6 +166,23 @@ static unsigned exam_factor2()
 	return result;
 }
 
+static unsigned exam_factor3()
+{
+	unsigned result = 0;
+	ex e;
+	symbol k("k"), n("n");
+	lst syms;
+	syms = k, n;
+	
+	e = ex("1/2*(-3+3*k-n)*(-2+3*k-n)*(-1+3*k-n)", syms);
+	result += check_factor(e);
+
+	e = ex("1/4*(2*k-n)*(-1+2*k-n)", syms);
+	result += check_factor(e);
+
+	return result;
+}
+
 unsigned exam_factor()
 {
 	unsigned result = 0;
@@ -174,6 +191,7 @@ unsigned exam_factor()
 
 	result += exam_factor1(); cout << '.' << flush;
 	result += exam_factor2(); cout << '.' << flush;
+	result += exam_factor3(); cout << '.' << flush;
 
 	return result;
 }
diff --git a/ginac/factor.cpp b/ginac/factor.cpp
index b7129135..da870d3d 100644
--- a/ginac/factor.cpp
+++ b/ginac/factor.cpp
@@ -92,7 +92,7 @@ ostream& operator<<(ostream& o, const vector<int>& v)
 	}
 	return o;
 }
-ostream& operator<<(ostream& o, const vector<cl_I>& v)
+static ostream& operator<<(ostream& o, const vector<cl_I>& v)
 {
 	vector<cl_I>::const_iterator i = v.begin(), end = v.end();
 	while ( i != end ) {
@@ -101,7 +101,7 @@ ostream& operator<<(ostream& o, const vector<cl_I>& v)
 	}
 	return o;
 }
-ostream& operator<<(ostream& o, const vector<cl_MI>& v)
+static ostream& operator<<(ostream& o, const vector<cl_MI>& v)
 {
 	vector<cl_MI>::const_iterator i = v.begin(), end = v.end();
 	while ( i != end ) {
@@ -2237,13 +2237,8 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 	const ex& x = *syms.begin();
 
 	// make polynomial primitive
-	ex p = poly.collect(x);
-	ex cont = p.lcoeff(x);
-	for ( int i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
-		cont = gcd(cont, p.coeff(x,i));
-		if ( cont == 1 ) break;
-	}
-	ex pp = expand(normal(p / cont));
+	ex unit, cont, pp;
+	poly.unitcontprim(x, unit, cont, pp);
 	if ( !is_a<numeric>(cont) ) {
 		return factor_sqrfree(cont) * factor_sqrfree(pp);
 	}
@@ -2427,7 +2422,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 		// try Hensel lifting
 		ex res = hensel_multivar(pp, x, epv, prime, l, modfactors, C);
 		if ( res != lst() ) {
-			ex result = cont;
+			ex result = cont * unit;
 			for ( size_t i=0; i<res.nops(); ++i ) {
 				result *= res.op(i).content(x) * res.op(i).unit(x);
 				result *= res.op(i).primpart(x);
-- 
2.44.0