X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Ffactor.cpp;h=2b77d1b0110f3e327216a11a20510b901a6b93ee;hb=709e61093f43462b5816f859a32a31cc00758da8;hp=fb73897218db06ed31b40f32bf52aabca7e33acf;hpb=2b0ad5c381dc081cc4066b0c5f939f5169ad9ab3;p=ginac.git

diff --git a/ginac/factor.cpp b/ginac/factor.cpp
index fb738972..2b77d1b0 100644
--- a/ginac/factor.cpp
+++ b/ginac/factor.cpp
@@ -1,16 +1,39 @@
 /** @file factor.cpp
  *
- *  Polynomial factorization code (implementation).
+ *  Polynomial factorization (implementation).
+ *
+ *  The interface function factor() at the end of this file is defined in the
+ *  GiNaC namespace. All other utility functions and classes are defined in an
+ *  additional anonymous namespace.
+ *
+ *  Factorization starts by doing a square free factorization and making the
+ *  coefficients integer. Then, depending on the number of free variables it
+ *  proceeds either in dedicated univariate or multivariate factorization code.
+ *
+ *  Univariate factorization does a modular factorization via Berlekamp's
+ *  algorithm and distinct degree factorization. Hensel lifting is used at the
+ *  end.
+ *  
+ *  Multivariate factorization uses the univariate factorization (applying a
+ *  evaluation homomorphism first) and Hensel lifting raises the answer to the
+ *  multivariate domain. The Hensel lifting code is completely distinct from the
+ *  code used by the univariate factorization.
  *
  *  Algorithms used can be found in
- *    [W1]  An Improved Multivariate Polynomial Factoring Algorithm,
- *          P.S.Wang, Mathematics of Computation, Vol. 32, No. 144 (1978) 1215--1231.
+ *    [Wan] An Improved Multivariate Polynomial Factoring Algorithm,
+ *          P.S.Wang,
+ *          Mathematics of Computation, Vol. 32, No. 144 (1978) 1215--1231.
  *    [GCL] Algorithms for Computer Algebra,
- *          K.O.Geddes, S.R.Czapor, G.Labahn, Springer Verlag, 1992.
+ *          K.O.Geddes, S.R.Czapor, G.Labahn,
+ *          Springer Verlag, 1992.
+ *    [Mig] Some Useful Bounds,
+ *          M.Mignotte, 
+ *          In "Computer Algebra, Symbolic and Algebraic Computation" (B.Buchberger et al., eds.),
+ *          pp. 259-263, Springer-Verlag, New York, 1982.
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2019 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
@@ -43,10 +66,10 @@
 #include "add.h"
 
 #include <algorithm>
-#include <cmath>
 #include <limits>
 #include <list>
 #include <vector>
+#include <stack>
 #ifdef DEBUGFACTOR
 #include <ostream>
 #endif
@@ -61,53 +84,57 @@ namespace GiNaC {
 #define DCOUT(str) cout << #str << endl
 #define DCOUTVAR(var) cout << #var << ": " << var << endl
 #define DCOUT2(str,var) cout << #str << ": " << var << endl
-#else
-#define DCOUT(str)
-#define DCOUTVAR(var)
-#define DCOUT2(str,var)
-#endif
-
-// anonymous namespace to hide all utility functions
-namespace {
-
-typedef vector<cl_MI> mvec;
-#ifdef DEBUGFACTOR
 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;
 }
-ostream& operator<<(ostream& o, const vector<cl_I>& v)
+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;
 	}
 	return o;
 }
-ostream& operator<<(ostream& o, const vector<cl_MI>& 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;
 	}
 	return o;
 }
-ostream& operator<<(ostream& o, const vector< vector<cl_MI> >& v)
+ostream& operator<<(ostream& o, const vector<numeric>& v)
+{
+	for ( size_t i=0; i<v.size(); ++i ) {
+		o << v[i] << " ";
+	}
+	return o;
+}
+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;
 	}
 	return o;
 }
-#endif
+#else
+#define DCOUT(str)
+#define DCOUTVAR(var)
+#define DCOUT2(str,var)
+#endif // def DEBUGFACTOR
+
+// anonymous namespace to hide all utility functions
+namespace {
 
 ////////////////////////////////////////////////////////////////////////////////
 // modular univariate polynomial code
@@ -192,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();
@@ -224,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();
@@ -444,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));
 	}
 }
 
@@ -621,9 +675,13 @@ static bool squarefree(const umodpoly& a)
 ////////////////////////////////////////////////////////////////////////////////
 // modular matrix
 
+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_)
 	{
@@ -635,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; }
@@ -770,6 +813,11 @@ ostream& operator<<(ostream& o, const modular_matrix& m)
 // END modular matrix
 ////////////////////////////////////////////////////////////////////////////////
 
+/** Calculates the Q matrix for a polynomial. Used by Berlekamp's algorithm.
+ *
+ *  @param[in]  a_  modular polynomial
+ *  @param[out] Q   Q matrix
+ */
 static void q_matrix(const umodpoly& a_, modular_matrix& Q)
 {
 	umodpoly a = a_;
@@ -793,6 +841,11 @@ static void q_matrix(const umodpoly& a_, modular_matrix& Q)
 	}
 }
 
