X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=5fda185981b13fdd82e800ff842db8f3ec2e1ce5;hp=462a8ad5192a8805dc93d20be4a363d49438fccd;hb=9f3e78a13c72ca73ade4d0829186edcdc1c1b2a0;hpb=9413cd14faaf2980de3884915e22a5beda068ecc

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 462a8ad5..5fda1859 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -6,7 +6,7 @@
  *  computation, square-free factorization and rational function normalization. */
 
 /*
- *  GiNaC Copyright (C) 1999-2016 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2019 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -146,7 +146,7 @@ struct sym_desc {
 	/** Maximum number of terms of leading coefficient of symbol in both polynomials */
 	size_t max_lcnops;
 
-	/** Commparison operator for sorting */
+	/** Comparison operator for sorting */
 	bool operator<(const sym_desc &x) const
 	{
 		if (max_deg == x.max_deg)
@@ -212,10 +212,10 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
 
 #if 0
 	std::clog << "Symbols:\n";
-	it = v.begin(); itend = v.end();
+	auto it = v.begin(), itend = v.end();
 	while (it != itend) {
-		std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << endl;
-		std::clog << "  lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << endl;
+		std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << std::endl;
+		std::clog << "  lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << std::endl;
 		++it;
 	}
 #endif
@@ -1470,7 +1470,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
 	}
 
 	// Some trivial cases
-	ex aex = a.expand(), bex = b.expand();
+	ex aex = a.expand();
 	if (aex.is_zero()) {
 		if (ca)
 			*ca = _ex0;
@@ -1478,6 +1478,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
 			*cb = _ex1;
 		return b;
 	}
+	ex bex = b.expand();
 	if (bex.is_zero()) {
 		if (ca)
 			*ca = _ex1;
@@ -1563,8 +1564,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
 			*cb = b;
 		return _ex1;
 	}
-	// move symbols which contained only in one of the polynomials
-	// to the end:
+	// move symbol contained only in one of the polynomials to the end:
 	rotate(sym_stats.begin(), vari, sym_stats.end());
 
 	sym_desc_vec::const_iterator var = sym_stats.begin();
@@ -1676,7 +1676,6 @@ static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb)
 			if (cb)
 				*cb = b;
 			return _ex1;
-			// XXX: do I need to check for p_gcd = -1?
 	}
 
 	// there are common factors:
@@ -1706,17 +1705,27 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
 	if (p.is_equal(b)) {
 		// a = p^n, b = p, gcd = p
 		if (ca)
-			*ca = pow(p, a.op(1) - 1);
+			*ca = pow(p, exp_a - 1);
 		if (cb)
 			*cb = _ex1;
 		return p;
-	} 
+	}
+	if (is_a<symbol>(p)) {
+		// Cancel trivial common factor
+		int ldeg_a = ex_to<numeric>(exp_a).to_int();
+		int ldeg_b = b.ldegree(p);
+		int min_ldeg = std::min(ldeg_a, ldeg_b);
+		if (min_ldeg > 0) {
+			ex common = pow(p, min_ldeg);
+			return gcd(pow(p, ldeg_a - min_ldeg), (b / common).expand(), ca, cb, false) * common;
+		}
+	}
 
 	ex p_co, bpart_co;
 	ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
 
