*/
/*
- * GiNaC Copyright (C) 1999-2022 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2024 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
return p[p.size() - 1];
}
+/** Make the polynomial unit normal (having unit normal leading coefficient).
+ *
+ * @param[in, out] a polynomial to make unit normal
+ * @return true if polynomial a was already unit normal, false otherwise
+ */
static bool normalize_in_field(umodpoly& a)
{
if (a.size() == 0)
return false;
}
+/** Remove leading zero coefficients from polynomial.
+ *
+ * @param[in, out] p polynomial from which the zero leading coefficients will be removed
+ * @param[in] hint assume all coefficients of order ≥ hint are zero
+ */
template<typename T> static void
canonicalize(T& p, const typename T::size_type hint = std::numeric_limits<typename T::size_type>::max())
{
- if (p.empty())
- return;
+ std::size_t i = min(p.size(), hint);
- std::size_t i = p.size() - 1;
- // Be fast if the polynomial is already canonicalized
- if (!zerop(p[i]))
- return;
-
- if (hint < p.size())
- i = hint;
+ while ( i-- && zerop(p[i]) ) { }
- bool is_zero = false;
- do {
- if (!zerop(p[i])) {
- ++i;
- break;
- }
- if (i == 0) {
- is_zero = true;
- break;
- }
- --i;
- } while (true);
-
- if (is_zero) {
- p.clear();
- return;
- }
-
- p.erase(p.begin() + i, p.end());
+ p.erase(p.begin() + i + 1, p.end());
}
// END COPY FROM UPOLY.H
-static void expt_pos(umodpoly& a, unsigned int q)
-{
- if ( a.empty() ) return;
- cl_MI zero = a[0].ring()->zero();
- int deg = degree(a);
- a.resize(degree(a)*q+1, zero);
- for ( int i=deg; i>0; --i ) {
- a[i*q] = a[i];
- a[i] = zero;
- }
-}
-
template<typename T> struct uvar_poly_p
{
static const bool value = false;
return e;
}
-/** Divides all coefficients of the polynomial a by the integer x.
+/** Divides all coefficients of the polynomial a by the positive integer x.
* All coefficients are supposed to be divisible by x. If they are not, the
- * the<cl_I> cast will raise an exception.
+ * division will raise an exception.
*
* @param[in,out] a polynomial of which the coefficients will be reduced by x
- * @param[in] x integer that divides the coefficients
+ * @param[in] x positive integer that divides the coefficients
*/
static void reduce_coeff(umodpoly& a, const cl_I& x)
{
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));
+ i = cl_MI(R, exquopos(c, x));
}
}
} while ( k-- );
fill(r.begin()+n, r.end(), a[0].ring()->zero());
- canonicalize(r);
+ canonicalize(r, n);
}
/** Calculates quotient of a/b.
static bool unequal_one(const umodpoly& a)
{
- if ( a.empty() ) return true;
return ( a.size() != 1 || a[0] != a[0].ring()->one() );
}
return equal_one(c);
}
+/** Computes w^q mod a.
+ * Uses theorem 2.1 from A.K.Lenstra's PhD thesis; see exercise 8.13 in [GCL].
+ *
+ * @param[in] w polynomial
+ * @param[in] a modulus polynomial
+ * @param[in] q common modulus of w and a
+ * @param[out] r result
+ */
+static void expt_pos_Q(const umodpoly& w, const umodpoly& a, unsigned int q, umodpoly& r)
+{
+ if ( w.empty() ) return;
+ cl_MI zero = w[0].ring()->zero();
+ int deg = degree(w);
+ umodpoly buf(deg*q+1, zero);
+ for ( size_t i=0; i<=deg; ++i ) {
+ buf[i*q] = w[i];
+ }
+ rem(buf, a, r);
+}
+
// END modular univariate polynomial code
////////////////////////////////////////////////////////////////////////////////
/** Calculates a^(1/prime).
*
- * @param[in] a polynomial
- * @param[in] prime prime number -> exponent 1/prime
- * @param[in] ap resulting polynomial
+ * @param[in] a polynomial
+ * @param[in] prime prime number -> exponent 1/prime
+ * @param[out] ap resulting polynomial
*/
static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
{
int nhalf = degree(a)/2;
int i = 1;
- umodpoly w(2);
- w[0] = R->zero();
- w[1] = R->one();
+ umodpoly w = {R->zero(), R->one()};
umodpoly x = w;
while ( i <= nhalf ) {
- expt_pos(w, q);
umodpoly buf;
- rem(w, a, buf);
+ expt_pos_Q(w, a, q, buf);
w = buf;
- umodpoly wx = w - x;
- gcd(a, wx, buf);
+ gcd(a, w - x, buf);
if ( unequal_one(buf) ) {
degrees.push_back(i);
ddfactors.push_back(buf);
- }
- if ( unequal_one(buf) ) {
umodpoly buf2;
div(a, buf, buf2);
a = buf2;
// calc bound B
int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
- cl_I maxmodulus = calc_bound(a) * ash(cl_I(1), maxdeg+1); // 2 * calc_bound(a) * 2^maxdeg
+ cl_I maxmodulus = ash(calc_bound(a), maxdeg+1); // = 2 * calc_bound(a) * 2^maxdeg
// step 1
cl_I alpha = lcoeff(a);
}
}
-/** Returns a new prime number.
+/** Returns a new small prime number.
*
- * @param[in] p prime number
- * @return next prime number after p
+ * @param[in] n an integer
+ * @return smallest prime greater than n
*/
-static unsigned int next_prime(unsigned int p)
-{
- static vector<unsigned int> primes;
- if (primes.empty()) {
- primes = {3, 5, 7};
- }
- 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;
- candidate += 2;
- i=-1;
- }
- primes.push_back(candidate);
- if (candidate > p)
+static unsigned int next_prime(unsigned int n)
+{
+ static vector<unsigned int> primes = {2, 3, 5, 7};
+ unsigned int candidate = primes.back();
+ while (primes.back() <= n) {
+ candidate += 2;
+ bool is_prime = true;
+ for (size_t i=1; primes[i]*primes[i]<=candidate; ++i) {
+ if (candidate % primes[i] == 0) {
+ is_prime = false;
break;
+ }
}
- return candidate;
+ if (is_prime)
+ primes.push_back(candidate);
}
for (auto & it : primes) {
- if ( it > p ) {
+ if ( it > n ) {
return it;
}
}
numeric prefac = ex_to<numeric>(ufaclst[i+1].lcoeff(x));
for ( int j=ftilde.size()-1; j>=0; --j ) {
int count = 0;
- while ( irem(prefac, ftilde[j]) == 0 ) {
- prefac = iquo(prefac, ftilde[j]);
+ numeric q;
+ while ( irem(prefac, ftilde[j], q) == 0 ) {
+ prefac = q;
++count;
}
if ( count ) {
numeric prefac = ex_to<numeric>(ufaclst[i+1].lcoeff(x));
for ( int j=ftilde.size()-1; j>=0; --j ) {
int count = 0;
- while ( irem(prefac, ftilde[j]) == 0 ) {
- prefac = iquo(prefac, ftilde[j]);
+ numeric q;
+ while ( irem(prefac, ftilde[j], q) == 0 ) {
+ prefac = q;
++count;
}
while ( irem(ex_to<numeric>(delta)*prefac, ftilde[j]) == 0 ) {
}
// set up evaluation points
- EvalPoint ep;
vector<EvalPoint> epv;
auto s = syms_wox.begin();
for ( size_t i=0; i<a.size(); ++i ) {
- ep.x = *s++;
- ep.evalpoint = a[i].to_int();
- epv.push_back(ep);
+ epv.emplace_back(EvalPoint{*s++, a[i].to_int()});
}
// calc bound p^l
maxdeg = ufaclst[i].degree(x);
}
}
- cl_I B = calc_bound(u, x) * ash(cl_I(1), maxdeg+1); // 2 * calc_bound(u,x) * 2^maxdeg
+ cl_I B = ash(calc_bound(u, x), maxdeg+1); // = 2 * calc_bound(u,x) * 2^maxdeg
cl_I l = 1;
cl_I pl = prime;
while ( pl < B ) {
copy_if(syms.begin(), syms.end(),
inserter(syms_wox, syms_wox.end()), [x](const ex& y){ return y != x; });
- factorization_ctx ctx = {.poly = poly, .x = x,
- .syms_wox = syms_wox};
+ factorization_ctx ctx{poly, x, syms_wox};
// make polynomial primitive
poly.unitcontprim(x, ctx.unit, ctx.cont, ctx.pp);
ctx.vn = ctx.pp.collect(x).lcoeff(x);
ctx.vnlst = put_factors_into_vec(factor(ctx.vn));
- ctx.modulus = (ctx.vnlst.size() > 3) ? ctx.vnlst.size() : 3;
+ ctx.modulus = (ctx.vnlst.size() > 3) ? ctx.vnlst.size() : numeric(3);
ctx_in_x.push_back(ctx);
}
git://www.ginac.de / ginac.git / blobdiff