X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=5d044fc6f40f6da06a6aaff3b12c670115cd7c3a;hp=b2a4e8c8d0cfdfe2204799fefe96b2710771cb24;hb=c77689e7ac8d8f4dbca0f337b6e9acf2419010ff;hpb=9c822131b5a057640b040af6df828e0d1ed0222e

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index b2a4e8c8..5d044fc6 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-2005 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2008 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
@@ -1268,17 +1268,19 @@ class gcdheu_failed {};
  *  polynomials and an iterator to the first element of the sym_desc vector
  *  passed in. This function is used internally by gcd().
  *
- *  @param a  first multivariate polynomial (expanded)
- *  @param b  second multivariate polynomial (expanded)
+ *  @param a  first integer multivariate polynomial (expanded)
+ *  @param b  second integer multivariate polynomial (expanded)
  *  @param ca  cofactor of polynomial a (returned), NULL to suppress
  *             calculation of cofactor
  *  @param cb  cofactor of polynomial b (returned), NULL to suppress
  *             calculation of cofactor
  *  @param var iterator to first element of vector of sym_desc structs
- *  @return the GCD as a new expression
+ *  @param res the GCD (returned)
+ *  @return true if GCD was computed, false otherwise.
  *  @see gcd
  *  @exception gcdheu_failed() */
-static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
+static bool heur_gcd_z(ex& res, const ex &a, const ex &b, ex *ca, ex *cb,
+	               sym_desc_vec::const_iterator var)
 {
 #if STATISTICS
 	heur_gcd_called++;
@@ -1286,7 +1288,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 
 	// Algorithm only works for non-vanishing input polynomials
 	if (a.is_zero() || b.is_zero())
-		return (new fail())->setflag(status_flags::dynallocated);
+		return false;
 
 	// GCD of two numeric values -> CLN
 	if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
@@ -1295,7 +1297,8 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 			*ca = ex_to<numeric>(a) / g;
 		if (cb)
 			*cb = ex_to<numeric>(b) / g;
-		return g;
+		res = g;
+		return true;
 	}
 
 	// The first symbol is our main variable
@@ -1325,9 +1328,13 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 
 		// Apply evaluation homomorphism and calculate GCD
 		ex cp, cq;
-		ex gamma = heur_gcd(p.subs(x == xi, subs_options::no_pattern), q.subs(x == xi, subs_options::no_pattern), &cp, &cq, var+1).expand();
-		if (!is_exactly_a<fail>(gamma)) {
-
+		ex gamma;
+		bool found = heur_gcd_z(gamma,
+			                p.subs(x == xi, subs_options::no_pattern),
+			                q.subs(x == xi, subs_options::no_pattern),
+				        &cp, &cq, var+1);
+		if (found) {
+			gamma = gamma.expand();
 			// Reconstruct polynomial from GCD of mapped polynomials
 			ex g = interpolate(gamma, xi, x, maxdeg);
 
@@ -1338,17 +1345,84 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 			ex dummy;
 			if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) {
 				g *= gc;
-				return g;
+				res = g;
+				return true;
 			}
 		}
 
 		// Next evaluation point
 		xi = iquo(xi * isqrt(isqrt(xi)) * numeric(73794), numeric(27011));
 	}
-	return (new fail())->setflag(status_flags::dynallocated);
+	return false;
+}
+
+/** Compute GCD of multivariate polynomials using the heuristic GCD algorithm.
+ *  get_symbol_stats() must have been called previously with the input
+ *  polynomials and an iterator to the first element of the sym_desc vector
+ *  passed in. This function is used internally by gcd().
+ *
+ *  @param a  first rational multivariate polynomial (expanded)
+ *  @param b  second rational multivariate polynomial (expanded)
+ *  @param ca  cofactor of polynomial a (returned), NULL to suppress
+ *             calculation of cofactor
+ *  @param cb  cofactor of polynomial b (returned), NULL to suppress
+ *             calculation of cofactor
+ *  @param var iterator to first element of vector of sym_desc structs
+ *  @param res the GCD (returned)
+ *  @return true if GCD was computed, false otherwise.
+ *  @see heur_gcd_z
+ *  @see gcd
+ */
+static bool heur_gcd(ex& res, const ex& a, const ex& b, ex *ca, ex *cb,
+	             sym_desc_vec::const_iterator var)
+{
+	if (a.info(info_flags::integer_polynomial) && 
+	    b.info(info_flags::integer_polynomial)) {
+		try {
+			return heur_gcd_z(res, a, b, ca, cb, var);
+		} catch (gcdheu_failed) {
+			return false;
+		}
+	}
+
+	// convert polynomials to Z[X]
+	const numeric a_lcm = lcm_of_coefficients_denominators(a);
+	const numeric ab_lcm = lcmcoeff(b, a_lcm);
+
+	const ex ai = a*ab_lcm;
+	const ex bi = b*ab_lcm;
+	if (!ai.info(info_flags::integer_polynomial))
+		throw std::logic_error("heur_gcd: not an integer polynomial [1]");
+
+	if (!bi.info(info_flags::integer_polynomial))
+		throw std::logic_error("heur_gcd: not an integer polynomial [2]");
+
+	bool found = false;
+	try {
+		found = heur_gcd_z(res, ai, bi, ca, cb, var);
+	} catch (gcdheu_failed) {
+		return false;
+	}
+	
+	// GCD is not unique, it's defined up to a unit (i.e. invertible
+	// element). If the coefficient ring is a field, every its element is
+	// invertible, so one can multiply the polynomial GCD with any element
+	// of the coefficient field. We use this ambiguity to make cofactors
+	// integer polynomials.
+	if (found)
+		res /= ab_lcm;
+	return found;
 }
 
 
