X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=66e682c9a7da88f244a8e6cb0592db81682b2829;hp=9e9bde18f0f65a4a9583091986a47250f18312a8;hb=47923e5885d0e437d6cbe257c25f9bca757ea001;hpb=da64e515abf7243bc4c84ca3631470931c4e6691

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 9e9bde18..66e682c9 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-2007 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
@@ -233,14 +233,14 @@ static numeric lcmcoeff(const ex &e, const numeric &l)
 	if (e.info(info_flags::rational))
 		return lcm(ex_to<numeric>(e).denom(), l);
 	else if (is_exactly_a<add>(e)) {
-		numeric c = _num1;
+		numeric c = *_num1_p;
 		for (size_t i=0; i<e.nops(); i++)
 			c = lcmcoeff(e.op(i), c);
 		return lcm(c, l);
 	} else if (is_exactly_a<mul>(e)) {
-		numeric c = _num1;
+		numeric c = *_num1_p;
 		for (size_t i=0; i<e.nops(); i++)
-			c *= lcmcoeff(e.op(i), _num1);
+			c *= lcmcoeff(e.op(i), *_num1_p);
 		return lcm(c, l);
 	} else if (is_exactly_a<power>(e)) {
 		if (is_a<symbol>(e.op(0)))
@@ -260,7 +260,7 @@ static numeric lcmcoeff(const ex &e, const numeric &l)
  *  @return LCM of denominators of coefficients */
 static numeric lcm_of_coefficients_denominators(const ex &e)
 {
-	return lcmcoeff(e, _num1);
+	return lcmcoeff(e, *_num1_p);
 }
 
 /** Bring polynomial from Q[X] to Z[X] by multiplying in the previously
@@ -273,9 +273,9 @@ static ex multiply_lcm(const ex &e, const numeric &lcm)
 	if (is_exactly_a<mul>(e)) {
 		size_t num = e.nops();
 		exvector v; v.reserve(num + 1);
-		numeric lcm_accum = _num1;
+		numeric lcm_accum = *_num1_p;
 		for (size_t i=0; i<num; i++) {
-			numeric op_lcm = lcmcoeff(e.op(i), _num1);
+			numeric op_lcm = lcmcoeff(e.op(i), *_num1_p);
 			v.push_back(multiply_lcm(e.op(i), op_lcm));
 			lcm_accum *= op_lcm;
 		}
@@ -310,7 +310,7 @@ numeric ex::integer_content() const
 
 numeric basic::integer_content() const
 {
-	return _num1;
+	return *_num1_p;
 }
 
 numeric numeric::integer_content() const
@@ -322,7 +322,7 @@ numeric add::integer_content() const
 {
 	epvector::const_iterator it = seq.begin();
 	epvector::const_iterator itend = seq.end();
-	numeric c = _num0, l = _num1;
+	numeric c = *_num0_p, l = *_num1_p;
 	while (it != itend) {
 		GINAC_ASSERT(!is_exactly_a<numeric>(it->rest));
 		GINAC_ASSERT(is_exactly_a<numeric>(it->coeff));
@@ -614,6 +614,72 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 	if (!get_first_symbol(a, x) && !get_first_symbol(b, x))
 		throw(std::invalid_argument("invalid expression in divide()"));
 
+	// Try to avoid expanding partially factored expressions.
+	if (is_exactly_a<mul>(b)) {
+	// Divide sequentially by each term
+		ex rem_new, rem_old = a;
+		for (size_t i=0; i < b.nops(); i++) {
+			if (! divide(rem_old, b.op(i), rem_new, false))
+				return false;
+			rem_old = rem_new;
+		}
+		q = rem_new;
+		return true;
+	} else if (is_exactly_a<power>(b)) {
+		const ex& bb(b.op(0));
+		int exp_b = ex_to<numeric>(b.op(1)).to_int();
+		ex rem_new, rem_old = a;
+		for (int i=exp_b; i>0; i--) {
+			if (! divide(rem_old, bb, rem_new, false))
+				return false;
+			rem_old = rem_new;
+		}
+		q = rem_new;
+		return true;
+	} 
+	
+	if (is_exactly_a<mul>(a)) {
+		// Divide sequentially each term. If some term in a is divisible 
+		// by b we are done... and if not, we can't really say anything.
+		size_t i;
+		ex rem_i;
+		bool divisible_p = false;
+		for (i=0; i < a.nops(); ++i) {
+			if (divide(a.op(i), b, rem_i, false)) {
+				divisible_p = true;
+				break;
+			}
+		}
+		if (divisible_p) {
+			exvector resv;
+			resv.reserve(a.nops());
+			for (size_t j=0; j < a.nops(); j++) {
+				if (j==i)
+					resv.push_back(rem_i);
+				else
+					resv.push_back(a.op(j));
+			}
+			q = (new mul(resv))->setflag(status_flags::dynallocated);
+			return true;
+		}
+	} else if (is_exactly_a<power>(a)) {
+		// The base itself might be divisible by b, in that case we don't
+		// need to expand a
+		const ex& ab(a.op(0));
+		int a_exp = ex_to<numeric>(a.op(1)).to_int();
+		ex rem_i;
+		if (divide(ab, b, rem_i, false)) {
+			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.
+	}
+	
 	// Polynomial long division (recursive)
 	ex r = a.expand();
 	if (r.is_zero()) {
@@ -714,6 +780,31 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 	}
 #endif
 
+	if (is_exactly_a<power>(b)) {
+		const ex& bb(b.op(0));
+		ex qbar = a;
+		int exp_b = ex_to<numeric>(b.op(1)).to_int();
+		for (int i=exp_b; i>0; i--) {
+			if (!divide_in_z(qbar, bb, q, var))
+				return false;
+			qbar = q;
+		}
+		return true;
+	}
+
+	if (is_exactly_a<mul>(b)) {
+		ex qbar = a;
+		for (const_iterator itrb = b.begin(); itrb != b.end(); ++itrb) {
+			sym_desc_vec sym_stats;
+			get_symbol_stats(a, *itrb, sym_stats);
+			if (!divide_in_z(qbar, *itrb, q, sym_stats.begin()))
+				return false;
+
+			qbar = q;
+		}
+		return true;
+	}
+
 	// Main symbol
 	const ex &x = var->sym;
 
@@ -730,24 +821,24 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
 	// Compute values at evaluation points 0..adeg
 	vector<numeric> alpha; alpha.reserve(adeg + 1);
 	exvector u; u.reserve(adeg + 1);
-	numeric point = _num0;
+	numeric point = *_num0_p;
 	ex c;
 	for (i=0; i<=adeg; i++) {
 		ex bs = b.subs(x == point, subs_options::no_pattern);
 		while (bs.is_zero()) {
-			point += _num1;
+			point += *_num1_p;
 			bs = b.subs(x == point, subs_options::no_pattern);
 		}
 		if (!divide_in_z(a.subs(x == point, subs_options::no_pattern), bs, c, var+1))
 			return false;
 		alpha.push_back(point);
 		u.push_back(c);
-		point += _num1;
+		point += *_num1_p;
 	}
 
 	// Compute inverses
 	vector<numeric> rcp; rcp.reserve(adeg + 1);
-	rcp.push_back(_num0);
+	rcp.push_back(*_num0_p);
 	for (k=1; k<=adeg; k++) {
 		numeric product = alpha[k] - alpha[0];
 		for (i=1; i<k; i++)
@@ -1062,7 +1153,7 @@ numeric ex::max_coefficient() const
  *  @see heur_gcd */
 numeric basic::max_coefficient() const
 {
-	return _num1;
+	return *_num1_p;
 }
 
 numeric numeric::max_coefficient() const
@@ -1222,9 +1313,9 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 	numeric mq = q.max_coefficient();
 	numeric xi;
 	if (mp > mq)
-		xi = mq * _num2 + _num2;
+		xi = mq * (*_num2_p) + (*_num2_p);
 	else
-		xi = mp * _num2 + _num2;
+		xi = mp * (*_num2_p) + (*_num2_p);
 
 	// 6 tries maximum
 	for (int t=0; t<6; t++) {
@@ -1247,11 +1338,7 @@ 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;
-				ex lc = g.lcoeff(x);
-				if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
-					return -g;
-				else
-					return g;
+				return g;
 			}
 		}
 
