From: Richard Kreckel Date: Mon, 14 Nov 2022 23:08:28 +0000 (+0100) Subject: Refactor internal calc_bound() routines. X-Git-Tag: release_1-8-5~8 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=commitdiff_plain;h=7ddcbce353f5811a8967297eece3fdff7ae04f50;p=ginac.git Refactor internal calc_bound() routines. Reduce them to computation of absolute values of the polynomial's roots and document that this is known as Landau's inequality. The multiplication with a power of two is done by the caller. Also remove the internal fallback computation using the maxcoeff since this cannot happen. (In the worst case, there is only one non-zero coefficient 'a'. Then we compute sqrt(a^2) which is not smaller than a. Otherwise, we compute sqrt(a^2 + x) > a.) --- diff --git a/ginac/factor.cpp b/ginac/factor.cpp index 9bd105ec..bbefcf3f 100644 --- a/ginac/factor.cpp +++ b/ginac/factor.cpp @@ -1187,38 +1187,30 @@ static upoly replace_lc(const upoly& poly, const cl_I& lc) return r; } -/** Calculates the bound for the modulus. - * See [Mig]. +/** Calculates bound for the product of absolute values (modulus) of the roots. + * Uses Landau's inequality, see [Mig]. */ -static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg) +static inline cl_I calc_bound(const ex& a, const ex& x) { - cl_I maxcoeff = 0; - cl_R coeff = 0; + cl_R radicand = 0; for ( int i=a.degree(x); i>=a.ldegree(x); --i ) { cl_I aa = abs(the(ex_to(a.coeff(x, i)).to_cl_N())); - if ( aa > maxcoeff ) maxcoeff = aa; - coeff = coeff + square(aa); + radicand = radicand + square(aa); } - cl_I coeffnorm = ceiling1(the(cln::sqrt(coeff))); - cl_I B = coeffnorm * ash(cl_I(1), maxdeg); // coeffnorm * 2^maxdeg - return ( B > maxcoeff ) ? B : maxcoeff; + return ceiling1(the(cln::sqrt(radicand))); } -/** Calculates the bound for the modulus. - * See [Mig]. +/** Calculates bound for the product of absolute values (modulus) of the roots. + * Uses Landau's inequality, see [Mig]. */ -static inline cl_I calc_bound(const upoly& a, int maxdeg) +static inline cl_I calc_bound(const upoly& a) { - cl_I maxcoeff = 0; - cl_R coeff = 0; + cl_R radicand = 0; for ( int i=degree(a); i>=0; --i ) { cl_I aa = abs(a[i]); - if ( aa > maxcoeff ) maxcoeff = aa; - coeff = coeff + square(aa); + radicand = radicand + square(aa); } - cl_I coeffnorm = ceiling1(the(cln::sqrt(coeff))); - cl_I B = coeffnorm * ash(cl_I(1), maxdeg); // coeffnorm * 2^maxdeg - return ( B > maxcoeff ) ? B : maxcoeff; + return ceiling1(the(cln::sqrt(radicand))); } /** Hensel lifting as used by factor_univariate(). @@ -1240,7 +1232,7 @@ static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, // calc bound B int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_); - cl_I maxmodulus = 2*calc_bound(a, maxdeg); + cl_I maxmodulus = calc_bound(a) * ash(cl_I(1), maxdeg+1); // 2 * calc_bound(a) * 2^maxdeg // step 1 cl_I alpha = lcoeff(a); @@ -1343,7 +1335,7 @@ static unsigned int next_prime(unsigned int p) throw logic_error("next_prime: should not reach this point!"); } -/** Manages the splitting a vector of of modular factors into two partitions. +/** Manages the splitting of a vector of modular factors into two partitions. */ class factor_partition { @@ -2394,7 +2386,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms) maxdeg = ufaclst[i].degree(x); } } - cl_I B = 2*calc_bound(u, x, maxdeg); + cl_I B = calc_bound(u, x) * ash(cl_I(1), maxdeg+1); // 2 * calc_bound(u,x) * 2^maxdeg cl_I l = 1; cl_I pl = prime; while ( pl < B ) {