From c434014237cb04e46af2e0a0c5a3e6d374121547 Mon Sep 17 00:00:00 2001
From: "Vladimir V. Kisil" <V.Kisilv@leeds.ac.uk>
Date: Thu, 23 Jul 2020 22:59:50 +0100
Subject: [PATCH 1/2] Fix of a bug in sqrfree_parfrac() related to Yun
 factorisation.

Yun's algorithm does not account for constant factors, thus
sqrfree_parfrac() drop such factors in partial fractions.

Signed-off-by: Vladimir V. Kisil <V.Kisilv@leeds.ac.uk>
---
 check/exam_paranoia.cpp | 16 ++++++++++++++++
 ginac/normal.cpp        | 20 +++++++++++++-------
 2 files changed, 29 insertions(+), 7 deletions(-)

diff --git a/check/exam_paranoia.cpp b/check/exam_paranoia.cpp
index d96cda59..52466b05 100644
--- a/check/exam_paranoia.cpp
+++ b/check/exam_paranoia.cpp
@@ -636,6 +636,21 @@ unsigned exam_paranoia25()
 	return 0;
 }
 
+// Bug in partial fraction expansion
+unsigned exam_paranoia26()
+{
+	symbol x("x");
+	ex ex1=pow(x,4)/(x-1)/4;
+	ex ex2=sqrfree_parfrac(ex1,x);
+	ex e = (ex1-ex2).normal();
+
+	if (! e.is_zero()) {
+		clog << "partial fraction expansion of " << ex1 << " produces error.\n";
+		return 1;
+	}
+	return 0;
+}
+
 unsigned exam_paranoia()
 {
 	unsigned result = 0;
@@ -668,6 +683,7 @@ unsigned exam_paranoia()
 	result += exam_paranoia23();  cout << '.' << flush;
 	result += exam_paranoia24();  cout << '.' << flush;
 	result += exam_paranoia25();  cout << '.' << flush;
+	result += exam_paranoia26();  cout << '.' << flush;
 	
 	return result;
 }
diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index c234a86f..784661dc 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1809,6 +1809,7 @@ static epvector sqrfree_yun(const ex &a, const symbol &x)
 	epvector results;
 	ex exponent = _ex0;
 	do {
+		//std::clog << "w=" << w << ", z=" << z << ", g=" << g << std::endl;
 		w = quo(w, g, x);
 		if (w.is_zero()) {
 			return results;
@@ -1936,14 +1937,19 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 	// Find numerator and denominator
 	ex nd = numer_denom(a);
 	ex numer = nd.op(0), denom = nd.op(1);
-//clog << "numer = " << numer << ", denom = " << denom << endl;
+//std::clog << "numer = " << numer << ", denom = " << denom << std::endl;
 
 	// Convert N(x)/D(x) -> Q(x) + R(x)/D(x), so degree(R) < degree(D)
 	ex red_poly = quo(numer, denom, x), red_numer = rem(numer, denom, x).expand();
-//clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << endl;
+//std::clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << std::endl;
 
 	// Factorize denominator and compute cofactors
 	epvector yun = sqrfree_yun(denom, x);
+	// Yun's algorithm does not account for constant factors.  (For univariate
+	// polynomials it works only in the monic case.)  We can correct this by
+	// inserting what has been lost back into the result.  For completeness
+	// we'll also have to recurse down that factor in the remaining variables.
+	ex lead_coeff  = quo(denom, mul(yun), x);
 	size_t yun_max_exponent = yun.empty() ? 0 : ex_to<numeric>(yun.back().coeff).to_int();
 	exvector factor, cofac;
 	for (size_t i=0; i<yun.size(); i++) {
@@ -1961,8 +1967,8 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 		}
 	}
 	size_t num_factors = factor.size();
-//clog << "factors  : " << exprseq(factor) << endl;
-//clog << "cofactors: " << exprseq(cofac) << endl;
+//std::clog << "factors  : " << exprseq(factor) << std::endl;
+//std::clog << "cofactors: " << exprseq(cofac) << std::endl;
 
 	// Construct coefficient matrix for decomposition
 	int max_denom_deg = denom.degree(x);
@@ -1973,8 +1979,8 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 			sys(i, j) = cofac[j].coeff(x, i);
 		rhs(i, 0) = red_numer.coeff(x, i);
 	}
-//clog << "coeffs: " << sys << endl;
-//clog << "rhs   : " << rhs << endl;
+//std::clog << "coeffs: " << sys << std::endl;
+//std::clog << "rhs   : " << rhs << std::endl;
 
 	// Solve resulting linear system
 	matrix vars(num_factors, 1);
@@ -1985,7 +1991,7 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 	// Sum up decomposed fractions
 	ex sum = 0;
 	for (size_t i=0; i<num_factors; i++)
-		sum += sol(i, 0) / factor[i];
+		sum += sol(i, 0) / factor[i] / lead_coeff;
 
 	return red_poly + sum;
 }
-- 
2.27.0

