ginac/polynomial/mgcd.cpp

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