Fix of a bug in sqrfree_parfrac() related to Yun factorisation.
authorVladimir V. Kisil <V.Kisilv@leeds.ac.uk>
Mon, 27 Jul 2020 18:58:27 +0000 (20:58 +0200)
committerRichard Kreckel <kreckel@ginac.de>
Mon, 27 Jul 2020 18:58:27 +0000 (20:58 +0200)
Yun's algorithm does not account for constant factors, thus
sqrfree_parfrac() dropped such factors in partial fractions.

check/exam_paranoia.cpp
ginac/normal.cpp

index d96cda59c035b0a37443b9719b370bd56cd75f0f..52466b059b055fe2884636de1aef5a063cfa72c8 100644 (file)
@@ -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;
 }
index 8e44bee24db29555bb4ad43794b560bcabc17ad3..60c09d9cc4669c0d4fa859b5c69f09c163fbf64c 100644 (file)
@@ -1936,14 +1936,18 @@ 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.
+       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 +1965,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 +1977,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 +1989,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;
 }