+/** Determine the nullspace of a matrix M-1.
+ *
+ *  @param[in,out] M      matrix, will be modified
+ *  @param[out]    basis  calculated nullspace of M-1
+ */
 static void nullspace(modular_matrix& M, vector<mvec>& basis)
 {
 	const size_t n = M.rowsize();
@@ -837,11 +890,20 @@ static void nullspace(modular_matrix& M, vector<mvec>& basis)
 	}
 }
 
+/** Berlekamp's modular factorization.
+ *  
+ *  The implementation follows the algorithm in chapter 8 of [GCL].
+ *
+ *  @param[in]  a    modular polynomial
+ *  @param[out] upv  vector containing modular factors. if upv was not empty the
+ *                   new elements are added at the end
+ */
 static void berlekamp(const umodpoly& a, upvec& upv)
 {
 	cl_modint_ring R = a[0].ring();
 	umodpoly one(1, R->one());
 
+	// find nullspace of Q matrix
 	modular_matrix Q(degree(a), degree(a), R->zero());
 	q_matrix(a, Q);
 	vector<mvec> nu;
@@ -849,17 +911,18 @@ static void berlekamp(const umodpoly& a, upvec& upv)
 
 	const unsigned int k = nu.size();
 	if ( k == 1 ) {
+		// irreducible
 		return;
 	}
 
-	list<umodpoly> factors;
-	factors.push_back(a);
+	list<umodpoly> factors = {a};
 	unsigned int size = 1;
 	unsigned int r = 1;
 	unsigned int q = cl_I_to_uint(R->modulus);
 
 	list<umodpoly>::iterator u = factors.begin();
 
+	// calculate all gcd's
 	while ( true ) {
 		for ( unsigned int s=0; s<q; ++s ) {
 			umodpoly nur = nu[r];
@@ -872,21 +935,18 @@ static void berlekamp(const umodpoly& a, upvec& upv)
 				div(*u, g, uo);
 				if ( equal_one(uo) ) {
 					throw logic_error("berlekamp: unexpected divisor.");
-				}
-				else {
+				} else {
 					*u = uo;
 				}
 				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;
 				}
@@ -899,6 +959,16 @@ static void berlekamp(const umodpoly& a, upvec& upv)
 	}
 }
 
+// modular square free factorization is not used at the moment so we deactivate
+// the code
+#if 0
+
+/** Calculates a^(1/prime).
+ *  
+ *  @param[in] a      polynomial
+ *  @param[in] prime  prime number -> exponent 1/prime
+ *  @param[in] ap     resulting polynomial
+ */
 static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
 {
 	size_t newdeg = degree(a)/prime;
@@ -909,6 +979,12 @@ static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
 	}
 }
 
+/** Modular square free factorization.
+ *
+ *  @param[in]  a        polynomial
+ *  @param[out] factors  modular factors
+ *  @param[out] mult     corresponding multiplicities (exponents)
+ */
 static void modsqrfree(const umodpoly& a, upvec& factors, vector<int>& mult)
 {
 	const unsigned int prime = cl_I_to_uint(a[0].ring()->modulus);
@@ -942,8 +1018,7 @@ static void modsqrfree(const umodpoly& a, upvec& factors, vector<int>& mult)
 				mult[i] *= prime;
 			}
 		}
-	}
-	else {
+	} else {
 		umodpoly ap;
 		expt_1_over_p(a, prime, ap);
 		size_t previ = mult.size();
@@ -954,6 +1029,18 @@ static void modsqrfree(const umodpoly& a, upvec& factors, vector<int>& mult)
 	}
 }
 
+#endif // deactivation of square free factorization
+
+/** Distinct degree factorization (DDF).
+ *  
+ *  The implementation follows the algorithm in chapter 8 of [GCL].
+ *
+ *  @param[in]  a_         modular polynomial
+ *  @param[out] degrees    vector containing the degrees of the factors of the
+ *                         corresponding polynomials in ddfactors.
+ *  @param[out] ddfactors  vector containing polynomials which factors have the
+ *                         degree given in degrees.
+ */
 static void distinct_degree_factor(const umodpoly& a_, vector<int>& degrees, upvec& ddfactors)
 {
 	umodpoly a = a_;
@@ -995,6 +1082,16 @@ static void distinct_degree_factor(const umodpoly& a_, vector<int>& degrees, upv
 	}
 }
 
