X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Ffactor.cpp;h=dbe3ed3ea3bf2c494650cd13e3b0a190c6f77dc0;hp=f70e41bd37ffc5ebbdbc6618af445c14d6e25608;hb=4eff40dfafb72aeb484910c9ecef791b95332e9f;hpb=761d1597532504762c1f9b438c7727c4f74d7da3

diff --git a/ginac/factor.cpp b/ginac/factor.cpp
index f70e41bd..dbe3ed3e 100644
--- a/ginac/factor.cpp
+++ b/ginac/factor.cpp
@@ -33,7 +33,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2017 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2020 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
@@ -66,7 +66,6 @@
 #include "add.h"
 
 #include <algorithm>
-#include <cmath>
 #include <limits>
 #include <list>
 #include <vector>
@@ -1492,6 +1491,9 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
 	poly.unitcontprim(x, unit, cont, prim_ex);
 	upoly prim;
 	upoly_from_ex(prim, prim_ex, x);
+	if (prim_ex.is_equal(1)) {
+		return poly;
+	}
 
 	// determine proper prime and minimize number of modular factors
 	prime = 3;
@@ -1874,7 +1876,8 @@ static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex&
 {
 	vector<ex> a = a_;
 
-	const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
+	const cl_I modulus = expt_pos(cl_I(p),k);
+	const cl_modint_ring R = find_modint_ring(modulus);
 	const size_t r = a.size();
 	const size_t nu = I.size() + 1;
 
@@ -1947,16 +1950,13 @@ static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex&
 			cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
 			upvec delta_s = univar_diophant(amod, x, m, p, k);
 			cl_MI modcm;
-			cl_I poscm = cm;
-			while ( poscm < 0 ) {
-				poscm = poscm + expt_pos(cl_I(p),k);
-			}
+			cl_I poscm = plusp(cm) ? cm : mod(cm, modulus);
 			modcm = cl_MI(R, poscm);
 			for ( size_t j=0; j<delta_s.size(); ++j ) {
 				delta_s[j] = delta_s[j] * modcm;
 				sigma[j] = sigma[j] + umodpoly_to_ex(delta_s[j], x);
 			}
-			if ( nterms > 1 ) {
+			if ( nterms > 1 && i+1 != nterms ) {
 				z = c.op(i+1);
 			}
 		}
@@ -2070,13 +2070,8 @@ static ex hensel_multivar(const ex& a, const ex& x, const vector<EvalPoint>& I,
 		acand *= U[i];
 	}
 	if ( expand(a-acand).is_zero() ) {
-		lst res;
-		for ( size_t i=0; i<U.size(); ++i ) {
-			res.append(U[i]);
-		}
-		return res;
+		return lst(U.begin(), U.end());
 	} else {
-		lst res;
 		return lst{};
 	}
 }
