ginac/polynomial/heur_gcd_uvar.h

   1 /** @file heur_gcd_uvar.h
   2  *
   3  *  Heuristic GCD code. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2017 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 #ifndef GINAC_UPOLY_HEUR_GCD_H
  24 #define GINAC_UPOLY_HEUR_GCD_H
  25
  26 #include "upoly.h"
  27 #include "ring_traits.h"
  28 #include "normalize.h"
  29 #include "remainder.h"
  30 #include "eval_uvar.h"
  31 #include "interpolate_padic_uvar.h"
  32
  33 #include <algorithm>
  34
  35 namespace GiNaC {
  36
  37 /// Compute GCD of primitive univariate polynomials.
  38 template<typename T> static bool
  39 heur_gcd_z_pp(T& g, const T& a, const T& b, unsigned max_tries = 66)
  40 {
  41         typedef typename T::value_type ring_t;
  42         const ring_t n73794 = get_ring_elt(b[0], 73794);
  43         const ring_t n27011 = get_ring_elt(b[0], 27011);
  44         const std::size_t maxdeg = std::max(degree(a), degree(b));
  45         T r, gg;
  46         gg.reserve(maxdeg);
  47         r.reserve(maxdeg);
  48
  49         // find the evaluation point
  50         ring_t xi = (std::min(max_coeff(a), max_coeff(b)) + 1) << 1;
  51
  52         do {
  53                 const ring_t av = eval(a, xi);
  54                 const ring_t bv = eval(b, xi);
  55                 const ring_t gamma = gcd(av, bv);
  56                 interpolate(gg, gamma, xi, maxdeg);
  57                 normalize_in_ring(gg);
  58                 remainder_in_ring(r, a, gg);
  59                 if (r.empty()) {
  60                         swap(g, gg);
  61                         return true;
  62                 }
  63                 // next evaluation point
  64                 xi = truncate1(xi*isqrt(isqrt(xi))*n73794, n27011);
  65         } while (--max_tries != 0);
  66         return false;
  67 }
  68
  69 template<typename T> static bool
  70 heur_gcd_z_priv(T& g, const T& a, const T& b, const unsigned max_tries = 66)
  71 {
  72         typedef typename T::value_type ring_t;
  73         ring_t acont, bcont;
  74         T a_(a), b_(b);
  75         normalize_in_ring(a_, &acont);
  76         normalize_in_ring(b_, &bcont);
  77         const ring_t gc = gcd(acont, bcont);
  78         bool found = heur_gcd_z_pp(g, a_, b_, max_tries);
  79         if (found) {
  80                 g *= gc;
  81                 return true;
  82         }
  83         return false;
  84 }
  85
  86 } // namespace GiNaC
  87
  88 #endif // ndef GINAC_UPOLY_HEUR_GCD_H