+/** Modular same degree factorization.
+ *  Same degree factorization is a kind of misnomer. It performs distinct degree
+ *  factorization, but instead of using the Cantor-Zassenhaus algorithm it
+ *  (sub-optimally) uses Berlekamp's algorithm for the factors of the same
+ *  degree.
+ *
+ *  @param[in]  a    modular polynomial
+ *  @param[out] upv  vector containing modular factors. if upv was not empty the
+ *                   new elements are added at the end
+ */
 static void same_degree_factor(const umodpoly& a, upvec& upv)
 {
 	cl_modint_ring R = a[0].ring();
@@ -1006,15 +1103,25 @@ static void same_degree_factor(const umodpoly& a, upvec& upv)
 	for ( size_t i=0; i<degrees.size(); ++i ) {
 		if ( degrees[i] == degree(ddfactors[i]) ) {
 			upv.push_back(ddfactors[i]);
-		}
-		else {
+		} else {
 			berlekamp(ddfactors[i], upv);
 		}
 	}
 }
 
+// Yes, we can (choose).
 #define USE_SAME_DEGREE_FACTOR
 
+/** Modular univariate factorization.
+ *
+ *  In principle, we have two algorithms at our disposal: Berlekamp's algorithm
+ *  and same degree factorization (SDF). SDF seems to be slightly faster in
+ *  almost all cases so it is activated as default.
+ *
+ *  @param[in]  p    modular polynomial
+ *  @param[out] upv  vector containing modular factors. if upv was not empty the
+ *                   new elements are added at the end
+ */
 static void factor_modular(const umodpoly& p, upvec& upv)
 {
 #ifdef USE_SAME_DEGREE_FACTOR
@@ -1024,7 +1131,7 @@ static void factor_modular(const umodpoly& p, upvec& upv)
 #endif
 }
 
-/** Calculates polynomials s and t such that a*s+b*t==1.
+/** Calculates modular polynomials s and t such that a*s+b*t==1.
  *  Assertion: a and b are relatively prime and not zero.
  *
  *  @param[in]  a  polynomial
@@ -1061,19 +1168,23 @@ 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);
 }
 
+/** Replaces the leading coefficient in a polynomial by a given number.
+ *
+ *  @param[in] poly  polynomial to change
+ *  @param[in] lc    new leading coefficient
+ *  @return          changed polynomial
+ */
 static upoly replace_lc(const upoly& poly, const cl_I& lc)
 {
 	if ( poly.empty() ) return poly;
@@ -1082,6 +1193,9 @@ static upoly replace_lc(const upoly& poly, const cl_I& lc)
 	return r;
 }
 
+/** Calculates the bound for the modulus.
+ *  See [Mig].
+ */
 static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg)
 {
 	cl_I maxcoeff = 0;
@@ -1096,6 +1210,9 @@ static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg)
 	return ( B > maxcoeff ) ? B : maxcoeff;
 }
 
+/** Calculates the bound for the modulus.
+ *  See [Mig].
+ */
 static inline cl_I calc_bound(const upoly& a, int maxdeg)
 {
 	cl_I maxcoeff = 0;
@@ -1110,6 +1227,18 @@ static inline cl_I calc_bound(const upoly& a, int maxdeg)
 	return ( B > maxcoeff ) ? B : maxcoeff;
 }
 
+/** Hensel lifting as used by factor_univariate().
+ *
+ *  The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ *  @param[in]  a_   primitive univariate polynomials
+ *  @param[in]  p    prime number that does not divide lcoeff(a)
+ *  @param[in]  u1_  modular factor of a (mod p)
+ *  @param[in]  w1_  modular factor of a (mod p), relatively prime to u1_,
+ *                   fulfilling  u1_*w1_ == a mod p
+ *  @param[out] u    lifted factor
+ *  @param[out] w    lifted factor, u*w = a
+ */
 static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, upoly& u, upoly& w)
 {
 	upoly a = a_;
@@ -1180,45 +1309,52 @@ static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_,
 		if ( alpha != 1 ) {
 			w = w / alpha;
 		}
-	}
-	else {
+	} else {
 		u.clear();
 	}
 }
 