+// gcd helper to handle partially factored polynomials (to avoid expanding
+// large expressions). At least one of the arguments should be a power.
+static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args);
+
+// gcd helper to handle partially factored polynomials (to avoid expanding
+// large expressions). At least one of the arguments should be a product.
+static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args);
+
 /** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X)
  *  and b(X) in Z[X]. Optionally also compute the cofactors of a and b,
  *  defined by a = ca * gcd(a, b) and b = cb * gcd(a, b).
@@ -1360,7 +1434,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return the GCD as a new expression */
-ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
+ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
 {
 #if STATISTICS
 	gcd_called++;
@@ -1391,48 +1465,182 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
 	}
 
 	// Partially factored cases (to avoid expanding large expressions)
-	if (is_exactly_a<mul>(a)) {
-		if (is_exactly_a<mul>(b) && b.nops() > a.nops())
-			goto factored_b;
-factored_a:
-		size_t num = a.nops();
-		exvector g; g.reserve(num);
-		exvector acc_ca; acc_ca.reserve(num);
-		ex part_b = b;
-		for (size_t i=0; i<num; i++) {
-			ex part_ca, part_cb;
-			g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args));
-			acc_ca.push_back(part_ca);
-			part_b = part_cb;
+	if (!(options & gcd_options::no_part_factored)) {
+		if (is_exactly_a<mul>(a) || is_exactly_a<mul>(b))
+			return gcd_pf_mul(a, b, ca, cb, check_args);
+#if FAST_COMPARE
+		if (is_exactly_a<power>(a) || is_exactly_a<power>(b))
+			return gcd_pf_pow(a, b, ca, cb, check_args);
+#endif
+	}
+
+	// Some trivial cases
+	ex aex = a.expand(), bex = b.expand();
+	if (aex.is_zero()) {
+		if (ca)
+			*ca = _ex0;
+		if (cb)
+			*cb = _ex1;
+		return b;
+	}
+	if (bex.is_zero()) {
+		if (ca)
+			*ca = _ex1;
+		if (cb)
+			*cb = _ex0;
+		return a;
+	}
+	if (aex.is_equal(_ex1) || bex.is_equal(_ex1)) {
+		if (ca)
+			*ca = a;
+		if (cb)
+			*cb = b;
+		return _ex1;
+	}
+#if FAST_COMPARE
+	if (a.is_equal(b)) {
+		if (ca)
+			*ca = _ex1;
+		if (cb)
+			*cb = _ex1;
+		return a;
+	}
+#endif
+
+	if (is_a<symbol>(aex)) {
+		if (! bex.subs(aex==_ex0, subs_options::no_pattern).is_zero()) {
+			if (ca)
+				*ca = a;
+			if (cb)
+				*cb = b;
+			return _ex1;
 		}
+	}
+
+	if (is_a<symbol>(bex)) {
+		if (! aex.subs(bex==_ex0, subs_options::no_pattern).is_zero()) {
+			if (ca)
+				*ca = a;
+			if (cb)
+				*cb = b;
+			return _ex1;
+		}
+	}
+
+	if (is_exactly_a<numeric>(aex)) {
+		numeric bcont = bex.integer_content();
+		numeric g = gcd(ex_to<numeric>(aex), bcont);
 		if (ca)
-			*ca = (new mul(acc_ca))->setflag(status_flags::dynallocated);
+			*ca = ex_to<numeric>(aex)/g;
 		if (cb)
-			*cb = part_b;
-		return (new mul(g))->setflag(status_flags::dynallocated);
-	} else if (is_exactly_a<mul>(b)) {
-		if (is_exactly_a<mul>(a) && a.nops() > b.nops())
-			goto factored_a;
-factored_b:
-		size_t num = b.nops();
-		exvector g; g.reserve(num);
-		exvector acc_cb; acc_cb.reserve(num);
-		ex part_a = a;
-		for (size_t i=0; i<num; i++) {
-			ex part_ca, part_cb;
-			g.push_back(gcd(part_a, b.op(i), &part_ca, &part_cb, check_args));
-			acc_cb.push_back(part_cb);
-			part_a = part_ca;
+	 		*cb = bex/g;
+		return g;
+	}
+
+	if (is_exactly_a<numeric>(bex)) {
+		numeric acont = aex.integer_content();
+		numeric g = gcd(ex_to<numeric>(bex), acont);
+		if (ca)
+			*ca = aex/g;
+		if (cb)
+			*cb = ex_to<numeric>(bex)/g;
+		return g;
+	}
+
+	// Gather symbol statistics
+	sym_desc_vec sym_stats;
+	get_symbol_stats(a, b, sym_stats);
+
+	// The symbol with least degree which is contained in both polynomials
+	// is our main variable
+	sym_desc_vec::iterator vari = sym_stats.begin();
+	while ((vari != sym_stats.end()) && 
+	       (((vari->ldeg_b == 0) && (vari->deg_b == 0)) ||
+	        ((vari->ldeg_a == 0) && (vari->deg_a == 0))))
+		vari++;
+
+	// No common symbols at all, just return 1:
+	if (vari == sym_stats.end()) {
+		// N.B: keep cofactors factored
+		if (ca)
+			*ca = a;
+		if (cb)
+			*cb = b;
+		return _ex1;
+	}
+	// move symbols which 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();
+	const ex &x = var->sym;
+
+	// Cancel trivial common factor
+	int ldeg_a = var->ldeg_a;
+	int ldeg_b = var->ldeg_b;
+	int min_ldeg = std::min(ldeg_a,ldeg_b);
+	if (min_ldeg > 0) {
+		ex common = power(x, min_ldeg);
+		return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common;
+	}
+
+	// Try to eliminate variables
+	if (var->deg_a == 0 && var->deg_b != 0 ) {
+		ex bex_u, bex_c, bex_p;
+		bex.unitcontprim(x, bex_u, bex_c, bex_p);
+		ex g = gcd(aex, bex_c, ca, cb, false);
+		if (cb)
+			*cb *= bex_u * bex_p;
+		return g;
+	} else if (var->deg_b == 0 && var->deg_a != 0) {
+		ex aex_u, aex_c, aex_p;
+		aex.unitcontprim(x, aex_u, aex_c, aex_p);
+		ex g = gcd(aex_c, bex, ca, cb, false);
+		if (ca)
+			*ca *= aex_u * aex_p;
+		return g;
+	}
+
+	// Try heuristic algorithm first, fall back to PRS if that failed
+	ex g;
+	if (!(options & gcd_options::no_heur_gcd)) {
+		bool found = heur_gcd(g, aex, bex, ca, cb, var);
+		if (found) {
+			// heur_gcd have already computed cofactors...
+			if (g.is_equal(_ex1)) {
+				// ... but we want to keep them factored if possible.
+				if (ca)
+					*ca = a;
+				if (cb)
+					*cb = b;
+			}
+			return g;
 		}
+#if STATISTICS
+		else {
+			heur_gcd_failed++;
+		}
+#endif
+	}
+
+	g = sr_gcd(aex, bex, var);
+	if (g.is_equal(_ex1)) {
+		// Keep cofactors factored if possible
 		if (ca)
-			*ca = part_a;
+			*ca = a;
 		if (cb)
-			*cb = (new mul(acc_cb))->setflag(status_flags::dynallocated);
-		return (new mul(g))->setflag(status_flags::dynallocated);
+			*cb = b;
+	} else {
+		if (ca)
+			divide(aex, g, *ca, false);
+		if (cb)
+			divide(bex, g, *cb, false);
 	}
+	return g;
+}
 
