X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=6870b77f4b7c8b589d05a1ceafd842cfadb82e88;hp=e97a7ce4a2c19c7f530881e5c1257e7562f20604;hb=acae7ab5a4dc94d1f54ba794f32f5764cdb4d704;hpb=1d09676118944532e6100c93383d1659b1253a60

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index e97a7ce4..6870b77f 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-2008 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2015 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
@@ -23,9 +23,6 @@
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#include <algorithm>
-#include <map>
-
 #include "normal.h"
 #include "basic.h"
 #include "ex.h"
@@ -44,6 +41,10 @@
 #include "pseries.h"
 #include "symbol.h"
 #include "utils.h"
+#include "polynomial/chinrem_gcd.h"
+
+#include <algorithm>
+#include <map>
 
 namespace GiNaC {
 
@@ -672,12 +673,13 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 			q = rem_i*power(ab, a_exp - 1);
 			return true;
 		}
-		for (int i=2; i < a_exp; i++) {
-			if (divide(power(ab, i), b, rem_i, false)) {
-				q = rem_i*power(ab, a_exp - i);
-				return true;
-			}
-		} // ... so we *really* need to expand expression.
+// code below is commented-out because it leads to a significant slowdown
+//		for (int i=2; i < a_exp; i++) {
+//			if (divide(power(ab, i), b, rem_i, false)) {
+//				q = rem_i*power(ab, a_exp - i);
+//				return true;
+//			}
+//		} // ... so we *really* need to expand expression.
 	}
 	
 	// Polynomial long division (recursive)
@@ -959,7 +961,7 @@ ex ex::content(const ex &x) const
 		return lcoeff * c / lcoeff.unit(x);
 	ex cont = _ex0;
 	for (int i=ldeg; i<=deg; i++)
-		cont = gcd(r.coeff(x, i), cont, NULL, NULL, false);
+		cont = gcd(r.coeff(x, i), cont, nullptr, nullptr, false);
 	return cont * c;
 }
 
@@ -1099,7 +1101,7 @@ static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
 	// Remove content from c and d, to be attached to GCD later
 	ex cont_c = c.content(x);
 	ex cont_d = d.content(x);
-	ex gamma = gcd(cont_c, cont_d, NULL, NULL, false);
+	ex gamma = gcd(cont_c, cont_d, nullptr, nullptr, false);
 	if (ddeg == 0)
 		return gamma;
 	c = c.primpart(x, cont_c);
@@ -1270,9 +1272,9 @@ class gcdheu_failed {};
  *
  *  @param a  first integer multivariate polynomial (expanded)
  *  @param b  second integer multivariate polynomial (expanded)
- *  @param ca  cofactor of polynomial a (returned), NULL to suppress
+ *  @param ca  cofactor of polynomial a (returned), nullptr to suppress
  *             calculation of cofactor
- *  @param cb  cofactor of polynomial b (returned), NULL to suppress
+ *  @param cb  cofactor of polynomial b (returned), nullptr to suppress
  *             calculation of cofactor
  *  @param var iterator to first element of vector of sym_desc structs
  *  @param res the GCD (returned)
@@ -1363,9 +1365,9 @@ static bool heur_gcd_z(ex& res, const ex &a, const ex &b, ex *ca, ex *cb,
  *
  *  @param a  first rational multivariate polynomial (expanded)
  *  @param b  second rational multivariate polynomial (expanded)
- *  @param ca  cofactor of polynomial a (returned), NULL to suppress
+ *  @param ca  cofactor of polynomial a (returned), nullptr to suppress
  *             calculation of cofactor
- *  @param cb  cofactor of polynomial b (returned), NULL to suppress
+ *  @param cb  cofactor of polynomial b (returned), nullptr to suppress
  *             calculation of cofactor
  *  @param var iterator to first element of vector of sym_desc structs
  *  @param res the GCD (returned)
@@ -1429,8 +1431,8 @@ static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb);
  *
  *  @param a  first multivariate polynomial
  *  @param b  second multivariate polynomial
- *  @param ca pointer to expression that will receive the cofactor of a, or NULL
- *  @param cb pointer to expression that will receive the cofactor of b, or NULL
+ *  @param ca pointer to expression that will receive the cofactor of a, or nullptr
+ *  @param cb pointer to expression that will receive the cofactor of b, or nullptr
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return the GCD as a new expression */
@@ -1533,7 +1535,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
 		if (ca)
 			*ca = ex_to<numeric>(aex)/g;
 		if (cb)
-	 		*cb = bex/g;
+			*cb = bex/g;
 		return g;
 	}
 
@@ -1622,8 +1624,15 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
 		}
 #endif
 	}
+	if (options & gcd_options::use_sr_gcd) {
+		g = sr_gcd(aex, bex, var);
+	} else {
+		exvector vars;
+		for (std::size_t n = sym_stats.size(); n-- != 0; )
+			vars.push_back(sym_stats[n].sym);
+		g = chinrem_gcd(aex, bex, vars);
+	}
 