+/** Returns a new prime number.
+ *
+ *  @param[in] p  prime number
+ *  @return       next prime number after p
+ */
 static unsigned int next_prime(unsigned int p)
 {
 	static vector<unsigned int> primes;
-	if ( primes.size() == 0 ) {
-		primes.push_back(3); primes.push_back(5); primes.push_back(7);
+	if (primes.empty()) {
+		primes = {3, 5, 7};
 	}
-	vector<unsigned int>::const_iterator it = primes.begin();
 	if ( p >= primes.back() ) {
 		unsigned int candidate = primes.back() + 2;
 		while ( true ) {
 			size_t n = primes.size()/2;
 			for ( size_t i=0; i<n; ++i ) {
-				if ( candidate % primes[i] ) continue;
+				if (candidate % primes[i])
+					continue;
 				candidate += 2;
 				i=-1;
 			}
 			primes.push_back(candidate);
-			if ( candidate > p ) break;
+			if (candidate > p)
+				break;
 		}
 		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!");
 }
 
+/** Manages the splitting a vector of of modular factors into two partitions.
+ */
 class factor_partition
 {
 public:
+	/** Takes the vector of modular factors and initializes the first partition */
 	factor_partition(const upvec& factors_) : factors(factors_)
 	{
 		n = factors.size();
@@ -1234,9 +1370,8 @@ public:
 	size_t size() const { return n; }
 	size_t size_left() const { return n-len; }
 	size_t size_right() const { return len; }
-#ifdef DEBUGFACTOR
-	void get() const { DCOUTVAR(k); }
-#endif
+	/** Initializes the next partition.
+	    Returns true, if there is one, false otherwise. */
 	bool next()
 	{
 		if ( last == n-1 ) {
@@ -1266,15 +1401,16 @@ public:
 			if ( len > n/2 ) return false;
 			fill(k.begin(), k.begin()+len, 1);
 			fill(k.begin()+len+1, k.end(), 0);
-		}
-		else {
+		} else {
 			k[last++] = 0;
 			k[last] = 1;
 		}
 		split();
 		return true;
 	}
+	/** Get first partition */
 	umodpoly& left() { return lr[0]; }
+	/** Get second partition */
 	umodpoly& right() { return lr[1]; }
 private:
 	void split_cached()
@@ -1288,8 +1424,7 @@ private:
 			if ( d ) {
 				if ( cache[pos].size() >= d ) {
 					lr[group] = lr[group] * cache[pos][d-1];
-				}
-				else {
+				} else {
 					if ( cache[pos].size() == 0 ) {
 						cache[pos].push_back(factors[pos] * factors[pos+1]);
 					}
@@ -1303,8 +1438,7 @@ private:
 					}
 					lr[group] = lr[group] * cache[pos].back();
 				}
-			}
-			else {
+			} else {
 				lr[group] = lr[group] * factors[pos];
 			}
 		} while ( i < n );
@@ -1315,8 +1449,7 @@ private:
 		lr[1] = one;
 		if ( n > 6 ) {
 			split_cached();
-		}
-		else {
+		} else {
 			for ( size_t i=0; i<n; ++i ) {
 				lr[k[i]] = lr[k[i]] * factors[i];
 			}
@@ -1324,7 +1457,7 @@ private:
 	}
 private:
 	umodpoly lr[2];
-	vector< vector<umodpoly> > cache;
+	vector<vector<umodpoly>> cache;
 	upvec factors;
 	umodpoly one;
 	size_t n;
@@ -1333,12 +1466,25 @@ private:
 	vector<int> k;
 };
 
+/** Contains a pair of univariate polynomial and its modular factors.
+ *  Used by factor_univariate().
+ */
 struct ModFactors
 {
 	upoly poly;
 	upvec factors;
 };
 
+/** Univariate polynomial factorization.
+ *
+ *  Modular factorization is tried for several primes to minimize the number of
+ *  modular factors. Then, Hensel lifting is performed.
+ *
+ *  @param[in]     poly   expanded square free univariate polynomial
+ *  @param[in]     x      symbol
+ *  @param[in,out] prime  prime number to start trying modular factorization with,
+ *                        output value is the prime number actually used
+ */
 static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
 {
 	ex unit, cont, prim_ex;
@@ -1352,7 +1498,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;
@@ -1378,8 +1534,7 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
 			minfactors = trialfactors.size();
 			lastp = prime;
 			trials = 1;
-		}
-		else {
+		} else {
 			++trials;
 		}
 	}
@@ -1398,8 +1553,10 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
 		const size_t n = tocheck.top().factors.size();
 		factor_partition part(tocheck.top().factors);
 		while ( true ) {
+			// call Hensel lifting
 			hensel_univar(tocheck.top().poly, prime, part.left(), part.right(), f1, f2);
 			if ( !f1.empty() ) {
+				// successful, update the stack and the result
 				if ( part.size_left() == 1 ) {
 					if ( part.size_right() == 1 ) {
 						result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
@@ -1431,15 +1588,13 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
 						}
 					}
 					break;
-				}
-				else {
+				} 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];
-						}
-						else {
+						} else {
 							*i1++ = tocheck.top().factors[i];
 						}
 					}
@@ -1451,9 +1606,11 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
 					tocheck.push(mf);
 					break;
 				}
-			}
-			else {
+			} else {
+				// not successful
 				if ( !part.next() ) {
+					// if no more combinations left, return polynomial as
+					// irreducible
 					result *= upoly_to_ex(tocheck.top().poly, x);
 					tocheck.pop();
 					break;
@@ -1465,21 +1622,51 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
 	return unit * cont * result;
 }
 
+/** Second interface to factor_univariate() to be used if the information about
+ *  the prime is not needed.
+ */
 static inline ex factor_univariate(const ex& poly, const ex& x)
 {
 	unsigned int prime;
 	return factor_univariate(poly, x, prime);
 }
 
+/** Represents an evaluation point (<symbol>==<integer>).
+ */
 struct EvalPoint
 {
 	ex x;
 	int evalpoint;
 };
 
+#ifdef DEBUGFACTOR
+ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
+{
+	for ( size_t i=0; i<v.size(); ++i ) {
+		o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
+	}
+	return o;
+}
+#endif // def DEBUGFACTOR
+
 // forward declaration
 static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex& c, const vector<EvalPoint>& I, unsigned int d, unsigned int p, unsigned int k);
 
