ginac/polynomial/sr_gcd_uvar.h

   1 #ifndef GINAC_UPOLY_SR_GCD_TCC
   2 #define GINAC_UPOLY_SR_GCD_TCC
   3 #include "upoly.hpp"
   4 #include "ring_traits.hpp"
   5 #include "normalize.tcc"
   6 #include "prem_uvar.h"
   7 #include <limits>
   8
   9 namespace GiNaC
  10 {
  11
  12 /// Calculate GCD of two univariate polynomials @a a and @a b using
  13 /// subresultant pseudo-remainder sequence method
  14 template<typename T> static bool
  15 sr_gcd_priv(T& g, T a, T b,
  16             unsigned tries = std::numeric_limits<unsigned>::max())
  17 {
  18         // handle zero polynomials
  19         if (a.empty()) {
  20                 g.clear();
  21                 return true;
  22         }
  23         if (b.empty()) {
  24                 g.clear();
  25                 return true;
  26         }
  27
  28         typedef typename T::value_type ring_t;
  29         if (degree(a) < degree(b))
  30                 swap(a, b);
  31         ring_t acont, bcont;
  32         normalize_in_ring(a, &acont);
  33         normalize_in_ring(b, &bcont);
  34
  35         const ring_t one = get_ring_elt(b[0], 1);
  36
  37         const ring_t gamma = gcd(acont, bcont);
  38         if (degree(b) == 0) {
  39                 g.resize(1);
  40                 g[0] = gamma;
  41                 return true;
  42         }
  43
  44         T r(std::min(degree(a), degree(b)));
  45
  46         ring_t ri = one, psi = one;
  47
  48         do {
  49                 std::size_t delta = degree(a) - degree(b);
  50                 pseudoremainder(r, a, b);
  51                 if (r.empty()) {
  52                         normalize_in_ring(b);
  53                         b *= gamma;
  54                         swap(b, g);
  55                         return true;
  56                 }
  57                 a = b;
  58
  59                 ring_t ri_psi_delta = delta > 0 ?  ri*expt_pos(psi, delta) : ri;
  60
  61                 bool divisible_p = divide(b, r, ri_psi_delta);
  62                 bug_on(!divisible_p, "division failed: r = " << r <<
  63                         ", ri = " << ri << ", psi = " << psi);
  64
  65                 if ((degree(b) == 0) && (degree(r) == 0)) {
  66                         g.resize(1);
  67                         g[0] = gamma;
  68                         return true;
  69                 }
  70                 if (degree(b) == 0) {
  71                         normalize_in_ring(r);
  72                         r *= gamma;
  73                         swap(r, g);
  74                         return true;
  75                 }
  76
  77                 ri = lcoeff(a);
  78                 if (delta == 1)
  79                         psi = ri;
  80                 else if (delta) {
  81                         const ring_t ri_delta = expt_pos(ri, delta);
  82                         const ring_t psi_delta_1 = expt_pos(psi, delta - 1);
  83                         bool sanity_check = div(psi, ri_delta, psi_delta_1);
  84                         bug_on(!sanity_check, "division failed: ri = " << ri
  85                                 << ", psi = " << psi << ", delta = " << delta);
  86                 }
  87                 if (--tries == 0)
  88                         return false;
  89         } while (true);
  90 }
  91
  92 }
  93
  94 #endif // GINAC_UPOLY_SR_GCD_TCC
  95