-	g = sr_gcd(aex, bex, var);
 	if (g.is_equal(_ex1)) {
 		// Keep cofactors factored if possible
 		if (ca)
@@ -1639,109 +1648,91 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
 	return g;
 }
 
-static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
+// gcd helper to handle partially factored polynomials (to avoid expanding
+// large expressions). Both arguments should be powers.
+static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb)
 {
-	if (is_exactly_a<power>(a)) {
-		ex p = a.op(0);
-		const ex& exp_a = a.op(1);
-		if (is_exactly_a<power>(b)) {
-			ex pb = b.op(0);
-			const ex& exp_b = b.op(1);
-			if (p.is_equal(pb)) {
-				// a = p^n, b = p^m, gcd = p^min(n, m)
-				if (exp_a < exp_b) {
-					if (ca)
-						*ca = _ex1;
-					if (cb)
-						*cb = power(p, exp_b - exp_a);
-					return power(p, exp_a);
-				} else {
-					if (ca)
-						*ca = power(p, exp_a - exp_b);
-					if (cb)
-						*cb = _ex1;
-					return power(p, exp_b);
-				}
-			} else {
-				ex p_co, pb_co;
-				ex p_gcd = gcd(p, pb, &p_co, &pb_co, false);
-				if (p_gcd.is_equal(_ex1)) {
-					// a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==>
-					// gcd(a,b) = 1
-					if (ca)
-						*ca = a;
-					if (cb)
-						*cb = b;
-					return _ex1;
-					// XXX: do I need to check for p_gcd = -1?
-				} else {
-					// there are common factors:
-					// a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
-					// gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
-					if (exp_a < exp_b) {
-						return power(p_gcd, exp_a)*
-							gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false);
-					} else {
-						return power(p_gcd, exp_b)*
-							gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false);
-					}
-				} // p_gcd.is_equal(_ex1)
-			} // p.is_equal(pb)
-
-		} else {
-			if (p.is_equal(b)) {
-				// a = p^n, b = p, gcd = p
-				if (ca)
-					*ca = power(p, a.op(1) - 1);
-				if (cb)
-					*cb = _ex1;
-				return p;
-			} 
-
-			ex p_co, bpart_co;
-			ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
-
-			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)
-					*cb = b;
-				return _ex1;
-			} else {
-				// a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
-				return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false);
-			}
-		} // is_exactly_a<power>(b)
+	ex p = a.op(0);
+	const ex& exp_a = a.op(1);
+	ex pb = b.op(0);
+	const ex& exp_b = b.op(1);
 
-	} else if (is_exactly_a<power>(b)) {
-		ex p = b.op(0);
-		if (p.is_equal(a)) {
-			// a = p, b = p^n, gcd = p
+	// a = p^n, b = p^m, gcd = p^min(n, m)
+	if (p.is_equal(pb)) {
+		if (exp_a < exp_b) {
 			if (ca)
 				*ca = _ex1;
 			if (cb)
-				*cb = power(p, b.op(1) - 1);
-			return p;
+				*cb = power(p, exp_b - exp_a);
+			return power(p, exp_a);
+		} else {
+			if (ca)
+				*ca = power(p, exp_a - exp_b);
+			if (cb)
+				*cb = _ex1;
+			return power(p, exp_b);
 		}
+	}
 
-		ex p_co, apart_co;
-		const ex& exp_b(b.op(1));
-		ex p_gcd = gcd(a, p, &apart_co, &p_co, false);
-		if (p_gcd.is_equal(_ex1)) {
-			// b=p(x)^n, gcd(a, p) = 1 ==> gcd(a, b) == 1
+	ex p_co, pb_co;
+	ex p_gcd = gcd(p, pb, &p_co, &pb_co, false);
+	// a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==> gcd(a,b) = 1
+	if (p_gcd.is_equal(_ex1)) {
 			if (ca)
 				*ca = a;
 			if (cb)
 				*cb = b;
 			return _ex1;
-		} else {
-			// there are common factors:
-			// a(x) = g(x) A(x), b(x) = g(x)^n B(x)^n ==> gcd = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
+			// XXX: do I need to check for p_gcd = -1?
+	}
+
+	// there are common factors:
+	// a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
+	// gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
+	if (exp_a < exp_b) {
+		ex pg =  gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false);
+		return power(p_gcd, exp_a)*pg;
+	} else {
+		ex pg = gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false);
+		return power(p_gcd, exp_b)*pg;
+	}
+}
+
+static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
+{
+	if (is_exactly_a<power>(a) && is_exactly_a<power>(b))
+		return gcd_pf_pow_pow(a, b, ca, cb);
+
+	if (is_exactly_a<power>(b) && (! is_exactly_a<power>(a)))
+		return gcd_pf_pow(b, a, cb, ca);
+
+	GINAC_ASSERT(is_exactly_a<power>(a));
+
+	ex p = a.op(0);
+	const ex& exp_a = a.op(1);
+	if (p.is_equal(b)) {
+		// a = p^n, b = p, gcd = p
+		if (ca)
+			*ca = power(p, a.op(1) - 1);
+		if (cb)
+			*cb = _ex1;
+		return p;
+	} 
 
-			return p_gcd*gcd(apart_co, power(p_gcd, exp_b-1)*power(p_co, exp_b), ca, cb, false);
-		} // p_gcd.is_equal(_ex1)
+	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)) {
+		if (ca)
+			*ca = a;
+		if (cb)
+			*cb = b;
+		return _ex1;
 	}