+/** Utility function for multivariate Hensel lifting.
+ *
+ *  Solves the equation
+ *    s_1*b_1 + ... + s_r*b_r == 1 mod p^k
+ *  with deg(s_i) < deg(a_i)
+ *  and with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r
+ *
+ *  The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ *  @param[in]  a   vector of modular univariate polynomials
+ *  @param[in]  x   symbol
+ *  @param[in]  p   prime number
+ *  @param[in]  k   p^k is modulus
+ *  @return         vector of polynomials (s_i)
+ */
 static upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned int k)
 {
 	const size_t r = a.size();
@@ -1508,20 +1695,35 @@ static upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, uns
 	return s;
 }
 
-/**
- *  Assert: a not empty.
+/** Changes the modulus of a modular polynomial. Used by eea_lift().
+ *
+ *  @param[in]     R  new modular ring
+ *  @param[in,out] a  polynomial to change (in situ)
  */
 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);
 }
 
+/** Utility function for multivariate Hensel lifting.
+ *
+ *  Solves  s*a + t*b == 1 mod p^k  given a,b.
+ *
+ *  The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ *  @param[in]  a   polynomial
+ *  @param[in]  b   polynomial
+ *  @param[in]  x   symbol
+ *  @param[in]  p   prime number
+ *  @param[in]  k   p^k is modulus
+ *  @param[out] s_  output polynomial
+ *  @param[out] t_  output polynomial
+ */
 static void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, unsigned int k, umodpoly& s_, umodpoly& t_)
 {
 	cl_modint_ring R = find_modint_ring(p);
@@ -1565,6 +1767,21 @@ static void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned
 	s_ = s; t_ = t;
 }
 
+/** Utility function for multivariate Hensel lifting.
+ *
+ *  Solves the equation
+ *    s_1*b_1 + ... + s_r*b_r == x^m mod p^k
+ *  with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r
+ *
+ *  The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ *  @param a  vector with univariate polynomials mod p^k
+ *  @param x  symbol
+ *  @param m  exponent of x^m in the equation to solve
+ *  @param p  prime number
+ *  @param k  p^k is modulus
+ *  @return   vector of polynomials (s_i)
+ */
 static upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
 {
 	cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
@@ -1579,8 +1796,7 @@ static upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsign
 			rem(bmod, a[j], buf);
 			result.push_back(buf);
 		}
-	}
-	else {
+	} else {
 		umodpoly s, t;
 		eea_lift(a[1], a[0], x, p, k, s, t);
 		umodpoly bmod = umodpoly_to_umodpoly(s, R, m);
@@ -1595,10 +1811,14 @@ static upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsign
 	return result;
 }
 
+/** Map used by function make_modular().
+ *  Finds every coefficient in a polynomial and replaces it by is value in the
+ *  given modular ring R (symmetric representation).
+ */
 struct make_modular_map : public map_function {
 	cl_modint_ring R;
 	make_modular_map(const cl_modint_ring& R_) : R(R_) { }
-	ex operator()(const ex& e)
+	ex operator()(const ex& e) override
 	{
 		if ( is_a<add>(e) || is_a<mul>(e) ) {
 			return e.map(*this);
@@ -1610,8 +1830,7 @@ struct make_modular_map : public map_function {
 			numeric n(R->retract(emod));
 			if ( n > halfmod ) {
 				return n-mod;
-			}
-			else {
+			} else {
 				return n;
 			}
 		}
@@ -1619,18 +1838,43 @@ struct make_modular_map : public map_function {
 	}
 };
 
+/** Helps mimicking modular multivariate polynomial arithmetic.
+ *
+ *  @param e  expression of which to make the coefficients equal to their value
+ *            in the modular ring R (symmetric representation)
+ *  @param R  modular ring
+ *  @return   resulting expression
+ */
 static ex make_modular(const ex& e, const cl_modint_ring& R)
 {
 	make_modular_map map(R);
 	return map(e.expand());
 }
 
+/** Utility function for multivariate Hensel lifting.
+ *
+ *  Returns the polynomials s_i that fulfill
+ *    s_1*b_1 + ... + s_r*b_r == c mod <I^(d+1),p^k>
+ *  with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r
+ *
+ *  The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ *  @param a_  vector of multivariate factors mod p^k
+ *  @param x   symbol (equiv. x_1 in [GCL])
+ *  @param c   polynomial mod p^k
+ *  @param I   vector of evaluation points
+ *  @param d   maximum total degree of result
+ *  @param p   prime number
+ *  @param k   p^k is modulus
+ *  @return    vector of polynomials (s_i)
+ */
 static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex& c, const vector<EvalPoint>& I,
                                     unsigned int d, unsigned int p, unsigned int k)
 {
 	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;
 
@@ -1680,8 +1924,7 @@ static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex&
 				e = make_modular(buf, R);
 			}
 		}
