G_do_hoelder: fix case with real x values which are not of type cl_R.
[ginac.git] / ginac / polynomial / divide_in_z_p.cpp
1 /** @file divide_in_z_p.cpp
2  *
3  *  Implementation of polynomial division in Z/Zp. */
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 #include "add.h"
24 #include "operators.h"
25 #include "power.h"
26 #include "smod_helpers.h"
27
28 namespace GiNaC {
29
30 /** 
31  * Exact polynomial division of a, b \in Z_p[x_0, \ldots, x_n]
32  * It doesn't check whether the inputs are proper polynomials, so be careful
33  * of what you pass in.
34  *  
35  * @param a  first multivariate polynomial (dividend)
36  * @param b  second multivariate polynomial (divisor)
37  * @param q  quotient (returned)
38  * @param var variables X iterator to first element of vector of symbols
39  *
40  * @return "true" when exact division succeeds (the quotient is returned in
41  *          q), "false" otherwise.
42  * @note @a p = 0 means the base ring is Z
43  */
44 bool divide_in_z_p(const ex &a, const ex &b, ex &q, const exvector& vars, const long p)
45 {
46         static const ex _ex1(1);
47         if (b.is_zero())
48                 throw(std::overflow_error("divide_in_z: division by zero"));
49         if (b.is_equal(_ex1)) {
50                 q = a;
51                 return true;
52         }
53         if (is_exactly_a<numeric>(a)) {
54                 if (is_exactly_a<numeric>(b)) {
55                         // p == 0 means division in Z
56                         if (p == 0) {
57                                 const numeric tmp = ex_to<numeric>(a/b);
58                                 if (tmp.is_integer()) {
59                                         q = tmp;
60                                         return true;
61                                 } else
62                                         return false;
63                         } else {
64                                 q = (a*recip(ex_to<numeric>(b), p)).smod(numeric(p));
65                                 return true;
66                         }
67                 } else
68                         return false;
69         }
70         if (a.is_equal(b)) {
71                 q = _ex1;
72                 return true;
73         }
74
75         // Main symbol
76         const ex &x = vars.back();
77
78         // Compare degrees
79         int adeg = a.degree(x), bdeg = b.degree(x);
80         if (bdeg > adeg)
81                 return false;
82
83         // Polynomial long division (recursive)
84         ex r = a.expand();
85         if (r.is_zero())
86                 return true;
87         int rdeg = adeg;
88         ex eb = b.expand();
89         ex blcoeff = eb.coeff(x, bdeg);
90         exvector v;
91         v.reserve(std::max(rdeg - bdeg + 1, 0));
92         exvector rest_vars(vars);
93         rest_vars.pop_back();
94         while (rdeg >= bdeg) {
95                 ex term, rcoeff = r.coeff(x, rdeg);
96                 if (!divide_in_z_p(rcoeff, blcoeff, term, rest_vars, p))
97                         break;
98                 term = (term*power(x, rdeg - bdeg)).expand();
99                 v.push_back(term);
100                 r = (r - term*eb).expand();
101                 if (p != 0)
102                         r = r.smod(numeric(p));
103                 if (r.is_zero()) {
104                         q = dynallocate<add>(v);
105                         return true;
106                 }
107                 rdeg = r.degree(x);
108         }
109         return false;
110 }
111
112 } // namespace GiNaC