-#if FAST_COMPARE
-	// Input polynomials of the form poly^n are sometimes also trivial
+static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args)
+{
 	if (is_exactly_a<power>(a)) {
 		ex p = a.op(0);
 		const ex& exp_a = a.op(1);
@@ -1534,133 +1742,51 @@ factored_b:
 			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)
 	}
-#endif
+}
 
-	// Some trivial cases
-	ex aex = a.expand(), bex = b.expand();
-	if (aex.is_zero()) {
-		if (ca)
-			*ca = _ex0;
-		if (cb)
-			*cb = _ex1;
-		return b;
-	}
-	if (bex.is_zero()) {
-		if (ca)
-			*ca = _ex1;
-		if (cb)
-			*cb = _ex0;
-		return a;
-	}
-	if (aex.is_equal(_ex1) || bex.is_equal(_ex1)) {
-		if (ca)
-			*ca = a;
-		if (cb)
-			*cb = b;
-		return _ex1;
-	}
-#if FAST_COMPARE
-	if (a.is_equal(b)) {
+static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args)
+{
+	if (is_exactly_a<mul>(a)) {
+		if (is_exactly_a<mul>(b) && b.nops() > a.nops())
+			goto factored_b;
+factored_a:
+		size_t num = a.nops();
+		exvector g; g.reserve(num);
+		exvector acc_ca; acc_ca.reserve(num);
+		ex part_b = b;
+		for (size_t i=0; i<num; i++) {
+			ex part_ca, part_cb;
+			g.push_back(gcd(a.op(i), part_b, &part_ca, &part_cb, check_args));
+			acc_ca.push_back(part_ca);
+			part_b = part_cb;
+		}
 		if (ca)
-			*ca = _ex1;
+			*ca = (new mul(acc_ca))->setflag(status_flags::dynallocated);
 		if (cb)
-			*cb = _ex1;
-		return a;
-	}
-#endif
-
-	if (is_a<symbol>(aex)) {
-		if (! bex.subs(aex==_ex0, subs_options::no_pattern).is_zero()) {
-			if (ca)
-				*ca = a;
-			if (cb)
-				*cb = b;
-			return _ex1;
-		}
-	}
-
-	if (is_a<symbol>(bex)) {
-		if (! aex.subs(bex==_ex0, subs_options::no_pattern).is_zero()) {
-			if (ca)
-				*ca = a;
-			if (cb)
-				*cb = b;
-			return _ex1;
+			*cb = part_b;
+		return (new mul(g))->setflag(status_flags::dynallocated);
+	} else if (is_exactly_a<mul>(b)) {
+		if (is_exactly_a<mul>(a) && a.nops() > b.nops())
+			goto factored_a;
+factored_b:
+		size_t num = b.nops();
+		exvector g; g.reserve(num);
+		exvector acc_cb; acc_cb.reserve(num);
+		ex part_a = a;
+		for (size_t i=0; i<num; i++) {
+			ex part_ca, part_cb;
+			g.push_back(gcd(part_a, b.op(i), &part_ca, &part_cb, check_args));
+			acc_cb.push_back(part_cb);
+			part_a = part_ca;
 		}
-	}
-
-	// Gather symbol statistics
-	sym_desc_vec sym_stats;
-	get_symbol_stats(a, b, sym_stats);
-
-	// The symbol with least degree is our main variable
-	sym_desc_vec::const_iterator var = sym_stats.begin();
-	const ex &x = var->sym;
-
-	// Cancel trivial common factor
-	int ldeg_a = var->ldeg_a;
-	int ldeg_b = var->ldeg_b;
-	int min_ldeg = std::min(ldeg_a,ldeg_b);
-	if (min_ldeg > 0) {
-		ex common = power(x, min_ldeg);
-		return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common;
-	}
-
-	// Try to eliminate variables
-	if (var->deg_a == 0) {
-		ex bex_u, bex_c, bex_p;
-		bex.unitcontprim(x, bex_u, bex_c, bex_p);
-		ex g = gcd(aex, bex_c, ca, cb, false);
-		if (cb)
-			*cb *= bex_u * bex_p;
-		return g;
-	} else if (var->deg_b == 0) {
-		ex aex_u, aex_c, aex_p;
-		aex.unitcontprim(x, aex_u, aex_c, aex_p);
-		ex g = gcd(aex_c, bex, ca, cb, false);
 		if (ca)
-			*ca *= aex_u * aex_p;
-		return g;
-	}
-
-	// Try heuristic algorithm first, fall back to PRS if that failed
-	ex g;
-	try {
-		g = heur_gcd(aex, bex, ca, cb, var);
-	} catch (gcdheu_failed) {
-		g = fail();
-	}
-	if (is_exactly_a<fail>(g)) {
-#if STATISTICS
-		heur_gcd_failed++;
-#endif
-		g = sr_gcd(aex, bex, var);
-		if (g.is_equal(_ex1)) {
-			// Keep cofactors factored if possible
-			if (ca)
-				*ca = a;
-			if (cb)
-				*cb = b;
-		} else {
-			if (ca)
-				divide(aex, g, *ca, false);
-			if (cb)
-				divide(bex, g, *cb, false);
-		}
-	} else {
-		if (g.is_equal(_ex1)) {
-			// Keep cofactors factored if possible
-			if (ca)
-				*ca = a;
-			if (cb)
-				*cb = b;
-		}
+			*ca = part_a;
+		if (cb)
+			*cb = (new mul(acc_cb))->setflag(status_flags::dynallocated);
+		return (new mul(g))->setflag(status_flags::dynallocated);
 	}
-
-	return g;
 }
 
-
 /** Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X].
  *
  *  @param a  first multivariate polynomial