ginac/polynomial/mgcd.cpp

   1 /** @file mgcd.cpp
   2  *
   3  *  Chinese remainder algorithm. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2010 Johannes Gutenberg University Mainz, Germany
   7  *
   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.
  12  *
  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.
  17  *
  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
  21  */
  22
  23 #include "operators.h"
  24 #include "chinrem_gcd.h"
  25 #include "pgcd.h"
  26 #include "collect_vargs.h"
  27 #include "primes_factory.h"
  28 #include "divide_in_z_p.h"
  29 #include "poly_cra.h"
  30 #include <numeric> // std::accumulate
  31
  32 #include <cln/integer.h>
  33
  34 namespace GiNaC {
  35
  36 static cln::cl_I extract_integer_content(ex& Apr, const ex& A)
  37 {
  38         static const cln::cl_I n1(1);
  39         const numeric icont_ = A.integer_content();
  40         const cln::cl_I icont = cln::the<cln::cl_I>(icont_.to_cl_N());
  41         if (icont != 1) {
  42                 Apr = (A/icont_).expand();
  43                 return icont;
  44         } else {
  45                 Apr = A;
  46                 return n1;
  47         }
  48 }
  49
  50 ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
  51 {
  52         ex A, B;
  53         const cln::cl_I a_icont = extract_integer_content(A, A_);
  54         const cln::cl_I b_icont = extract_integer_content(B, B_);
  55         const cln::cl_I c = cln::gcd(a_icont, b_icont);
  56
  57         const cln::cl_I a_lc = integer_lcoeff(A, vars);
  58         const cln::cl_I b_lc = integer_lcoeff(B, vars);
  59         const cln::cl_I g_lc = cln::gcd(a_lc, b_lc);
  60
  61         const ex& x(vars.back());
  62         exp_vector_t n = std::min(degree_vector(A, vars), degree_vector(B, vars));
  63         const int nTot = std::accumulate(n.begin(), n.end(), 0);
  64         const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient());
  65         const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient());
  66
  67         const cln::cl_I lcoeff_limit = (cln::cl_I(1) << nTot)*cln::abs(g_lc)*
  68                 std::min(A_max_coeff, B_max_coeff);
  69
  70
  71         cln::cl_I q = 0;
  72         ex H = 0;
  73
  74         long p;
  75         primes_factory pfactory;
  76         while (true) {
  77                 bool has_primes = pfactory(p, g_lc);
  78                 if (!has_primes)
  79                         throw chinrem_gcd_failed();
  80
  81                 const numeric pnum(p);
  82                 ex Ap = A.smod(pnum);
  83                 ex Bp = B.smod(pnum);
  84                 ex Cp = pgcd(Ap, Bp, vars, p);
  85
  86                 const cln::cl_I g_lcp = smod(g_lc, p);
  87                 const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars);
  88                 const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p);
  89                 Cp = (Cp*numeric(nlc)).expand().smod(pnum);
  90                 exp_vector_t cp_deg = degree_vector(Cp, vars);
  91                 if (zerop(cp_deg))
  92                         return numeric(g_lc);
  93                 if (zerop(q)) {
  94                         H = Cp;
  95                         n = cp_deg;
  96                         q = p;
  97                 } else {
  98                         if (cp_deg == n) {
  99                                 ex H_next = chinese_remainder(H, q, Cp, p);
 100                                 q = q*cln::cl_I(p);
 101                                 H = H_next;
 102                         } else if (cp_deg < n) {
 103                                 // all previous homomorphisms are unlucky
 104                                 q = p;
 105                                 H = Cp;
 106                                 n = cp_deg;
 107                         } else {
 108                                 // dp_deg > d_deg: current prime is bad
 109                         }
 110                 }
 111                 if (q < lcoeff_limit)
 112                         continue; // don't bother to do division checks
 113                 ex C, dummy1, dummy2;
 114                 extract_integer_content(C, H);
 115                 if (divide_in_z_p(A, C, dummy1, vars, 0) &&
 116                                 divide_in_z_p(B, C, dummy2, vars, 0))
 117                         return (numeric(c)*C).expand();
 118                 // else: try more primes
 119         }
 120 }
 121
 122 } // namespace GiNaC