3 * Chinese remainder algorithm. */
6 * GiNaC Copyright (C) 1999-2010 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 #include "operators.h"
24 #include "chinrem_gcd.h"
26 #include "collect_vargs.h"
27 #include "primes_factory.h"
28 #include "divide_in_z_p.h"
30 #include <numeric> // std::accumulate
32 #include <cln/integer.h>
33 #include <cln/integer_ring.h>
34 #include <cln/rational.h>
35 #include <cln/rational_ring.h>
39 static cln::cl_I extract_integer_content(ex& Apr, const ex& A)
41 static const cln::cl_I n1(1);
42 const numeric icont_ = A.integer_content();
43 if (cln::instanceof(icont_.to_cl_N(), cln::cl_I_ring)) {
44 const cln::cl_I icont = cln::the<cln::cl_I>(icont_.to_cl_N());
46 Apr = (A/icont_).expand();
53 if (cln::instanceof(icont_.to_cl_N(), cln::cl_RA_ring)) {
54 Apr = (A/icont_).expand();
55 // A is a polynomail over rationals, so GCD is defined
56 // up to arbitrary rational number.
59 GINAC_ASSERT(NULL == "expected polynomial over integers or rationals");
62 ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
65 const cln::cl_I a_icont = extract_integer_content(A, A_);
66 const cln::cl_I b_icont = extract_integer_content(B, B_);
67 const cln::cl_I c = cln::gcd(a_icont, b_icont);
69 const cln::cl_I a_lc = integer_lcoeff(A, vars);
70 const cln::cl_I b_lc = integer_lcoeff(B, vars);
71 const cln::cl_I g_lc = cln::gcd(a_lc, b_lc);
73 const ex& x(vars.back());
74 exp_vector_t n = std::min(degree_vector(A, vars), degree_vector(B, vars));
75 const int nTot = std::accumulate(n.begin(), n.end(), 0);
76 const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient());
77 const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient());
79 const cln::cl_I lcoeff_limit = (cln::cl_I(1) << nTot)*cln::abs(g_lc)*
80 std::min(A_max_coeff, B_max_coeff);
87 primes_factory pfactory;
89 bool has_primes = pfactory(p, g_lc);
91 throw chinrem_gcd_failed();
93 const numeric pnum(p);
96 ex Cp = pgcd(Ap, Bp, vars, p);
98 const cln::cl_I g_lcp = smod(g_lc, p);
99 const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars);
100 const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p);
101 Cp = (Cp*numeric(nlc)).expand().smod(pnum);
102 exp_vector_t cp_deg = degree_vector(Cp, vars);
104 return numeric(g_lc);
111 ex H_next = chinese_remainder(H, q, Cp, p);
114 } else if (cp_deg < n) {
115 // all previous homomorphisms are unlucky
120 // dp_deg > d_deg: current prime is bad
123 if (q < lcoeff_limit)
124 continue; // don't bother to do division checks
125 ex C, dummy1, dummy2;
126 extract_integer_content(C, H);
127 if (divide_in_z_p(A, C, dummy1, vars, 0) &&
128 divide_in_z_p(B, C, dummy2, vars, 0))
129 return (numeric(c)*C).expand();
130 // else: try more primes