+	// a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
+	ex rg = gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false);
+	return p_gcd*rg;
 }
 
 static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb)
@@ -1921,7 +1912,7 @@ ex sqrfree(const ex &a, const lst &l)
 	else
 		result *= quo(tmp, result, x);
 
-	// Put in the reational overall factor again and return
+	// Put in the rational overall factor again and return
 	return result *	lcm.inverse();
 }
 
@@ -2015,16 +2006,18 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
  *  @see ex::normal */
 static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup)
 {
+	// Since the repl contains replaced expressions we should search for them
+	ex e_replaced = e.subs(repl, subs_options::no_pattern);
+
 	// Expression already replaced? Then return the assigned symbol
-	exmap::const_iterator it = rev_lookup.find(e);
+	exmap::const_iterator it = rev_lookup.find(e_replaced);
 	if (it != rev_lookup.end())
 		return it->second;
-	
+
 	// Otherwise create new symbol and add to list, taking care that the
 	// replacement expression doesn't itself contain symbols from repl,
 	// because subs() is not recursive
 	ex es = (new symbol)->setflag(status_flags::dynallocated);
-	ex e_replaced = e.subs(repl, subs_options::no_pattern);
 	repl.insert(std::make_pair(es, e_replaced));
 	rev_lookup.insert(std::make_pair(e_replaced, es));
 	return es;
@@ -2037,16 +2030,18 @@ static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup)
  *  @see basic::to_polynomial */
 static ex replace_with_symbol(const ex & e, exmap & repl)
 {
+	// Since the repl contains replaced expressions we should search for them
+	ex e_replaced = e.subs(repl, subs_options::no_pattern);
+
 	// Expression already replaced? Then return the assigned symbol
 	for (exmap::const_iterator it = repl.begin(); it != repl.end(); ++it)
-		if (it->second.is_equal(e))
+		if (it->second.is_equal(e_replaced))
 			return it->first;
-	
+
 	// Otherwise create new symbol and add to list, taking care that the
 	// replacement expression doesn't itself contain symbols from repl,
 	// because subs() is not recursive
 	ex es = (new symbol)->setflag(status_flags::dynallocated);
-	ex e_replaced = e.subs(repl, subs_options::no_pattern);
 	repl.insert(std::make_pair(es, e_replaced));
 	return es;
 }
@@ -2217,7 +2212,7 @@ ex add::normal(exmap & repl, exmap & rev_lookup, int level) const
 			num_it++; den_it++;
 		}
 
-		// Additiion of two fractions, taking advantage of the fact that
+		// Addition of two fractions, taking advantage of the fact that
 		// the heuristic GCD algorithm computes the cofactors at no extra cost
 		ex co_den1, co_den2;
 		ex g = gcd(den, next_den, &co_den1, &co_den2, false);
@@ -2403,7 +2398,7 @@ ex ex::denom() const
 		return e.op(1).subs(repl, subs_options::no_pattern);
 }
 
-/** Get numerator and denominator of an expression. If the expresison is not
+/** Get numerator and denominator of an expression. If the expression is not
  *  of the normal form "numerator/denominator", it is first converted to this
  *  form and then a list [numerator, denominator] is returned.
  *
@@ -2592,10 +2587,10 @@ ex expairseq::to_rational(exmap & repl) const
 	}
 	ex oc = overall_coeff.to_rational(repl);
 	if (oc.info(info_flags::numeric))
-		return thisexpairseq(s, overall_coeff);
+		return thisexpairseq(std::move(s), overall_coeff);
 	else
 		s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1));
-	return thisexpairseq(s, default_overall_coeff());
+	return thisexpairseq(std::move(s), default_overall_coeff());
 }
 
 /** Implementation of ex::to_polynomial() for expairseqs. */
@@ -2610,10 +2605,10 @@ ex expairseq::to_polynomial(exmap & repl) const
 	}
 	ex oc = overall_coeff.to_polynomial(repl);
 	if (oc.info(info_flags::numeric))
-		return thisexpairseq(s, overall_coeff);
+		return thisexpairseq(std::move(s), overall_coeff);
 	else
 		s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1));
-	return thisexpairseq(s, default_overall_coeff());
+	return thisexpairseq(std::move(s), default_overall_coeff());
 }