ginac/polynomial/prem_uvar.h

   1 /** @file prem_uvar.h
   2  *
   3  *  Function to calculate the pseudo-remainder. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2022 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_POLYNOMIAL_PREM_H
  24 #define GINAC_POLYNOMIAL_PREM_H
  25
  26 #include "upoly.h"
  27 #include "debug.h"
  28 #include "remainder.h"
  29
  30 namespace GiNaC {
  31
  32 /// Compute the pseudo-remainder of univariate polynomials @a a and @a b
  33 /// Pseudo remainder \f$r(x)\f$ is defined as
  34 /// \f$\beta^l a(x) = b(x) q(x) + r(x) \f$, where \f$\beta\f$ is leading
  35 /// coefficient of \f$b(x)\f$ and \f$l = degree(a) - degree(b) + 1\f$
  36 /// FIXME: this implementation is extremely dumb.
  37 template<typename T> bool pseudoremainder(T& r, const T& a, const T& b)
  38 {
  39         typedef typename T::value_type ring_t;
  40         bug_on(b.size() == 0, "division by zero");
  41         if (a.size() == 1 && b.size() == 1) {
  42                 r.clear();
  43                 return true;
  44         }
  45         if (a.size() == 1) {
  46                 r = b;
  47                 return false;
  48         }
  49         if (degree(b) > degree(a)) {
  50                 r = b;
  51                 return false;
  52         }
  53
  54         const ring_t one = get_ring_elt(b[0], 1);
  55         const std::size_t l = degree(a) - degree(b) + 1;
  56         const ring_t blcoeff = lcoeff(b);
  57         const ring_t b_lth = expt_pos(blcoeff, l);
  58         if (b_lth == one)
  59                 return remainder_in_ring(r, a, b);
  60
  61         T a_(a);
  62         a_ *= b_lth;
  63         return remainder_in_ring(r, a_, b);
  64 }
  65
  66 } // namespace GiNaC
  67
  68 #endif // GINAC_POLYNOMIAL_PREM_H