Fixed compile errors introduced in previous commit
[ginac.git] / ginac / polynomial / mgcd.cpp
1 #include "operators.h"
2 #include "chinrem_gcd.h"
3 #include <cln/integer.h>
4 #include "pgcd.h"
5 #include "collect_vargs.h"
6 #include "primes_factory.h"
7 #include "divide_in_z_p.h"
8 #include "poly_cra.h"
9
10 namespace GiNaC
11 {
12
13 static cln::cl_I extract_integer_content(ex& Apr, const ex& A)
14 {
15         static const cln::cl_I n1(1);
16         const numeric icont_ = A.integer_content();
17         const cln::cl_I icont = cln::the<cln::cl_I>(icont_.to_cl_N());
18         if (icont != 1) {
19                 Apr = (A/icont_).expand();
20                 return icont;
21         } else {
22                 Apr = A;
23                 return n1;
24         }
25 }
26
27 ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
28 {
29         ex A, B;
30         const cln::cl_I a_icont = extract_integer_content(A, A_);
31         const cln::cl_I b_icont = extract_integer_content(B, B_);
32         const cln::cl_I c = cln::gcd(a_icont, b_icont);
33
34         const cln::cl_I a_lc = integer_lcoeff(A, vars);
35         const cln::cl_I b_lc = integer_lcoeff(B, vars);
36         const cln::cl_I g_lc = cln::gcd(a_lc, b_lc);
37
38         const ex& x(vars.back());
39         int n = std::min(A.degree(x), B.degree(x));
40         const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient()); 
41         const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient());
42         const cln::cl_I lcoeff_limit = (cln::cl_I(1) << n)*cln::abs(g_lc)*
43                 std::min(A_max_coeff, B_max_coeff);
44
45         cln::cl_I q = 0;
46         ex H = 0;
47
48         long p;
49         primes_factory pfactory;
50         while (true) {
51                 bool has_primes = pfactory(p, g_lc);
52                 if (!has_primes)
53                         throw chinrem_gcd_failed();
54
55                 const numeric pnum(p);
56                 ex Ap = A.smod(pnum);
57                 ex Bp = B.smod(pnum);
58                 ex Cp = pgcd(Ap, Bp, vars, p);
59
60                 const cln::cl_I g_lcp = smod(g_lc, p); 
61                 const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars);
62                 const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p);
63                 Cp = (Cp*numeric(nlc)).expand().smod(pnum);
64                 int cp_deg = Cp.degree(x);
65                 if (cp_deg == 0)
66                         return numeric(g_lc);
67                 if (zerop(q)) {
68                         H = Cp;
69                         n = cp_deg;
70                         q = p;
71                 } else {
72                         if (cp_deg == n) {
73                                 ex H_next = chinese_remainder(H, q, Cp, p);
74                                 q = q*cln::cl_I(p);
75                                 H = H_next;
76                         } else if (cp_deg < n) {
77                                 // all previous homomorphisms are unlucky
78                                 q = p;
79                                 H = Cp;
80                                 n = cp_deg;
81                         } else {
82                                 // dp_deg > d_deg: current prime is bad
83                         }
84                 }
85                 if (q < lcoeff_limit)
86                         continue; // don't bother to do division checks
87                 ex C, dummy1, dummy2;
88                 extract_integer_content(C, H);
89                 if (divide_in_z_p(A, C, dummy1, vars, 0) && 
90                                 divide_in_z_p(B, C, dummy2, vars, 0))
91                         return (numeric(c)*C).expand();
92                 // else: try more primes
93         }
94 }
95
96 } // namespace GiNaC
97