@@ -1348,10 +1435,12 @@ factored_b:
 	// Input polynomials of the form poly^n are sometimes also trivial
 	if (is_exactly_a<power>(a)) {
 		ex p = a.op(0);
+		const ex& exp_a = a.op(1);
 		if (is_exactly_a<power>(b)) {
-			if (p.is_equal(b.op(0))) {
+			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)
-				ex exp_a = a.op(1), exp_b = b.op(1);
 				if (exp_a < exp_b) {
 					if (ca)
 						*ca = _ex1;
@@ -1365,7 +1454,32 @@ factored_b:
 						*cb = _ex1;
 					return power(p, exp_b);
 				}
-			}
+			} else {
+				ex p_co, pb_co;
+				ex p_gcd = gcd(p, pb, &p_co, &pb_co, check_args);
+				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
@@ -1374,8 +1488,24 @@ factored_b:
 				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)
+
 	} else if (is_exactly_a<power>(b)) {
 		ex p = b.op(0);
 		if (p.is_equal(a)) {
@@ -1386,6 +1516,23 @@ factored_b:
 				*cb = power(p, b.op(1) - 1);
 			return p;
 		}
+
+		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
+			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))
+
+			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
 
@@ -1422,11 +1569,71 @@ factored_b:
 	}
 #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 = ex_to<numeric>(aex)/g;
