ginac/polynomial/heur_gcd_uvar.h

   1 #ifndef GINAC_UPOLY_HEUR_GCD
   2 #define GINAC_UPOLY_HEUR_GCD
   3 #include "upoly.hpp"
   4 #include "ring_traits.hpp"
   5 #include "normalize.tcc"
   6 #include "remainder.tcc"
   7 #include "eval_uvar.h"
   8 #include "interpolate_padic_uvar.h"
   9 #include <algorithm>
  10
  11 namespace GiNaC
  12 {
  13
  14 /// Compute GCD of primitive univariate polynomials.
  15 template<typename T> static bool
  16 heur_gcd_z_pp(T& g, const T& a, const T& b, unsigned max_tries = 66)
  17 {
  18         typedef typename T::value_type ring_t;
  19         const ring_t n73794 = get_ring_elt(b[0], 73794);
  20         const ring_t n27011 = get_ring_elt(b[0], 27011);
  21         const std::size_t maxdeg = std::max(degree(a), degree(b));
  22         T r, gg;
  23         gg.reserve(maxdeg);
  24         r.reserve(maxdeg);
  25
  26         // find the evaluation point
  27         ring_t xi = (std::min(max_coeff(a), max_coeff(b)) + 1) << 1;
  28
  29         do {
  30                 const ring_t av = eval(a, xi);
  31                 const ring_t bv = eval(b, xi);
  32                 const ring_t gamma = gcd(av, bv);
  33                 interpolate(gg, gamma, xi, maxdeg);
  34                 normalize_in_ring(gg);
  35                 remainder_in_ring(r, a, gg);
  36                 if (r.empty()) {
  37                         swap(g, gg);
  38                         return true;
  39                 }
  40                 // next evaluation point
  41                 xi = truncate1(xi*isqrt(isqrt(xi))*n73794, n27011);
  42         } while (--max_tries != 0);
  43         return false;
  44 }
  45
  46 template<typename T> static bool
  47 heur_gcd_z_priv(T& g, const T& a, const T& b, const unsigned max_tries = 66)
  48 {
  49         typedef typename T::value_type ring_t;
  50         ring_t acont, bcont;
  51         T a_(a), b_(b);
  52         normalize_in_ring(a_, &acont);
  53         normalize_in_ring(b_, &bcont);
  54         const ring_t gc = gcd(acont, bcont);
  55         bool found = heur_gcd_z_pp(g, a_, b_, max_tries);
  56         if (found) {
  57                 g *= gc;
  58                 return true;
  59         }
  60         return false;
  61 }
  62
  63 } // namespace GiNaC
  64
  65 #endif // GINAC_UPOLY_HEUR_GCD
  66