ginac/polynomial/mgcd.cpp

   1 #include "chinrem_gcd.h"
   2 #include <cln/integer.h>
   3 #include "pgcd.h"
   4 #include "collect_vargs.h"
   5 #include "primes_factory.h"
   6 #include "divide_in_z_p.h"
   7 #include "poly_cra.h"
   8
   9 namespace GiNaC
  10 {
  11
  12 static cln::cl_I extract_integer_content(ex& Apr, const ex& A)
  13 {
  14         static const cln::cl_I n1(1);
  15         const numeric icont_ = A.integer_content();
  16         const cln::cl_I icont = cln::the<cln::cl_I>(icont_.to_cl_N());
  17         if (icont != 1) {
  18                 Apr = (A/icont_).expand();
  19                 return icont;
  20         } else {
  21                 Apr = A;
  22                 return n1;
  23         }
  24 }
  25
  26 ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
  27 {
  28         ex A, B;
  29         const cln::cl_I a_icont = extract_integer_content(A, A_);
  30         const cln::cl_I b_icont = extract_integer_content(B, B_);
  31         const cln::cl_I c = cln::gcd(a_icont, b_icont);
  32
  33         const cln::cl_I a_lc = integer_lcoeff(A, vars);
  34         const cln::cl_I b_lc = integer_lcoeff(B, vars);
  35         const cln::cl_I g_lc = cln::gcd(a_lc, b_lc);
  36
  37         const ex& x(vars.back());
  38         int n = std::min(A.degree(x), B.degree(x));
  39         const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient());
  40         const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient());
  41         const cln::cl_I lcoeff_limit = (cln::cl_I(1) << n)*cln::abs(g_lc)*
  42                 std::min(A_max_coeff, B_max_coeff);
  43
  44         cln::cl_I q = 0;
  45         ex H = 0;
  46
  47         long p;
  48         primes_factory pfactory;
  49         while (true) {
  50                 bool has_primes = pfactory(p, g_lc);
  51                 if (!has_primes)
  52                         throw chinrem_gcd_failed();
  53
  54                 const numeric pnum(p);
  55                 ex Ap = A.smod(pnum);
  56                 ex Bp = B.smod(pnum);
  57                 ex Cp = pgcd(Ap, Bp, vars, p);
  58
  59                 const cln::cl_I g_lcp = smod(g_lc, p);
  60                 const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars);
  61                 const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p);
  62                 Cp = (Cp*numeric(nlc)).expand().smod(pnum);
  63                 int cp_deg = Cp.degree(x);
  64                 if (cp_deg == 0)
  65                         return numeric(g_lc);
  66                 if (zerop(q)) {
  67                         H = Cp;
  68                         n = cp_deg;
  69                         q = p;
  70                 } else {
  71                         if (cp_deg == n) {
  72                                 ex H_next = chinese_remainder(H, q, Cp, p);
  73                                 q = q*cln::cl_I(p);
  74                                 H = H_next;
  75                         } else if (cp_deg < n) {
  76                                 // all previous homomorphisms are unlucky
  77                                 q = p;
  78                                 H = Cp;
  79                                 n = cp_deg;
  80                         } else {
  81                                 // dp_deg > d_deg: current prime is bad
  82                         }
  83                 }
  84                 if (q < lcoeff_limit)
  85                         continue; // don't bother to do division checks
  86                 ex C, dummy1, dummy2;
  87                 extract_integer_content(C, H);
  88                 if (divide_in_z_p(A, C, dummy1, vars, 0) &&
  89                                 divide_in_z_p(B, C, dummy2, vars, 0))
  90                         return (numeric(c)*C).expand();
  91                 // else: try more primes
  92         }
  93 }
  94
  95 } // namespace GiNaC
  96