ginac/polynomial/upoly.h

   1 /** @file upoly.h
   2  *
   3  *  Interface to polynomials with integer and modular coefficients. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2015 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_H
  24 #define GINAC_UPOLY_H
  25
  26 #include "ring_traits.h"
  27 #include "debug.h"
  28 #include "compiler.h"
  29
  30 #include <cln/integer.h>
  31 #include <cln/integer_io.h>
  32 #include <cln/modinteger.h>
  33 #include <cstddef>
  34 #include <iterator>
  35 #include <limits>
  36 #include <stdexcept>
  37 #include <vector>
  38
  39 namespace GiNaC {
  40
  41 typedef std::vector<cln::cl_I> upoly;
  42 typedef std::vector<cln::cl_MI> umodpoly;
  43
  44 template<typename T> static std::size_t degree(const T& p)
  45 {
  46         return p.size() - 1;
  47 }
  48
  49 template<typename T> static typename T::value_type lcoeff(const T& p)
  50 {
  51         bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
  52         return p[p.size() - 1];
  53 }
  54
  55 template<typename T> static typename T::reference lcoeff(T& p)
  56 {
  57         bug_on(p.empty(), "lcoeff of a zero polynomial is undefined");
  58         return p[p.size() - 1];
  59 }
  60
  61 template<typename T> static typename T::value_type max_coeff(const T& p)
  62 {
  63         bug_on(p.empty(), "max_coeff of a zero polynomial is undefined");
  64         typename T::value_type curr = p[p.size() - 1];
  65         for (std::size_t i = p.size(); i-- != 0; ) {
  66                 if (p[i] > curr)
  67                         curr = p[i];
  68         }
  69         return curr;
  70 }
  71
  72
  73
  74 // Canonicalize the polynomial, i.e. remove leading zero coefficients
  75 template<typename T> static void
  76 canonicalize(T& p, const typename T::size_type hint =
  77                    std::numeric_limits<typename T::size_type>::max())
  78 {
  79         if (p.empty())
  80                 return;
  81
  82         std::size_t i = p.size() - 1;
  83         // Be fast if the polynomial is already canonicalized
  84         if (!zerop(p[i]))
  85                 return;
  86
  87         if (hint < p.size())
  88                 i = hint;
  89
  90         bool is_zero = false;
  91         do {
  92                 if (!zerop(p[i])) {
  93                         ++i;
  94                         break;
  95                 }
  96                 if (i == 0) {
  97                         is_zero = true;
  98                         break;
  99                 }
 100                 --i;
 101         } while (true);
 102
 103         if (is_zero) {
 104                 p.clear();
 105                 return;
 106         }
 107
 108         bug_on(!zerop(p.at(i)), "p[" << i << "] = " << p[i] << " != 0 would be erased.");
 109
 110         typename T::const_iterator it = p.begin() + i;
 111         for (std::size_t k = i; it != p.end(); ++it, ++k) {
 112                 bug_on(!zerop(*it), "p[" << k << "] = " << p[k] << " != 0 "
 113                                    "would be erased.");
 114         }
 115
 116         p.erase(p.begin() + i, p.end());
 117
 118         bug_on(!p.empty() && zerop(lcoeff(p)), "oops, lcoeff(p) = 0");
 119 }
 120
 121 // Multiply a univariate polynomial p \in R[x] with a constant c \in R
 122 template<typename T> static T&
 123 operator*=(T& p, const typename T::value_type& c)
 124 {
 125         if (p.empty())
 126                 return p;
 127         if (zerop(c)) {
 128                 p.clear();
 129                 return p;
 130         }
 131         if (c == the_one(p[0]))
 132                 return p;
 133
 134         for (std::size_t i = p.size(); i-- != 0; )
 135                 p[i] = p[i]*c;
 136         canonicalize(p);
 137         return p;
 138 }
 139
 140 /// Divide the polynomial @a p by the ring element @a c, put the result
 141 /// into @a r. If @a p is not divisible by @a c @a r is undefined.
 142 template<typename T> bool divide(T& r, const T& p, const typename T::value_type& c)
 143 {
 144         if (p.empty()) {
 145                 r.clear();
 146                 return true;
 147         }
 148         if (r.size() != p.size())
 149                 r.resize(p.size());
 150         bool divisible = true;
 151         for (std::size_t i = p.size(); i-- != 0; ) {
 152                 divisible = div(r[i], p[i], c);
 153                 if (!divisible)
 154                         break;
 155         }
 156         return divisible;
 157 }
 158
 159 template<typename T> bool divide(T& p, const typename T::value_type& c)
 160 {
 161         if (p.empty())
 162                 return true;
 163         T r(p.size());
 164         bool divisible = divide(r, p, c);
 165         if (divisible)
 166                 swap(p, r);
 167         return divisible;
 168 }
 169
 170 // Convert Z[x] -> Z/p[x]
 171
 172 static inline void
 173 make_umodpoly(umodpoly& up, const upoly& p, const cln::cl_modint_ring& R)
 174 {
 175         for (std::size_t i = p.size(); i-- != 0; )
 176                 up[i] = R->canonhom(p[i]);
 177         canonicalize(up);
 178 }
 179
 180 } // namespace GiNaC
 181
 182 #endif // GINAC_UPOLY_H