Update copyright notice.
[ginac.git] / ginac / polynomial / mgcd.cpp
1 /** @file mgcd.cpp
2  *
3  *  Chinese remainder algorithm. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2010 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 "operators.h"
24 #include "chinrem_gcd.h"
25 #include "pgcd.h"
26 #include "collect_vargs.h"
27 #include "primes_factory.h"
28 #include "divide_in_z_p.h"
29 #include "poly_cra.h"
30
31 #include <cln/integer.h>
32
33 namespace GiNaC {
34
35 static cln::cl_I extract_integer_content(ex& Apr, const ex& A)
36 {
37         static const cln::cl_I n1(1);
38         const numeric icont_ = A.integer_content();
39         const cln::cl_I icont = cln::the<cln::cl_I>(icont_.to_cl_N());
40         if (icont != 1) {
41                 Apr = (A/icont_).expand();
42                 return icont;
43         } else {
44                 Apr = A;
45                 return n1;
46         }
47 }
48
49 ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
50 {
51         ex A, B;
52         const cln::cl_I a_icont = extract_integer_content(A, A_);
53         const cln::cl_I b_icont = extract_integer_content(B, B_);
54         const cln::cl_I c = cln::gcd(a_icont, b_icont);
55
56         const cln::cl_I a_lc = integer_lcoeff(A, vars);
57         const cln::cl_I b_lc = integer_lcoeff(B, vars);
58         const cln::cl_I g_lc = cln::gcd(a_lc, b_lc);
59
60         const ex& x(vars.back());
61         int n = std::min(A.degree(x), B.degree(x));
62         const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient()); 
63         const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient());
64         const cln::cl_I lcoeff_limit = (cln::cl_I(1) << n)*cln::abs(g_lc)*
65                 std::min(A_max_coeff, B_max_coeff);
66
67         cln::cl_I q = 0;
68         ex H = 0;
69
70         long p;
71         primes_factory pfactory;
72         while (true) {
73                 bool has_primes = pfactory(p, g_lc);
74                 if (!has_primes)
75                         throw chinrem_gcd_failed();
76
77                 const numeric pnum(p);
78                 ex Ap = A.smod(pnum);
79                 ex Bp = B.smod(pnum);
80                 ex Cp = pgcd(Ap, Bp, vars, p);
81
82                 const cln::cl_I g_lcp = smod(g_lc, p); 
83                 const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars);
84                 const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p);
85                 Cp = (Cp*numeric(nlc)).expand().smod(pnum);
86                 int cp_deg = Cp.degree(x);
87                 if (cp_deg == 0)
88                         return numeric(g_lc);
89                 if (zerop(q)) {
90                         H = Cp;
91                         n = cp_deg;
92                         q = p;
93                 } else {
94                         if (cp_deg == n) {
95                                 ex H_next = chinese_remainder(H, q, Cp, p);
96                                 q = q*cln::cl_I(p);
97                                 H = H_next;
98                         } else if (cp_deg < n) {
99                                 // all previous homomorphisms are unlucky
100                                 q = p;
101                                 H = Cp;
102                                 n = cp_deg;
103                         } else {
104                                 // dp_deg > d_deg: current prime is bad
105                         }
106                 }
107                 if (q < lcoeff_limit)
108                         continue; // don't bother to do division checks
109                 ex C, dummy1, dummy2;
110                 extract_integer_content(C, H);
111                 if (divide_in_z_p(A, C, dummy1, vars, 0) && 
112                                 divide_in_z_p(B, C, dummy2, vars, 0))
113                         return (numeric(c)*C).expand();
114                 // else: try more primes
115         }
116 }
117
118 } // namespace GiNaC