ginac/polynomial/mod_gcd.cpp

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