ginac/polynomial/prem_uvar.h

   1 #ifndef GINAC_POLYNOMIAL_PREM_TCC
   2 #define GINAC_POLYNOMIAL_PREM_TCC
   3 #include "upoly.hpp"
   4 #include "debug.hpp"
   5 #include "remainder.tcc"
   6
   7 namespace GiNaC
   8 {
   9
  10 /// Compute the pseudo-remainder of univariate polynomials @a a and @a b
  11 /// Pseudo remainder \f$r(x)\f$ is defined as
  12 /// \f$\beta^l a(x) = b(x) q(x) + r(x) \f$, where \f$\beta\f$ is leading
  13 /// coefficient of \f$b(x)\f$ and \f$l = degree(a) - degree(b) + 1\f$
  14 /// FIXME: this implementation is extremely dumb.
  15 template<typename T> bool pseudoremainder(T& r, const T& a, const T& b)
  16 {
  17         typedef typename T::value_type ring_t;
  18         bug_on(b.size() == 0, "division by zero");
  19         if (a.size() == 1 && b.size() == 1) {
  20                 r.clear();
  21                 return true;
  22         }
  23         if (a.size() == 1) {
  24                 r = b;
  25                 return false;
  26         }
  27         if (degree(b) > degree(a)) {
  28                 r = b;
  29                 return false;
  30         }
  31
  32         const ring_t one = get_ring_elt(b[0], 1);
  33         const std::size_t l = degree(a) - degree(b) + 1;
  34         const ring_t blcoeff = lcoeff(b);
  35         const ring_t b_lth = expt_pos(blcoeff, l);
  36         if (b_lth == one)
  37                 return remainder_in_ring(r, a, b);
  38
  39         T a_(a);
  40         a_ *= b_lth;
  41         return remainder_in_ring(r, a_, b);
  42 }
  43
  44 }
  45
  46 #endif // GINAC_POLYNOMIAL_PREM_TCC
  47