+		if (cb)
+	 		*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 is our main variable
+	// 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;
 
@@ -1440,14 +1647,14 @@ factored_b:
 	}
 
 	// Try to eliminate variables
-	if (var->deg_a == 0) {
+	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) {
+	} 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);
@@ -1842,7 +2049,7 @@ static ex frac_cancel(const ex &n, const ex &d)
 {
 	ex num = n;
 	ex den = d;
-	numeric pre_factor = _num1;
+	numeric pre_factor = *_num1_p;
 
 //std::clog << "frac_cancel num = " << num << ", den = " << den << std::endl;
 
@@ -2287,6 +2494,17 @@ ex power::to_polynomial(exmap & repl) const
 {
 	if (exponent.info(info_flags::posint))
 		return power(basis.to_rational(repl), exponent);
+	else if (exponent.info(info_flags::negint))
+	{
+		ex basis_pref = collect_common_factors(basis);
+		if (is_exactly_a<mul>(basis_pref) || is_exactly_a<power>(basis_pref)) {
+			// (A*B)^n will be automagically transformed to A^n*B^n
+			ex t = power(basis_pref, exponent);
+			return t.to_polynomial(repl);
+		}
+		else
+			return power(replace_with_symbol(power(basis, _ex_1), repl), -exponent);
+	} 
 	else
 		return replace_with_symbol(*this, repl);
 }
@@ -2344,7 +2562,7 @@ static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
 		for (size_t i=0; i<num; i++) {
 			ex x = e.op(i).to_polynomial(repl);
 
-			if (is_exactly_a<add>(x) || is_exactly_a<mul>(x)) {
+			if (is_exactly_a<add>(x) || is_exactly_a<mul>(x) || is_a<power>(x)) {
 				ex f = 1;
 				x = find_common_factor(x, f, repl);
 				x *= f;
@@ -2403,8 +2621,16 @@ term_done:	;
 		return (new mul(v))->setflag(status_flags::dynallocated);
 
 	} else if (is_exactly_a<power>(e)) {
-
-		return e.to_polynomial(repl);
+		const ex e_exp(e.op(1));
+		if (e_exp.info(info_flags::integer)) {
+			ex eb = e.op(0).to_polynomial(repl);
+			ex factor_local(_ex1);
+			ex pre_res = find_common_factor(eb, factor_local, repl);
+			factor *= power(factor_local, e_exp);
+			return power(pre_res, e_exp);
+			
+		} else
+			return e.to_polynomial(repl);
 
 	} else
 		return e;
@@ -2415,7 +2641,7 @@ term_done:	;
  *  'a*(b*x+b*y)' to 'a*b*(x+y)'. */
 ex collect_common_factors(const ex & e)
 {
-	if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
+	if (is_exactly_a<add>(e) || is_exactly_a<mul>(e) || is_exactly_a<power>(e)) {
 
 		exmap repl;
 		ex factor = 1;