-	}
-	else {
+	} else {
 		upvec amod;
 		for ( size_t i=0; i<a.size(); ++i ) {
 			umodpoly up;
@@ -1695,8 +1938,7 @@ static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex&
 		if ( is_a<add>(c) ) {
 			nterms = c.nops();
 			z = c.op(0);
-		}
-		else {
+		} else {
 			nterms = 1;
 			z = c;
 		}
@@ -1705,16 +1947,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);
 			}
 		}
@@ -1727,17 +1966,24 @@ static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex&
 	return sigma;
 }
 
-#ifdef DEBUGFACTOR
-ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
-{
-	for ( size_t i=0; i<v.size(); ++i ) {
-		o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
-	}
-	return o;
-}
-#endif // def DEBUGFACTOR
-
-static ex hensel_multivar(const ex& a, const ex& x, const vector<EvalPoint>& I, unsigned int p, const cl_I& l, const upvec& u, const vector<ex>& lcU)
+/** Multivariate Hensel lifting.
+ *  The implementation follows the algorithm in chapter 6 of [GCL].
+ *  Since we don't have a data type for modular multivariate polynomials, the
+ *  respective operations are done in a GiNaC::ex and the function
+ *  make_modular() is then called to make the coefficient modular p^l.
+ *
+ *  @param a    multivariate polynomial primitive in x
+ *  @param x    symbol (equiv. x_1 in [GCL])
+ *  @param I    vector of evaluation points (x_2==a_2,x_3==a_3,...)
+ *  @param p    prime number (should not divide lcoeff(a mod I))
+ *  @param l    p^l is the modulus of the lifted univariate field
+ *  @param u    vector of modular (mod p^l) factors of a mod I
+ *  @param lcU  correct leading coefficient of the univariate factors of a mod I
+ *  @return     list GiNaC::lst with lifted factors (multivariate factors of a),
+ *              empty if Hensel lifting did not succeed
+ */
+static ex hensel_multivar(const ex& a, const ex& x, const vector<EvalPoint>& I,
+                          unsigned int p, const cl_I& l, const upvec& u, const vector<ex>& lcU)
 {
 	const size_t nu = I.size() + 1;
 	const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
@@ -1821,22 +2067,19 @@ 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;
-	}
-	else {
-		lst res;
-		return lst();
+		return lst(U.begin(), U.end());
+	} else {
+		return lst{};
 	}
 }
 
+/** Takes a factorized expression and puts the factors in a lst. The exponents
+ *  of the factors are discarded, e.g. 7*x^2*(y+1)^4 --> {7,x,y+1}. The first
+ *  element of the list is always the numeric coefficient.
+ */
 static ex put_factors_into_lst(const ex& e)
 {
 	lst result;
-
 	if ( is_a<numeric>(e) ) {
 		result.append(e);
 		return result;
@@ -1847,8 +2090,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) ) {
@@ -1871,20 +2115,11 @@ static ex put_factors_into_lst(const ex& e)
 	throw runtime_error("put_factors_into_lst: bad term.");
 }
 
-#ifdef DEBUGFACTOR
-ostream& operator<<(ostream& o, const vector<numeric>& v)
-{
-	for ( size_t i=0; i<v.size(); ++i ) {
-		o << v[i] << " ";
-	}
-	return o;
-}
-#endif // def DEBUGFACTOR
-
-/** Checks whether in a set of numbers each has a unique prime factor.
+/** Checks a set of numbers for whether each number has a unique prime factor.
  *
  *  @param[in]  f  list of numbers to check
- *  @return        true: if number set is bad, false: otherwise
+ *  @return        true: if number set is bad, false: if set is okay (has unique
+ *                 prime factors)
  */
 static bool checkdivisors(const lst& f)
 {
@@ -1914,7 +2149,7 @@ static bool checkdivisors(const lst& f)
  *  1. lcoeff(evaluated_polynomial) does not vanish
  *  2. factors of lcoeff(evaluated_polynomial) have each a unique prime factor
  *  3. evaluated_polynomial is square free
- *  See [W1] for more details.
+ *  See [Wan] for more details.
  *
  *  @param[in]     u        multivariate polynomial to be factored
  *  @param[in]     vn       leading coefficient of u in x (x==first symbol in syms)
@@ -1931,7 +2166,7 @@ static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst
 	const ex& x = *syms.begin();
 	while ( true ) {
 		++modulus;
-		/* generate a set of integers ... */
+		// generate a set of integers ...
 		u0 = u;
 		ex vna = vn;
 		ex vnatry;
@@ -1941,19 +2176,19 @@ static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst
 			do {
 				a[i] = mod(numeric(rand()), 2*modulus) - modulus;
 				vnatry = vna.subs(*s == a[i]);
-				/* ... for which the leading coefficient doesn't vanish ... */
+				// ... for which the leading coefficient doesn't vanish ...
 			} while ( vnatry == 0 );
 			vna = vnatry;
 			u0 = u0.subs(*s == a[i]);
 			++s;
 		}
-		/* ... for which u0 is square free ... */
+		// ... for which u0 is square free ...
 		ex g = gcd(u0, u0.diff(ex_to<symbol>(x)));
 		if ( !is_a<numeric>(g) ) {
 			continue;
 		}
 		if ( !is_a<numeric>(vn) ) {
-			/* ... and for which the evaluated factors have each an unique prime factor */
+			// ... and for which the evaluated factors have each an unique prime factor
 			lst fnum = f;
 			fnum.let_op(0) = fnum.op(0) * u0.content(x);
 			for ( size_t i=1; i<fnum.nops(); ++i ) {
@@ -1969,7 +2204,7 @@ static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst
 				continue;
 			}
 		}
-		/* ok, we have a valid set now */
+		// ok, we have a valid set now
 		return;
 	}
 }
