Make sqrfree_parfrac() work in the general case and make it supported.

diff --git a/check/exam_sqrfree.cpp b/check/exam_sqrfree.cpp
index 79d776f689e54c4edd595c262718f41630eccb77..3e94ed67e68afb4efc53c605f0381d65f6b25be2 100644 (file)
--- 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 c840fb69668c30f68c92e6d90d9306fe2c74370f..40368aef43dd22fc6665d9cd18b3694f5650193e 100644 (file)
--- 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 383115f0ce7fe911100aed8b57d6389dcefb04f9..cad9d455bcf9ace6e4db1fbe57d5a37346498358 100644 (file)
--- 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 617416836b1e327ee5872170f9ad63deaad9c4a9..37aedd38ef8298569df194206cfa5717f791a4ef 100644 (file)
--- 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 88d1ce787b12686597937c46c4a4fd5073e4b28d..8fb9f13d4f2f6c5bd5107c343be917dee76d44f3 100644 (file)
--- 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},