G_do_hoelder: fix case with real x values which are not of type cl_R.
[ginac.git] / ginac / polynomial / sr_gcd_uvar.h
1 /** @file sr_gcd_uvar.h
2  *
3  *  GCD function for univariate polynomials using PRS method. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2018 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_SR_GCD_H
24 #define GINAC_UPOLY_SR_GCD_H
25
26 #include "upoly.h"
27 #include "ring_traits.h"
28 #include "normalize.h"
29 #include "prem_uvar.h"
30
31 #include <limits>
32
33 namespace GiNaC {
34
35 /// Calculate GCD of two univariate polynomials @a a and @a b using
36 /// subresultant pseudo-remainder sequence method
37 template<typename T> static bool
38 sr_gcd_priv(T& g, T a, T b,
39             unsigned tries = std::numeric_limits<unsigned>::max())
40 {
41         // handle zero polynomials
42         if (a.empty()) {
43                 g.clear();
44                 return true;
45         }
46         if (b.empty()) {
47                 g.clear();
48                 return true;
49         }
50
51         typedef typename T::value_type ring_t;
52         if (degree(a) < degree(b))
53                 swap(a, b);
54         ring_t acont, bcont;
55         normalize_in_ring(a, &acont);
56         normalize_in_ring(b, &bcont);
57
58         const ring_t one = get_ring_elt(b[0], 1);
59
60         const ring_t gamma = gcd(acont, bcont);
61         if (degree(b) == 0) {
62                 g.resize(1);
63                 g[0] = gamma;
64                 return true;
65         }
66
67         T r(std::min(degree(a), degree(b)));
68
69         ring_t ri = one, psi = one;
70
71         do {
72                 std::size_t delta = degree(a) - degree(b);
73                 pseudoremainder(r, a, b);
74                 if (r.empty()) {
75                         normalize_in_ring(b);
76                         b *= gamma;
77                         swap(b, g);
78                         return true;
79                 }
80                 a = b;
81
82                 ring_t ri_psi_delta = delta > 0 ?  ri*expt_pos(psi, delta) : ri;
83
84                 bool divisible_p = divide(b, r, ri_psi_delta);
85                 bug_on(!divisible_p, "division failed: r = " << r <<
86                         ", ri = " << ri << ", psi = " << psi);
87
88                 if ((degree(b) == 0) && (degree(r) == 0)) {
89                         g.resize(1);
90                         g[0] = gamma;
91                         return true;
92                 }
93                 if (degree(b) == 0) {
94                         normalize_in_ring(r);
95                         r *= gamma;
96                         swap(r, g);
97                         return true;
98                 }
99
100                 ri = lcoeff(a);
101                 if (delta == 1)
102                         psi = ri;
103                 else if (delta) {
104                         const ring_t ri_delta = expt_pos(ri, delta);
105                         const ring_t psi_delta_1 = expt_pos(psi, delta - 1);
106                         bool sanity_check = div(psi, ri_delta, psi_delta_1);
107                         bug_on(!sanity_check, "division failed: ri = " << ri
108                                 << ", psi = " << psi << ", delta = " << delta);
109                 }
110                 if (--tries == 0)
111                         return false;
112         } while (true);
113 }
114
115 } // namespace GiNaC
116
117 #endif // GINAC_UPOLY_SR_GCD_H