*/
/*
- * GiNaC Copyright (C) 1999-2010 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
#include "add.h"
#include <algorithm>
-#include <cmath>
#include <limits>
#include <list>
#include <vector>
+#include <stack>
#ifdef DEBUGFACTOR
#include <ostream>
#endif
#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;
}
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<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;
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));
}
}
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_)
{
}
bool is_col_zero(size_t col) const
{
- mvec::const_iterator i = m.begin() + col;
for ( size_t rr=0; rr<r; ++rr ) {
std::size_t i = col + rr*c;
if ( !zerop(m[i]) ) {
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);
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;
}
mult[i] *= prime;
}
}
- }
- else {
+ } else {
umodpoly ap;
expt_1_over_p(a, prime, ap);
size_t previ = mult.size();
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);
}
}
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);
}
if ( alpha != 1 ) {
w = w / alpha;
}
- }
- else {
+ } else {
u.clear();
}
}
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!");
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;
}
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]);
}
}
lr[group] = lr[group] * cache[pos].back();
}
- }
- else {
+ } else {
lr[group] = lr[group] * factors[pos];
}
} while ( i < n );
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];
}
}
private:
umodpoly lr[2];
- vector< vector<umodpoly> > cache;
+ vector<vector<umodpoly>> cache;
upvec factors;
umodpoly one;
size_t n;
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;
minfactors = trialfactors.size();
lastp = prime;
trials = 1;
- }
- else {
+ } else {
++trials;
}
}
}
}
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];
}
}
tocheck.push(mf);
break;
}
- }
- else {
+ } else {
// not successful
if ( !part.next() ) {
// if no more combinations left, return polynomial as
{
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);
}
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);
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);
numeric n(R->retract(emod));
if ( n > halfmod ) {
return n-mod;
- }
- else {
+ } else {
return n;
}
}
{
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;
e = make_modular(buf, R);
}
}
- }
- else {
+ } else {
upvec amod;
for ( size_t i=0; i<a.size(); ++i ) {
umodpoly up;
if ( is_a<add>(c) ) {
nterms = c.nops();
z = c.op(0);
- }
- else {
+ } else {
nterms = 1;
z = c;
}
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;
- }
- else {
- lst res;
- return lst();
+ return lst(U.begin(), U.end());
+ } else {
+ return lst{};
}
}
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) ) {
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);
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].
}
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 ) {
// try Hensel lifting
ex res = hensel_multivar(pp, x, epv, prime, l, modfactors, C);
- if ( res != lst() ) {
+ 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);
*/
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);
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;
}
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);
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));
+ }
+ }
+}
-/** 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) ) {
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;
+ 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);
}
- else {
- res *= t;
- }
- }
- 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
git://www.ginac.de / ginac.git / blobdiff