-	// a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
 	if (p_gcd.is_equal(_ex1)) {
+		// a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
 		if (ca)
 			*ca = a;
 		if (cb)
@@ -1785,32 +1794,39 @@ ex lcm(const ex &a, const ex &b, bool check_args)
  *  @param a  multivariate polynomial over Z[X], treated here as univariate
  *            polynomial in x (needs not be expanded).
  *  @param x  variable to factor in
- *  @return   vector of factors sorted in ascending degree */
-static exvector sqrfree_yun(const ex &a, const symbol &x)
+ *  @return   vector of expairs (factor, exponent), sorted by exponents */
+static epvector sqrfree_yun(const ex &a, const symbol &x)
 {
-	exvector res;
 	ex w = a;
 	ex z = w.diff(x);
 	ex g = gcd(w, z);
 	if (g.is_zero()) {
-		return res;
+		return epvector{};
 	}
 	if (g.is_equal(_ex1)) {
-		res.push_back(a);
-		return res;
+		return epvector{expair(a, _ex1)};
 	}
-	ex y;
+	epvector results;
+	ex exponent = _ex0;
 	do {
 		w = quo(w, g, x);
 		if (w.is_zero()) {
-			return res;
+			return results;
+		}
+		z = quo(z, g, x) - w.diff(x);
+		exponent = exponent + 1;
+		if (w.is_equal(x)) {
+			// shortcut for x^n with n ∈ ℕ
+			exponent += quo(z, w.diff(x), x);
+			results.push_back(expair(w, exponent));
+			break;
 		}
-		y = quo(z, g, x);
-		z = y - w.diff(x);
 		g = gcd(w, z);
-		res.push_back(g);
+		if (!g.is_equal(_ex1)) {
+			results.push_back(expair(g, exponent));
+		}
 	} while (!z.is_zero());
-	return res;
+	return results;
 }
 
 
@@ -1878,30 +1894,28 @@ ex sqrfree(const ex &a, const lst &l)
 	const ex tmp = multiply_lcm(a,lcm);
 
 	// find the factors
-	exvector factors = sqrfree_yun(tmp, x);
+	epvector factors = sqrfree_yun(tmp, x);
 
-	// construct the next list of symbols with the first element popped
-	lst newargs = args;
-	newargs.remove_first();
+	// remove symbol x and proceed recursively with the remaining symbols
+	args.remove_first();
 
 	// recurse down the factors in remaining variables
-	if (newargs.nops()>0) {
+	if (args.nops()>0) {
 		for (auto & it : factors)
-			it = sqrfree(it, newargs);
+			it.rest = sqrfree(it.rest, args);
 	}
 
 	// Done with recursion, now construct the final result
 	ex result = _ex1;
-	int p = 1;
 	for (auto & it : factors)
-		result *= pow(it, p++);
+		result *= pow(it.rest, it.coeff);
 
 	// 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.
-	if (newargs.nops()>0)
-		result *= sqrfree(quo(tmp, result, x), newargs);
+	if (args.nops()>0)
+		result *= sqrfree(quo(tmp, result, x), args);
 	else
 		result *= quo(tmp, result, x);
 
@@ -1929,24 +1943,21 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
 //clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << endl;
 
 	// Factorize denominator and compute cofactors
-	exvector yun = sqrfree_yun(denom, x);
-//clog << "yun factors: " << exprseq(yun) << endl;
-	size_t num_yun = yun.size();
-	exvector factor; factor.reserve(num_yun);
-	exvector cofac; cofac.reserve(num_yun);
-	for (size_t i=0; i<num_yun; i++) {
-		if (!yun[i].is_equal(_ex1)) {
-			for (size_t j=0; j<=i; j++) {
-				factor.push_back(pow(yun[i], j+1));
-				ex prod = _ex1;
-				for (size_t k=0; k<num_yun; k++) {
-					if (k == i)
-						prod *= pow(yun[k], i-j);
-					else
-						prod *= pow(yun[k], k+1);
-				}
-				cofac.push_back(prod.expand());
+	epvector yun = sqrfree_yun(denom, 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++) {
+		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));
+			ex prod = _ex1;
+			for (size_t k=0; k<yun.size(); k++) {
+				if (yun[k].coeff == i_exponent)
+					prod *= pow(yun[k].rest, i_exponent-1-j);
+				else
+					prod *= pow(yun[k].rest, yun[k].coeff);
 			}
+			cofac.push_back(prod.expand());
 		}
 	}
 	size_t num_factors = factor.size();