ginac/polynomial/mod_gcd.cpp

   1 #include "upoly.hpp"
   2 #include "gcd_euclid.tcc"
   3 #include "cra_garner.hpp"
   4 #include <cln/random.h>
   5 #include <cln/numtheory.h>
   6 #include <stdexcept>
   7 #include "debug.hpp"
   8
   9 namespace GiNaC
  10 {
  11 /**
  12  * @brief Remove the integer content from univariate polynomials A and B.
  13  *
  14  * As a byproduct compute the GCD of contents.
  15  */
  16 static void remove_content(upoly& A, upoly& B, upoly::value_type& c)
  17 {
  18         // remove the integer
  19         upoly::value_type acont, bcont;
  20         normalize_in_ring(A, &acont);
  21         normalize_in_ring(B, &bcont);
  22         c = gcd(acont, bcont);
  23 };
  24
  25 /// Check if @a candidate divides both @a A and @a B
  26 static bool
  27 do_division_check(const upoly& A, const upoly& B, const upoly& candidate)
  28 {
  29         upoly r1;
  30         remainder_in_ring(r1, A, candidate);
  31         if (r1.size() != 0)
  32                 return false;
  33
  34         upoly r2;
  35         remainder_in_ring(r2, B, candidate);
  36         if (r2.size() != 0)
  37                 return false;
  38
  39         return true;
  40 }
  41
  42 /**
  43  * Given two GCD candidates H \in Z/q[x], and C \in Z/p[x] (where p is a prime)
  44  * compute the candidate in Z/(q*p) with CRA (chinise remainder algorithm)
  45  *
  46  * @param H \in Z/q[x] GCD candidate, will be updated by this function
  47  * @param q modulus of H, will NOT be updated by this function
  48  * @param C \in Z/p[x] GCD candidate, will be updated by this function
  49  * @param p modulus of C
  50  */
  51 static void
  52 update_the_candidate(upoly& H, const upoly::value_type& q,
  53                      const umodpoly& C,
  54                      const upoly::value_type& p,
  55                      const cln::cl_modint_ring& R)
  56 {
  57         typedef upoly::value_type ring_t;
  58         std::vector<ring_t> moduli(2);
  59         moduli[0] = q;
  60         moduli[1] = p;
  61         if (H.size() < C.size())
  62                 H.resize(C.size());
  63
  64         for (std::size_t  i = C.size(); i-- != 0; ) {
  65                 std::vector<ring_t> coeffs(2);
  66                 coeffs[0] = H[i];
  67                 coeffs[1] = R->retract(C[i]);
  68                 H[i] = integer_cra(coeffs, moduli);
  69         }
  70 }
  71
  72 /// Convert Z/p[x] -> Z[x]
  73 static void retract(upoly& p, const umodpoly& cp,
  74                     const cln::cl_modint_ring& Rp)
  75 {
  76         p.resize(cp.size());
  77         for (std::size_t i = cp.size(); i-- != 0; )
  78                 p[i] = Rp->retract(cp[i]);
  79 }
  80
  81
  82 /// Find the prime which is > p, and does NOT divide g
  83 static void find_next_prime(cln::cl_I& p, const cln::cl_I& g)
  84 {
  85         do {
  86                 ++p;
  87                 p = nextprobprime(p);
  88         } while (zerop(mod(g, p)));
  89 }
  90
  91 /// Compute the GCD of univariate polynomials A, B \in Z[x]
  92 void mod_gcd(upoly& result, upoly A, upoly B)
  93 {
  94         typedef upoly::value_type ring_t;
  95         ring_t content_gcd;
  96         // remove the integer content
  97         remove_content(A, B, content_gcd);
  98
  99         // compute the coefficient bound for GCD(a, b)
 100         ring_t g = gcd(lcoeff(A), lcoeff(B));
 101         std::size_t max_gcd_degree = std::min(degree(A), degree(B));
 102         ring_t limit = (ring_t(1) << max_gcd_degree)*g*
 103                        std::min(max_coeff(A), max_coeff(B));
 104         ring_t q(0);
 105         upoly H;
 106
 107         int count = 0;
 108         const ring_t p_threshold = ring_t(1) << (8*sizeof(void *));
 109         ring_t p = isqrt(std::min(max_coeff(A), max_coeff(B)));
 110         while (true) {
 111                 if (count >= 8) {
 112                         count = 0;
 113                         if (p < p_threshold)
 114                                 p <<= 1;
 115                         else
 116                                 p = p + (p >> 4);
 117                 } else
 118                         ++count;
 119                 find_next_prime(p, g);
 120
 121                 // Map the polynomials onto Z/p[x]
 122                 cln::cl_modint_ring Rp = cln::find_modint_ring(p);
 123                 cln::cl_MI gp = Rp->canonhom(g);
 124                 umodpoly ap(A.size()), bp(B.size());
 125                 make_umodpoly(ap, A, Rp);
 126                 make_umodpoly(bp, B, Rp);
 127
 128                 // Compute the GCD in Z/p[x]
 129                 umodpoly cp;
 130                 gcd_euclid(cp, ap, bp);
 131                 bug_on(cp.size() == 0, "gcd(ap, bp) = 0, with ap = " <<
 132                                         ap << ", and bp = " << bp);
 133
 134
 135                 // Normalize the candidate so that its leading coefficient
 136                 // is g mod p
 137                 umodpoly::value_type norm_factor = gp*recip(lcoeff(cp));
 138                 bug_on(zerop(norm_factor), "division in a field give 0");
 139
 140                 lcoeff(cp) = gp;
 141                 for (std::size_t k = cp.size() - 1; k-- != 0; )
 142                         cp[k] = cp[k]*norm_factor;
 143
 144
 145                 // check for unlucky homomorphisms
 146                 if (degree(cp) < max_gcd_degree) {
 147                         q = p;
 148                         max_gcd_degree = degree(cp);
 149                         retract(H, cp, Rp);
 150                 } else {
 151                         update_the_candidate(H, q, cp, p, Rp);
 152                         q = q*p;
 153                 }
 154
 155                 if (q > limit) {
 156                         result = H;
 157                         normalize_in_ring(result);
 158                         // if H divides both A and B it's a GCD
 159                         if (do_division_check(A, B, result)) {
 160                                 result *= content_gcd;
 161                                 break;
 162                         }
 163                         // H does not divide A and/or B, look for
 164                         // another one
 165                 } else if (degree(cp) == 0) {
 166                         // Polynomials are relatively prime
 167                         result.resize(1);
 168                         result[0] = content_gcd;
 169                         break;
 170                 }
 171         }
 172 }
 173
 174 } // namespace GiNaC
 175