Rewritten heuristic and PRS GCD for univariate polynomials, added benchmark.
[ginac.git] / ginac / polynomial / sr_gcd_uvar.h
1 #ifndef GINAC_UPOLY_SR_GCD_TCC
2 #define GINAC_UPOLY_SR_GCD_TCC
3 #include "upoly.hpp"
4 #include "ring_traits.hpp"
5 #include "normalize.tcc"
6 #include "prem_uvar.h"
7 #include <limits>
8
9 namespace GiNaC
10 {
11
12 /// Calculate GCD of two univariate polynomials @a a and @a b using
13 /// subresultant pseudo-remainder sequence method
14 template<typename T> static bool
15 sr_gcd_priv(T& g, T a, T b,
16             unsigned tries = std::numeric_limits<unsigned>::max())
17 {
18         // handle zero polynomials
19         if (a.empty()) {
20                 g.clear();
21                 return true;
22         }
23         if (b.empty()) {
24                 g.clear();
25                 return true;
26         }
27
28         typedef typename T::value_type ring_t;
29         if (degree(a) < degree(b))
30                 swap(a, b);
31         ring_t acont, bcont;
32         normalize_in_ring(a, &acont);
33         normalize_in_ring(b, &bcont);
34
35         const ring_t one = get_ring_elt(b[0], 1);
36
37         const ring_t gamma = gcd(acont, bcont);
38         if (degree(b) == 0) {
39                 g.resize(1);
40                 g[0] = gamma;
41                 return true;
42         }
43
44         T r(std::min(degree(a), degree(b)));
45
46         ring_t ri = one, psi = one;
47
48         do {
49                 std::size_t delta = degree(a) - degree(b);
50                 pseudoremainder(r, a, b);
51                 if (r.empty()) {
52                         normalize_in_ring(b);
53                         b *= gamma;
54                         swap(b, g);
55                         return true;
56                 }
57                 a = b;
58
59                 ring_t ri_psi_delta = delta > 0 ?  ri*expt_pos(psi, delta) : ri;
60
61                 bool divisible_p = divide(b, r, ri_psi_delta);
62                 bug_on(!divisible_p, "division failed: r = " << r <<
63                         ", ri = " << ri << ", psi = " << psi);
64
65                 if ((degree(b) == 0) && (degree(r) == 0)) {
66                         g.resize(1);
67                         g[0] = gamma;
68                         return true;
69                 }
70                 if (degree(b) == 0) {
71                         normalize_in_ring(r);
72                         r *= gamma;
73                         swap(r, g);
74                         return true;
75                 }
76
77                 ri = lcoeff(a);
78                 if (delta == 1)
79                         psi = ri;
80                 else if (delta) {
81                         const ring_t ri_delta = expt_pos(ri, delta);
82                         const ring_t psi_delta_1 = expt_pos(psi, delta - 1);
83                         bool sanity_check = div(psi, ri_delta, psi_delta_1);
84                         bug_on(!sanity_check, "division failed: ri = " << ri
85                                 << ", psi = " << psi << ", delta = " << delta);
86                 }
87                 if (--tries == 0)
88                         return false;
89         } while (true);
90 }
91
92 }
93
94 #endif // GINAC_UPOLY_SR_GCD_TCC
95