X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Ffactor.cpp;h=8f8c87e858340e625df38a3bffe1fb3cfc0356fa;hp=da870d3daf079a2d202376177074dbbc897e84e2;hb=852fb174c43b49a3ce1b568f959d23e8a1959ee1;hpb=a12c1d88131c5a301d35767f0c4c947064703418

diff --git a/ginac/factor.cpp b/ginac/factor.cpp
index da870d3d..8f8c87e8 100644
--- a/ginac/factor.cpp
+++ b/ginac/factor.cpp
@@ -33,7 +33,7 @@
  */
 
 /*
- *  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
@@ -86,15 +86,16 @@ namespace GiNaC {
 #define DCOUT2(str,var) cout << #str << ": " << var << endl
 ostream& operator<<(ostream& o, const vector<int>& v)
 {
-	vector<int>::const_iterator i = v.begin(), end = v.end();
+	auto i = v.begin(), end = v.end();
 	while ( i != end ) {
-		o << *i++ << " ";
+		o << *i << " ";
+		++i;
 	}
 	return o;
 }
 static ostream& operator<<(ostream& o, const vector<cl_I>& v)
 {
-	vector<cl_I>::const_iterator i = v.begin(), end = v.end();
+	auto i = v.begin(), end = v.end();
 	while ( i != end ) {
 		o << *i << "[" << i-v.begin() << "]" << " ";
 		++i;
@@ -103,7 +104,7 @@ static ostream& operator<<(ostream& o, const vector<cl_I>& v)
 }
 static ostream& operator<<(ostream& o, const vector<cl_MI>& v)
 {
-	vector<cl_MI>::const_iterator i = v.begin(), end = v.end();
+	auto i = v.begin(), end = v.end();
 	while ( i != end ) {
 		o << *i << "[" << i-v.begin() << "]" << " ";
 		++i;
@@ -117,9 +118,9 @@ ostream& operator<<(ostream& o, const vector<numeric>& v)
 	}
 	return o;
 }
-ostream& operator<<(ostream& o, const vector< vector<cl_MI> >& v)
+ostream& operator<<(ostream& o, const vector<vector<cl_MI>>& v)
 {
-	vector< vector<cl_MI> >::const_iterator i = v.begin(), end = v.end();
+	auto i = v.begin(), end = v.end();
 	while ( i != end ) {
 		o << i-v.begin() << ": " << *i << endl;
 		++i;
@@ -218,8 +219,32 @@ static void expt_pos(umodpoly& a, unsigned int q)
 	}
 }
 
+template<bool COND, typename T = void> struct enable_if
+{
+	typedef T type;
+};
+
+template<typename T> struct enable_if<false, T> { /* empty */ };
+
+template<typename T> struct uvar_poly_p
+{
+	static const bool value = false;
+};
+
+template<> struct uvar_poly_p<upoly>
+{
+	static const bool value = true;
+};
+
+template<> struct uvar_poly_p<umodpoly>
+{
+	static const bool value = true;
+};
+
 template<typename T>
-static T operator+(const T& a, const T& b)
+// Don't define this for anything but univariate polynomials.
+static typename enable_if<uvar_poly_p<T>::value, T>::type
+operator+(const T& a, const T& b)
 {
 	int sa = a.size();
 	int sb = b.size();
@@ -250,7 +275,11 @@ static T operator+(const T& a, const T& b)
 }
 
 template<typename T>
