]> www.ginac.de Git - ginac.git/blob - ginac/polynomial/prem_uvar.h
Fix power::subs() in some special cases.
[ginac.git] / 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