Rewritten heuristic and PRS GCD for univariate polynomials, added benchmark.
[ginac.git] / 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