From e5fedf6b3da7dc11b558628e1280cb252eafde04 Mon Sep 17 00:00:00 2001
From: Richard Kreckel <kreckel@ginac.de>
Date: Sun, 18 Sep 2022 21:36:31 +0200
Subject: [PATCH] Make sqrfree_parfrac() work in the general case and make it
 supported.

---
 check/exam_sqrfree.cpp  | 41 ++++++++++++++++++++++++++-
 doc/tutorial/ginac.texi | 21 ++++++++++++++
 ginac/normal.cpp        | 62 +++++++++++++++++++++++++++++------------
 ginsh/ginsh.1.in        |  3 ++
 ginsh/ginsh_parser.ypp  |  6 ++++
 5 files changed, 114 insertions(+), 19 deletions(-)

diff --git a/check/exam_sqrfree.cpp b/check/exam_sqrfree.cpp
index 79d776f6..3e94ed67 100644
--- a/check/exam_sqrfree.cpp
+++ b/check/exam_sqrfree.cpp
@@ -197,7 +197,46 @@ unsigned exam_sqrfree()
 	return result;
 }
 
+unsigned exam_sqrfree_parfrac()
+{
+	symbol x("x");
+	// (expression, number of terms after partial fraction decomposition)
+	vector<pair<ex, unsigned>> exams = {
+		{ex("(x - 1) / (x*(x^2 + 2))", lst{x}), 2},
+		{ex("(1 - x^10) / x", lst{x}), 2},
+		{ex("(2*x^3 + x + 3) / ((x^2 + 1)^2)", lst{x}), 2},
+		{ex("1 / (x * (x+1) * (x+2))", lst{x}), 3},
+		{ex("(x*x + 3*x - 1) / (x^2*(x^2 + 2)^2)", lst{x}), 4},
+		{ex("(1 - x^10) / (x + 2)", lst{x}), 11},
+		{ex("(1 - x + 3*x^2) / (x^3 * (2+x)^2)", lst{x}), 5},
+		{ex("(1 - x) / (x^4 * (x - 2)^3)", lst{x}), 6},
+		{ex("(1 - 2*x + x^9) / (x^5 * (1 - x + x^2)^5)", lst{x}), 12}
+	};
+	unsigned result = 0;
+
+	cout << "\n"
+	     << "examining square-free partial fraction decomposition" << flush;
+	for (auto e: exams) {
+		ex e1 = e.first;
+		ex e2 = sqrfree_parfrac(e1, x);
+		if (e2.nops() != e.second &&
+		    !normal(e1-e2).is_zero()) {
+			clog << "sqrfree_parfrac(" << e1 << ", " << x << ") erroneously returned "
+			     << e2 << endl;
+			++result;
+		}
+		cout << '.' << flush;
+	}
+
+	return result;
+}
+
 int main(int argc, char** argv)
 {
-	return exam_sqrfree();
+	unsigned result = 0;
+
+	result += exam_sqrfree();
+	result += exam_sqrfree_parfrac();
+
+	return result;
 }
diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index c840fb69..40368aef 100644
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -5379,6 +5379,27 @@ some care with subsequent processing of the result:
 Note also, how factors with the same exponents are not fully factorized
 with this method.
 
+@subsection Square-free partial fraction decomposition
+@cindex square-free partial fraction decomposition
+@cindex partial fraction decomposition
+@cindex @code{sqrfree_parfrac()}
+
+GiNaC also supports square-free partial fraction decomposition of
+rational functions:
+@example
+ex sqrfree_parfrac(const ex & a, const symbol & x);
+@end example
+It is called square-free because it assumes a square-free
+factorization of the input's denominator:
+@example
+    ...
+    symbol x("x");
+
+    ex rat = (x-4)/(pow(x,2)*(x+2));
+    cout << sqrfree_parfrac(rat, x) << endl;
+     // -> -2*x^(-2)+3/2*x^(-1)-3/2*(2+x)^(-1)
+@end example
+
 @subsection Polynomial factorization
 @cindex factorization
 @cindex polynomial factorization
diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 383115f0..cad9d455 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1958,10 +1958,12 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 	// Factorize denominator and compute cofactors
 	epvector yun = sqrfree_yun(denom, x);
 	exvector factor, cofac;
