- if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) {
- numeric g = gcd(ex_to_numeric(a), ex_to_numeric(b));
- numeric rg;
- if (ca || cb)
- rg = g.inverse();
- if (ca)
- *ca = ex_to_numeric(a).mul(rg);
- if (cb)
- *cb = ex_to_numeric(b).mul(rg);
- return g;
- }
-
- // The first symbol is our main variable
- const symbol *x = var->sym;
-
- // Remove integer content
- numeric gc = gcd(a.integer_content(), b.integer_content());
- numeric rgc = gc.inverse();
- ex p = a * rgc;
- ex q = b * rgc;
- int maxdeg = max(p.degree(*x), q.degree(*x));
-
- // Find evaluation point
- numeric mp = p.max_coefficient(), mq = q.max_coefficient();
- numeric xi;
- if (mp > mq)
- xi = mq * _num2() + _num2();
- else
- xi = mp * _num2() + _num2();
-
- // 6 tries maximum
- for (int t=0; t<6; t++) {
- if (xi.int_length() * maxdeg > 100000) {
-//clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << endl;
- throw gcdheu_failed();
+ if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
+ numeric g = gcd(ex_to<numeric>(a), ex_to<numeric>(b));
+ if (ca)
+ *ca = ex_to<numeric>(a) / g;
+ if (cb)
+ *cb = ex_to<numeric>(b) / g;
+ return g;
+ }
+
+ // The first symbol is our main variable
+ const symbol &x = *(var->sym);
+
+ // Remove integer content
+ numeric gc = gcd(a.integer_content(), b.integer_content());
+ numeric rgc = gc.inverse();
+ ex p = a * rgc;
+ ex q = b * rgc;
+ int maxdeg = std::max(p.degree(x), q.degree(x));
+
+ // Find evaluation point
+ numeric mp = p.max_coefficient();
+ numeric mq = q.max_coefficient();
+ numeric xi;
+ if (mp > mq)
+ xi = mq * _num2 + _num2;
+ else
+ xi = mp * _num2 + _num2;
+
+ // 6 tries maximum
+ for (int t=0; t<6; t++) {
+ if (xi.int_length() * maxdeg > 100000) {
+//std::clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << std::endl;
+ throw gcdheu_failed();
+ }
+
+ // Apply evaluation homomorphism and calculate GCD
+ ex cp, cq;
+ ex gamma = heur_gcd(p.subs(x == xi), q.subs(x == xi), &cp, &cq, var+1).expand();
+ if (!is_exactly_a<fail>(gamma)) {
+
+ // Reconstruct polynomial from GCD of mapped polynomials
+ ex g = interpolate(gamma, xi, x, maxdeg);
+
+ // Remove integer content
+ g /= g.integer_content();
+
+ // If the calculated polynomial divides both p and q, this is the GCD
+ ex dummy;
+ if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) {
+ g *= gc;
+ ex lc = g.lcoeff(x);
+ if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
+ return -g;
+ else
+ return g;
+ }
+#if 0
+ cp = interpolate(cp, xi, x);
+ if (divide_in_z(cp, p, g, var)) {
+ if (divide_in_z(g, q, cb ? *cb : dummy, var)) {
+ g *= gc;
+ if (ca)
+ *ca = cp;
+ ex lc = g.lcoeff(x);
+ if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
+ return -g;
+ else
+ return g;
+ }
+ }
+ cq = interpolate(cq, xi, x);
+ if (divide_in_z(cq, q, g, var)) {
+ if (divide_in_z(g, p, ca ? *ca : dummy, var)) {
+ g *= gc;
+ if (cb)
+ *cb = cq;
+ ex lc = g.lcoeff(x);
+ if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
+ return -g;
+ else
+ return g;
+ }
+ }
+#endif