ginac/polynomial/normalize.tcc

   1 #ifndef GINAC_UPOLY_NORMALIZE_TCC
   2 #define GINAC_UPOLY_NORMALIZE_TCC
   3 #include "upoly.hpp"
   4 #include "ring_traits.hpp"
   5 #include "debug.hpp"
   6
   7 namespace GiNaC
   8 {
   9
  10 /// Make the univariate polynomial @a a \in F[x] unit normal.
  11 /// F should be a field.
  12 /// Returns true if the polynomial @x is already unit normal, and false
  13 /// otherwise.
  14 static bool normalize_in_field(umodpoly& a, cln::cl_MI* content_ = 0)
  15 {
  16         if (a.size() == 0)
  17                 return true;
  18         if (lcoeff(a) == the_one(a[0])) {
  19                 if (content_)
  20                         *content_ = the_one(a[0]);
  21                 return true;
  22         }
  23
  24         const cln::cl_MI lc_1 = recip(lcoeff(a));
  25         for (std::size_t k = a.size(); k-- != 0; )
  26                 a[k] = a[k]*lc_1;
  27         if (content_)
  28                 *content_ = lc_1;
  29         return false;
  30 }
  31
  32 /// Make the univariate polynomial @a x unit normal. This version is used
  33 /// for rings which are not fields.
  34 /// Returns true if the polynomial @x is already unit normal, and false
  35 /// otherwise.
  36 template<typename T> bool
  37 normalize_in_ring(T& x, typename T::value_type* content_ = 0, int* unit_ = 0)
  38 {
  39         typedef typename T::value_type ring_t;
  40         static const ring_t one(1);
  41         if (x.empty())
  42                 return true;
  43
  44         bool something_changed = false;
  45         if (minusp(lcoeff(x))) {
  46                 something_changed = true;
  47                 if (unit_)
  48                         *unit_ = -1;
  49                 for (std::size_t i = x.size(); i-- != 0; )
  50                         x[i] = - x[i];
  51         }
  52
  53         if (degree(x) == 0) {
  54                 if (content_)
  55                         *content_ = x[0];
  56                 if (x[0] == one)
  57                         return something_changed;
  58                 x[0] = one;
  59                 return false; // initial polynomial was unit normal
  60         }
  61
  62         // Compute the gcd of coefficients
  63         ring_t content = lcoeff(x);
  64         // We want this function to be fast when applied to unit normal
  65         // polynomials. Hence we start from the leading coefficient.
  66         for (std::size_t i = x.size() - 1; i-- != 0; ) {
  67                 if (content == one) {
  68                         if (content_)
  69                                 *content_ = one;
  70                         return something_changed;
  71                 }
  72                 content = gcd(x[i], content);
  73         }
  74
  75         if (content == one) {
  76                 if (content_)
  77                         *content_ = one;
  78                 return something_changed;
  79         }
  80
  81         lcoeff(x) = one;
  82         for (std::size_t i = x.size() - 1; i-- != 0; )
  83                 x[i] = exquo(x[i], content);
  84
  85         if (content_)
  86                 *content_ = content;
  87         return false; // initial polynomial was not unit normal
  88 }
  89
  90 } // namespace GiNaC
  91
  92 #endif // GINAC_UPOLY_NORMALIZE_TCC
  93