+	int dim = 0;
 	for (size_t i=0; i<yun.size(); i++) {
 		numeric i_exponent = ex_to<numeric>(yun[i].coeff);
 		for (size_t j=0; j<i_exponent; j++) {
 			factor.push_back(pow(yun[i].rest, j+1));
+			dim += degree(yun[i].rest, x);
 			ex prod = _ex1;
 			for (size_t k=0; k<yun.size(); k++) {
 				if (yun[k].coeff == i_exponent)
@@ -1972,34 +1974,58 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 			cofac.push_back(prod.expand());
 		}
 	}
-	size_t num_factors = factor.size();
 //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);
-	matrix sys(max_denom_deg + 1, num_factors);
-	matrix rhs(max_denom_deg + 1, 1);
-	for (int i=0; i<=max_denom_deg; i++) {
-		for (size_t j=0; j<num_factors; j++)
-			sys(i, j) = cofac[j].coeff(x, i);
-		rhs(i, 0) = red_numer.coeff(x, i);
+	// Construct linear system for decomposition
+	matrix sys(dim, dim);
+	matrix rhs(dim, 1);
+	matrix vars(dim, 1);
+	for (size_t i=0, n=0, f=0; i<yun.size(); i++) {
+		int i_expo = to_int(ex_to<numeric>(yun[i].coeff));
+		for (size_t j=0; j<i_expo; j++) {
+			for (size_t k=0; k<degree(yun[i].rest, x); k++) {
+				GINAC_ASSERT(n < dim  &&  f < factor.size());
+
+				// column n of coefficient matrix
+				for (size_t r=0; r+k<dim; r++) {
+					sys(r+k, n) = cofac[f].coeff(x, r);
+				}
+
+				// element n of right hand side vector
+				rhs(n, 0) = red_numer.coeff(x, n);
+
+				// element n of free variables vector
+				vars(n, 0) = symbol();
+
+				n++;
+			}
+			f++;
+		}
 	}
 //std::clog << "coeffs: " << sys << std::endl;
 //std::clog << "rhs   : " << rhs << std::endl;
 
-	// Solve resulting linear system
-	matrix vars(num_factors, 1);
-	for (size_t i=0; i<num_factors; i++)
-		vars(i, 0) = symbol();
+	// Solve resulting linear system and sum up decomposed fractions
 	matrix sol = sys.solve(vars, rhs);
+//std::clog << "sol   : " << sol << std::endl;
+	ex sum = red_poly;
+	for (size_t i=0, n=0, f=0; i<yun.size(); i++) {
+		int i_expo = to_int(ex_to<numeric>(yun[i].coeff));
+		for (size_t j=0; j<i_expo; j++) {
+			ex frac_numer = 0;
+			for (size_t k=0; k<degree(yun[i].rest, x); k++) {
+				GINAC_ASSERT(n < dim  &&  f < factor.size());
+				frac_numer += sol(n, 0) * pow(x, k);
+				n++;
+			}
+			sum += frac_numer / factor[f];
 
-	// Sum up decomposed fractions
-	ex sum = 0;
-	for (size_t i=0; i<num_factors; i++)
-		sum += sol(i, 0) / factor[i];
+			f++;
+		}
+	}
 
-	return red_poly + sum;
+	return sum;
 }
 
 
diff --git a/ginsh/ginsh.1.in b/ginsh/ginsh.1.in
index 61741683..37aedd38 100644
--- a/ginsh/ginsh.1.in
+++ b/ginsh/ginsh.1.in
@@ -377,6 +377,9 @@ detail here. Please refer to the GiNaC documentation.
 .BI sqrfree( "expression [" ", " symbol-list] )
 \- square-free factorization of a polynomial
 .br
+.BI sqrfree_parfrac( expression ", " symbol )
+\- square-free partial fraction decomposition of rational function
+.br
 .BI sqrt( expression )
 \- square root
 .br
diff --git a/ginsh/ginsh_parser.ypp b/ginsh/ginsh_parser.ypp
index 88d1ce78..8fb9f13d 100644
--- a/ginsh/ginsh_parser.ypp
+++ b/ginsh/ginsh_parser.ypp
@@ -547,6 +547,11 @@ static ex f_sqrfree2(const exprseq &e)
 	return sqrfree(e[0], ex_to<lst>(e[1]));
 }
 
+static ex f_sqrfree_parfrac(const exprseq &e)
+{
+	return sqrfree_parfrac(e[0], ex_to<symbol>(e[1]));
+}
+
 static ex f_subs3(const exprseq &e)
 {
 	CHECK_ARG(1, lst, subs);
@@ -747,6 +752,7 @@ static const fcn_init builtin_fcns[] = {
 	{"sprem", f_sprem, 3},
 	{"sqrfree", f_sqrfree1, 1},
 	{"sqrfree", f_sqrfree2, 2},
+	{"sqrfree_parfrac", f_sqrfree_parfrac, 2},
 	{"sqrt", f_sqrt, 1},
 	{"subs", f_subs2, 2},
 	{"subs", f_subs3, 3},
-- 
2.49.0