@@ -2237,14 +2232,13 @@ static ex factor_sqrfree(const ex& poly);
  */
 static ex factor_multivariate(const ex& poly, const exset& syms)
 {
-	exset::const_iterator s;
 	const ex& x = *syms.begin();
 
 	// make polynomial primitive
 	ex unit, cont, pp;
 	poly.unitcontprim(x, unit, cont, pp);
 	if ( !is_a<numeric>(cont) ) {
-		return factor_sqrfree(cont) * factor_sqrfree(pp);
+		return unit * factor_sqrfree(cont) * factor_sqrfree(pp);
 	}
 
 	// factor leading coefficient
@@ -2317,7 +2311,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 			vector<numeric> ftilde(vnlst.nops()-1);
 			for ( size_t i=0; i<ftilde.size(); ++i ) {
 				ex ft = vnlst.op(i+1);
-				s = syms.begin();
+				auto s = syms.begin();
 				++s;
 				for ( size_t j=0; j<a.size(); ++j ) {
 					ft = ft.subs(*s == a[j]);
@@ -2391,7 +2385,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 		// set up evaluation points
 		EvalPoint ep;
 		vector<EvalPoint> epv;
-		s = syms.begin();
+		auto s = syms.begin();
 		++s;
 		for ( size_t i=0; i<a.size(); ++i ) {
 			ep.x = *s++;
@@ -2463,9 +2457,9 @@ static ex factor_sqrfree(const ex& poly)
 	if ( findsymbols.syms.size() == 1 ) {
 		// univariate case
 		const ex& x = *(findsymbols.syms.begin());
-		if ( poly.ldegree(x) > 0 ) {
+		int ld = poly.ldegree(x);
+		if ( ld > 0 ) {
 			// pull out direct factors
-			int ld = poly.ldegree(x);
 			ex res = factor_univariate(expand(poly/pow(x, ld)), x);
 			return res * pow(x,ld);
 		} else {
@@ -2505,13 +2499,41 @@ struct apply_factor_map : public map_function {
 	}
 };
 
-} // anonymous namespace
+/** Iterate through explicit factors of e, call yield(f, k) for
+ *  each factor of the form f^k.
+ *
+ *  Note that this function doesn't factor e itself, it only
+ *  iterates through the factors already explicitly present.
+ */
+template <typename F> void
+factor_iter(const ex &e, F yield)
+{
+	if (is_a<mul>(e)) {
+		for (const auto &f : e) {
+			if (is_a<power>(f)) {
+				yield(f.op(0), f.op(1));
+			} else {
+				yield(f, ex(1));
+			}
+		}
+	} else {
+		if (is_a<power>(e)) {
+			yield(e.op(0), e.op(1));
+		} else {
+			yield(e, ex(1));
+		}
+	}
+}
 
-/** Interface function to the outside world. It checks the arguments, tries a
- *  square free factorization, and then calls factor_sqrfree to do the hard
- *  work.
+/** This function factorizes a polynomial. It checks the arguments,
+ *  tries a square free factorization, and then calls factor_sqrfree
+ *  to do the hard work.
+ *
+ *  This function expands its argument, so for polynomials with
+ *  explicit factors it's better to call it on each one separately
+ *  (or use factor() which does just that).
  */
-ex factor(const ex& poly, unsigned options)
+static ex factor1(const ex& poly, unsigned options)
 {
 	// check arguments
 	if ( !poly.info(info_flags::polynomial) ) {
@@ -2538,48 +2560,35 @@ ex factor(const ex& poly, unsigned options)
 	ex sfpoly = sqrfree(poly.expand(), syms);
 
 	// factorize the square free components
-	if ( is_a<power>(sfpoly) ) {
-		// case: (polynomial)^exponent
-		const ex& base = sfpoly.op(0);
-		if ( !is_a<add>(base) ) {
-			// simple case: (monomial)^exponent
-			return sfpoly;
-		}
-		ex f = factor_sqrfree(base);
-		return pow(f, sfpoly.op(1));
-	}
-	if ( is_a<mul>(sfpoly) ) {
-		// case: multiple factors
-		ex res = 1;
-		for ( size_t i=0; i<sfpoly.nops(); ++i ) {
-			const ex& t = sfpoly.op(i);
-			if ( is_a<power>(t) ) {
-				const ex& base = t.op(0);
-				if ( !is_a<add>(base) ) {
-					res *= t;
-				} else {
-					ex f = factor_sqrfree(base);
-					res *= pow(f, t.op(1));
-				}
-			} else if ( is_a<add>(t) ) {
-				ex f = factor_sqrfree(t);
-				res *= f;
+	ex res = 1;
+	factor_iter(sfpoly,
+		[&](const ex &f, const ex &k) {
+			if ( is_a<add>(f) ) {
+				res *= pow(factor_sqrfree(f), k);
 			} else {
-				res *= t;
+				// simple case: (monomial)^exponent
+				res *= pow(f, k);
 			}
-		}
-		return res;
-	}
-	if ( is_a<symbol>(sfpoly) ) {
-		return poly;
-	}
-	// case: (polynomial)
-	ex f = factor_sqrfree(sfpoly);
-	return f;
+		});
+	return res;
 }
 
-} // namespace GiNaC
+} // anonymous namespace
 
-#ifdef DEBUGFACTOR
-#include "test.h"
-#endif
+/** Interface function to the outside world. It uses factor1()
+ *  on each of the explicitly present factors of poly.
+ */
+ex factor(const ex& poly, unsigned options)
+{
+	ex result = 1;
+	factor_iter(poly,
+		[&](const ex &f1, const ex &k1) {
+			factor_iter(factor1(f1, options),
+				[&](const ex &f2, const ex &k2) {
+					result *= pow(f2, k1*k2);
+				});
+		});
+	return result;
+}
+
+} // namespace GiNaC