-static T operator-(const T& a, const T& b)
+// Don't define this for anything but univariate polynomials. Otherwise
+// overload resolution might fail (this actually happens when compiling
+// GiNaC with g++ 3.4).
+static typename enable_if<uvar_poly_p<T>::value, T>::type
+operator-(const T& a, const T& b)
 {
 	int sa = a.size();
 	int sb = b.size();
@@ -470,11 +499,10 @@ static void reduce_coeff(umodpoly& a, const cl_I& x)
 	if ( a.empty() ) return;
 
 	cl_modint_ring R = a[0].ring();
-	umodpoly::iterator i = a.begin(), end = a.end();
-	for ( ; i!=end; ++i ) {
+	for (auto & i : a) {
 		// cln cannot perform this division in the modular field
-		cl_I c = R->retract(*i);
-		*i = cl_MI(R, the<cl_I>(c / x));
+		cl_I c = R->retract(i);
+		i = cl_MI(R, the<cl_I>(c / x));
 	}
 }
 
@@ -651,7 +679,9 @@ typedef vector<cl_MI> mvec;
 
 class modular_matrix
 {
+#ifdef DEBUGFACTOR
 	friend ostream& operator<<(ostream& o, const modular_matrix& m);
+#endif
 public:
 	modular_matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_)
 	{
@@ -663,90 +693,75 @@ public:
 	cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; }
 	void mul_col(size_t col, const cl_MI x)
 	{
-		mvec::iterator i = m.begin() + col;
 		for ( size_t rc=0; rc<r; ++rc ) {
-			*i = *i * x;
-			i += c;
+			std::size_t i = c*rc + col;
+			m[i] = m[i] * x;
 		}
 	}
 	void sub_col(size_t col1, size_t col2, const cl_MI fac)
 	{
-		mvec::iterator i1 = m.begin() + col1;
-		mvec::iterator i2 = m.begin() + col2;
 		for ( size_t rc=0; rc<r; ++rc ) {
-			*i1 = *i1 - *i2 * fac;
-			i1 += c;
-			i2 += c;
+			std::size_t i1 = col1 + c*rc;
+			std::size_t i2 = col2 + c*rc;
+			m[i1] = m[i1] - m[i2]*fac;
 		}
 	}
 	void switch_col(size_t col1, size_t col2)
 	{
-		cl_MI buf;
-		mvec::iterator i1 = m.begin() + col1;
-		mvec::iterator i2 = m.begin() + col2;
 		for ( size_t rc=0; rc<r; ++rc ) {
-			buf = *i1; *i1 = *i2; *i2 = buf;
-			i1 += c;
-			i2 += c;
+			std::size_t i1 = col1 + rc*c;
+			std::size_t i2 = col2 + rc*c;
+			std::swap(m[i1], m[i2]);
 		}
 	}
 	void mul_row(size_t row, const cl_MI x)
 	{
-		vector<cl_MI>::iterator i = m.begin() + row*c;
 		for ( size_t cc=0; cc<c; ++cc ) {
-			*i = *i * x;
-			++i;
+			std::size_t i = row*c + cc; 
+			m[i] = m[i] * x;
 		}
 	}
 	void sub_row(size_t row1, size_t row2, const cl_MI fac)
 	{
-		vector<cl_MI>::iterator i1 = m.begin() + row1*c;
-		vector<cl_MI>::iterator i2 = m.begin() + row2*c;
 		for ( size_t cc=0; cc<c; ++cc ) {
-			*i1 = *i1 - *i2 * fac;
-			++i1;
-			++i2;
+			std::size_t i1 = row1*c + cc;
+			std::size_t i2 = row2*c + cc;
+			m[i1] = m[i1] - m[i2]*fac;
 		}
 	}
 	void switch_row(size_t row1, size_t row2)
 	{
-		cl_MI buf;
-		vector<cl_MI>::iterator i1 = m.begin() + row1*c;
-		vector<cl_MI>::iterator i2 = m.begin() + row2*c;
 		for ( size_t cc=0; cc<c; ++cc ) {
-			buf = *i1; *i1 = *i2; *i2 = buf;
-			++i1;
-			++i2;
+			std::size_t i1 = row1*c + cc;
+			std::size_t i2 = row2*c + cc;
+			std::swap(m[i1], m[i2]);
 		}
 	}
 	bool is_col_zero(size_t col) const
 	{
-		mvec::const_iterator i = m.begin() + col;
 		for ( size_t rr=0; rr<r; ++rr ) {
-			if ( !zerop(*i) ) {
+			std::size_t i = col + rr*c;
+			if ( !zerop(m[i]) ) {
 				return false;
 			}
-			i += c;
 		}
 		return true;
 	}
 	bool is_row_zero(size_t row) const
 	{
-		mvec::const_iterator i = m.begin() + row*c;
 		for ( size_t cc=0; cc<c; ++cc ) {
-			if ( !zerop(*i) ) {
+			std::size_t i = row*c + cc;
+			if ( !zerop(m[i]) ) {
 				return false;
 			}
-			++i;
 		}
 		return true;
 	}
 	void set_row(size_t row, const vector<cl_MI>& newrow)
 	{
-		mvec::iterator i1 = m.begin() + row*c;
-		mvec::const_iterator i2 = newrow.begin(), end = newrow.end();
-		for ( ; i2 != end; ++i1, ++i2 ) {
-			*i1 = *i2;
+		for (std::size_t i2 = 0; i2 < newrow.size(); ++i2) {
+			std::size_t i1 = row*c + i2;
+			m[i1] = newrow[i2];
 		}
 	}
 	mvec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; }
@@ -927,15 +942,13 @@ static void berlekamp(const umodpoly& a, upvec& upv)
 				}
 				factors.push_back(g);
 				size = 0;
-				list<umodpoly>::const_iterator i = factors.begin(), end = factors.end();
-				while ( i != end ) {
-					if ( degree(*i) ) ++size; 
-					++i;
+				for (auto & i : factors) {
+					if (degree(i))
+						++size;
 				}
 				if ( size == k ) {
-					list<umodpoly>::const_iterator i = factors.begin(), end = factors.end();
-					while ( i != end ) {
-						upv.push_back(*i++);
+					for (auto & i : factors) {
+						upv.push_back(i);
 					}
 					return;
 				}
@@ -1159,15 +1172,13 @@ static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& s, umodpol
 		d2 = r2;
 	}
 	cl_MI fac = recip(lcoeff(a) * lcoeff(c));
-	umodpoly::iterator i = s.begin(), end = s.end();
-	for ( ; i!=end; ++i ) {
-		*i = *i * fac;
+	for (auto & i : s) {
+		i = i * fac;
 	}
 	canonicalize(s);
 	fac = recip(lcoeff(b) * lcoeff(c));
-	i = t.begin(), end = t.end();
-	for ( ; i!=end; ++i ) {
-		*i = *i * fac;
+	for (auto & i : t) {
+		i = i * fac;
 	}
 	canonicalize(t);
 }
@@ -1319,7 +1330,6 @@ static unsigned int next_prime(unsigned int p)
 	if ( primes.size() == 0 ) {
 		primes.push_back(3); primes.push_back(5); primes.push_back(7);
 	}
-	vector<unsigned int>::const_iterator it = primes.begin();
 	if ( p >= primes.back() ) {
 		unsigned int candidate = primes.back() + 2;
 		while ( true ) {
@@ -1334,10 +1344,9 @@ static unsigned int next_prime(unsigned int p)
 		}
 		return candidate;
 	}
-	vector<unsigned int>::const_iterator end = primes.end();
-	for ( ; it!=end; ++it ) {
-		if ( *it > p ) {
-			return *it;
+	for (auto & it : primes) {
+		if ( it > p ) {
+			return it;
 		}
 	}
 	throw logic_error("next_prime: should not reach this point!");
@@ -1455,7 +1464,7 @@ private:
 	}
 private:
 	umodpoly lr[2];
-	vector< vector<umodpoly> > cache;
+	vector<vector<umodpoly>> cache;
 	upvec factors;
 	umodpoly one;
 	size_t n;
@@ -1496,7 +1505,17 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
 	cl_modint_ring R;
 	unsigned int trials = 0;
 	unsigned int minfactors = 0;
-	cl_I lc = lcoeff(prim) * the<cl_I>(ex_to<numeric>(cont).to_cl_N());
+
+	const numeric& cont_n = ex_to<numeric>(cont);
+	cl_I i_cont;
+	if (cont_n.is_integer()) {
+		i_cont = the<cl_I>(cont_n.to_cl_N());
+	} else {
+		// poly \in Q[x] => poly = q ipoly, ipoly \in Z[x], q \in Q
+		// factor(poly) \equiv q factor(ipoly)
+		i_cont = cl_I(1);
+	}
+	cl_I lc = lcoeff(prim)*i_cont;
 	upvec factors;
 	while ( trials < 2 ) {
 		umodpoly modpoly;
@@ -1580,7 +1599,7 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
 				}
 				else {
 					upvec newfactors1(part.size_left()), newfactors2(part.size_right());
-					upvec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
+					auto i1 = newfactors1.begin(), i2 = newfactors2.begin();
 					for ( size_t i=0; i<n; ++i ) {
 						if ( part[i] ) {
 							*i2++ = tocheck.top().factors[i];
@@ -1696,9 +1715,8 @@ static void change_modulus(const cl_modint_ring& R, umodpoly& a)
 {
 	if ( a.empty() ) return;
 	cl_modint_ring oldR = a[0].ring();
-	umodpoly::iterator i = a.begin(), end = a.end();
-	for ( ; i!=end; ++i ) {
-		*i = R->canonhom(oldR->retract(*i));
+	for (auto & i : a) {
+		i = R->canonhom(oldR->retract(i));
 	}
 	canonicalize(a);
 }
@@ -2095,8 +2113,9 @@ static ex put_factors_into_lst(const ex& e)
 		return result;
 	}
 	if ( is_a<symbol>(e) || is_a<add>(e) ) {
-		result.append(1);
-		result.append(e);
+		ex icont(e.integer_content());
+		result.append(icont);
+		result.append(e/icont);
 		return result;
 	}
 	if ( is_a<mul>(e) ) {
@@ -2532,9 +2551,8 @@ ex factor(const ex& poly, unsigned options)
 		return poly;
 	}
 	lst syms;
-	exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
-	for ( ; i!=end; ++i ) {
-		syms.append(*i);
+	for (auto & i : findsymbols.syms ) {
+		syms.append(i);
 	}
 
 	// make poly square free