@@ -1977,31 +2212,38 @@ static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst
 // forward declaration
 static ex factor_sqrfree(const ex& poly);
 
-/**
- *  ASSERT: poly is expanded
+/** Multivariate factorization.
+ *  
+ *  The implementation is based on the algorithm described in [Wan].
+ *  An evaluation homomorphism (a set of integers) is determined that fulfills
+ *  certain criteria. The evaluated polynomial is univariate and is factorized
+ *  by factor_univariate(). The main work then is to find the correct leading
+ *  coefficients of the univariate factors. They have to correspond to the
+ *  factors of the (multivariate) leading coefficient of the input polynomial
+ *  (as defined for a specific variable x). After that the Hensel lifting can be
+ *  performed.
+ *
+ *  @param[in] poly  expanded, square free polynomial
+ *  @param[in] syms  contains the symbols in the polynomial
+ *  @return          factorized polynomial
  */
 static ex factor_multivariate(const ex& poly, const exset& syms)
 {
 	exset::const_iterator s;
 	const ex& x = *syms.begin();
 
-	/* make polynomial primitive */
-	ex p = poly.collect(x);
-	ex cont = p.lcoeff(x);
-	for ( int i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
-		cont = gcd(cont, p.coeff(x,i));
-		if ( cont == 1 ) break;
-	}
-	ex pp = expand(normal(p / cont));
+	// 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);
 	}
 
-	/* factor leading coefficient */
+	// factor leading coefficient
 	ex vn = pp.collect(x).lcoeff(x);
 	ex vnlst;
 	if ( is_a<numeric>(vn) ) {
-		vnlst = lst(vn);
+		vnlst = lst{vn};
 	}
 	else {
 		ex vnfactors = factor(vn);
@@ -2012,7 +2254,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 	numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3;
 	vector<numeric> a(syms.size()-1, 0);
 
-	/* try now to factorize until we are successful */
+	// try now to factorize until we are successful
 	while ( true ) {
 
 		unsigned int trialcount = 0;
@@ -2022,10 +2264,10 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 		ex u, delta;
 		ex ufac, ufaclst;
 
-		/* try several evaluation points to reduce the number of modular factors */
+		// try several evaluation points to reduce the number of factors
 		while ( trialcount < maxtrials ) {
 
-			/* generate a set of valid evaluation points */
+			// generate a set of valid evaluation points
 			generate_set(pp, vn, syms, ex_to<lst>(vnlst), modulus, u, a);
 
 			ufac = factor_univariate(u, x, prime);
@@ -2034,19 +2276,19 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 			delta = ufaclst.op(0);
 
 			if ( factor_count <= 1 ) {
-				/* irreducible */
+				// irreducible
 				return poly;
 			}
 			if ( min_factor_count < 0 ) {
-				/* first time here */
+				// first time here
 				min_factor_count = factor_count;
 			}
 			else if ( min_factor_count == factor_count ) {
-				/* one less to try */
+				// one less to try
 				++trialcount;
 			}
 			else if ( min_factor_count > factor_count ) {
-				/* new minimum, reset trial counter */
+				// new minimum, reset trial counter
 				min_factor_count = factor_count;
 				trialcount = 0;
 			}
@@ -2055,11 +2297,15 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 		// determine true leading coefficients for the Hensel lifting
 		vector<ex> C(factor_count);
 		if ( is_a<numeric>(vn) ) {
+			// easy case
 			for ( size_t i=1; i<ufaclst.nops(); ++i ) {
 				C[i-1] = ufaclst.op(i).lcoeff(x);
 			}
-		}
-		else {
+		} else {
+			// difficult case.
+			// we use the property of the ftilde having a unique prime factor.
+			// details can be found in [Wan].
+			// calculate ftilde
 			vector<numeric> ftilde(vnlst.nops()-1);
 			for ( size_t i=0; i<ftilde.size(); ++i ) {
 				ex ft = vnlst.op(i+1);
@@ -2071,7 +2317,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 				}
 				ftilde[i] = ex_to<numeric>(ft);
 			}
-
+			// calculate D and C
 			vector<bool> used_flag(ftilde.size(), false);
 			vector<ex> D(factor_count, 1);
 			if ( delta == 1 ) {
@@ -2090,8 +2336,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 					}
 					C[i] = D[i] * prefac;
 				}
-			}
-			else {
+			} else {
 				for ( int i=0; i<factor_count; ++i ) {
 					numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(x));
 					for ( int j=ftilde.size()-1; j>=0; --j ) {
@@ -2115,7 +2360,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 					C[i] = D[i] * prefac;
 				}
 			}
-
+			// check if something went wrong
 			bool some_factor_unused = false;
 			for ( size_t i=0; i<used_flag.size(); ++i ) {
 				if ( !used_flag[i] ) {
@@ -2127,12 +2372,15 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 				continue;
 			}
 		}
-
+		
+		// multiply the remaining content of the univariate polynomial into the
+		// first factor
 		if ( delta != 1 ) {
 			C[0] = C[0] * delta;
 			ufaclst.let_op(1) = ufaclst.op(1) * delta;
 		}
 
+		// set up evaluation points
 		EvalPoint ep;
 		vector<EvalPoint> epv;
 		s = syms.begin();
@@ -2157,15 +2405,18 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 			l = l + 1;
 			pl = pl * prime;
 		}
