*/
/*
- * GiNaC Copyright (C) 1999-2015 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
#include "add.h"
#include <algorithm>
-#include <cmath>
#include <limits>
#include <list>
#include <vector>
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;
{
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;
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);
}
}
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{};
}
}
*/
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
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]);
// 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++;
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 {
}
};
-} // 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) ) {
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
git://www.ginac.de / ginac.git / blobdiff