+		
+		// set up modular factors (mod p^l)
 		cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l));
 		upvec modfactors(ufaclst.nops()-1);
 		for ( size_t i=1; i<ufaclst.nops(); ++i ) {
 			umodpoly_from_ex(modfactors[i-1], ufaclst.op(i), x, R);
 		}
 
+		// try Hensel lifting
 		ex res = hensel_multivar(pp, x, epv, prime, l, modfactors, C);
-		if ( res != lst() ) {
-			ex result = cont;
+		if ( res != lst{} ) {
+			ex result = cont * unit;
 			for ( size_t i=0; i<res.nops(); ++i ) {
 				result *= res.op(i).content(x) * res.op(i).unit(x);
 				result *= res.op(i).primpart(x);
@@ -2175,9 +2426,11 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
 	}
 }
 
+/** Finds all symbols in an expression. Used by factor_sqrfree() and factor().
+ */
 struct find_symbols_map : public map_function {
 	exset syms;
-	ex operator()(const ex& e)
+	ex operator()(const ex& e) override
 	{
 		if ( is_a<symbol>(e) ) {
 			syms.insert(e);
@@ -2187,6 +2440,9 @@ struct find_symbols_map : public map_function {
 	}
 };
 
+/** Factorizes a polynomial that is square free. It calls either the univariate
+ *  or the multivariate factorization functions.
+ */
 static ex factor_sqrfree(const ex& poly)
 {
 	// determine all symbols in poly
@@ -2200,11 +2456,11 @@ static ex factor_sqrfree(const ex& poly)
 		// univariate case
 		const ex& x = *(findsymbols.syms.begin());
 		if ( poly.ldegree(x) > 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 {
+		} else {
 			ex res = factor_univariate(poly, x);
 			return res;
 		}
@@ -2215,10 +2471,13 @@ static ex factor_sqrfree(const ex& poly)
 	return res;
 }
 
+/** Map used by factor() when factor_options::all is given to access all
+ *  subexpressions and to call factor() on them.
+ */
 struct apply_factor_map : public map_function {
 	unsigned options;
 	apply_factor_map(unsigned options_) : options(options_) { }
-	ex operator()(const ex& e)
+	ex operator()(const ex& e) override
 	{
 		if ( e.info(info_flags::polynomial) ) {
 			return factor(e, options);
@@ -2228,22 +2487,51 @@ struct apply_factor_map : public map_function {
 			for ( size_t i=0; i<e.nops(); ++i ) {
 				if ( e.op(i).info(info_flags::polynomial) ) {
 					s1 += e.op(i);
-				}
-				else {
+				} else {
 					s2 += e.op(i);
 				}
 			}
-			s1 = s1.eval();
-			s2 = s2.eval();
 			return factor(s1, options) + s2.map(*this);
 		}
 		return e.map(*this);
 	}
 };
 
-} // 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));
+		}
+	}
+}
 
-ex factor(const ex& poly, unsigned options)
+/** 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).
+ */
+static ex factor1(const ex& poly, unsigned options)
 {
 	// check arguments
 	if ( !poly.info(info_flags::polynomial) ) {
@@ -2262,56 +2550,43 @@ 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
 	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;
-			}
-			else {
-				res *= t;
+	ex res = 1;
+	factor_iter(sfpoly,
+		[&](const ex &f, const ex &k) {
+			if ( is_a<add>(f) ) {
+				res *= pow(factor_sqrfree(f), k);
+			} else {
+				// 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;
+}
+
+} // anonymous namespace
+
+/** 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