+//clog << "gcd(" << a << "," << b << ")\n";
+#if STATISTICS
+ gcd_called++;
+#endif
+
+ // GCD of numerics -> CLN
+ 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));
+ if (ca)
+ *ca = ex_to_numeric(a) / g;
+ if (cb)
+ *cb = ex_to_numeric(b) / g;
+ return g;
+ }
+
+ // Check arguments
+ if (check_args && !a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) {
+ throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals"));
+ }
+
+ // Partially factored cases (to avoid expanding large expressions)
+ if (is_ex_exactly_of_type(a, mul)) {
+ if (is_ex_exactly_of_type(b, mul) && b.nops() > a.nops())
+ goto factored_b;
+factored_a:
+ ex g = _ex1();
+ ex acc_ca = _ex1();
+ ex part_b = b;
+ for (unsigned i=0; i<a.nops(); i++) {
+ ex part_ca, part_cb;
+ g *= gcd(a.op(i), part_b, &part_ca, &part_cb, check_args);
+ acc_ca *= part_ca;
+ part_b = part_cb;
+ }
+ if (ca)
+ *ca = acc_ca;
+ if (cb)
+ *cb = part_b;
+ return g;
+ } else if (is_ex_exactly_of_type(b, mul)) {
+ if (is_ex_exactly_of_type(a, mul) && a.nops() > b.nops())
+ goto factored_a;
+factored_b:
+ ex g = _ex1();
+ ex acc_cb = _ex1();
+ ex part_a = a;
+ for (unsigned i=0; i<b.nops(); i++) {
+ ex part_ca, part_cb;
+ g *= gcd(part_a, b.op(i), &part_ca, &part_cb, check_args);
+ acc_cb *= part_cb;
+ part_a = part_ca;
+ }
+ if (ca)
+ *ca = part_a;
+ if (cb)
+ *cb = acc_cb;
+ return g;
+ }
+
+#if FAST_COMPARE
+ // Input polynomials of the form poly^n are sometimes also trivial
+ if (is_ex_exactly_of_type(a, power)) {
+ ex p = a.op(0);
+ if (is_ex_exactly_of_type(b, power)) {
+ if (p.is_equal(b.op(0))) {
+ // a = p^n, b = p^m, gcd = p^min(n, m)
+ ex exp_a = a.op(1), exp_b = b.op(1);
+ if (exp_a < exp_b) {
+ if (ca)
+ *ca = _ex1();
+ if (cb)
+ *cb = power(p, exp_b - exp_a);
+ return power(p, exp_a);
+ } else {
+ if (ca)
+ *ca = power(p, exp_a - exp_b);
+ if (cb)
+ *cb = _ex1();
+ return power(p, exp_b);
+ }
+ }
+ } else {
+ if (p.is_equal(b)) {
+ // a = p^n, b = p, gcd = p
+ if (ca)
+ *ca = power(p, a.op(1) - 1);
+ if (cb)
+ *cb = _ex1();
+ return p;
+ }
+ }
+ } else if (is_ex_exactly_of_type(b, power)) {
+ ex p = b.op(0);
+ if (p.is_equal(a)) {
+ // a = p, b = p^n, gcd = p
+ if (ca)
+ *ca = _ex1();
+ if (cb)
+ *cb = power(p, b.op(1) - 1);
+ return p;
+ }
+ }
+#endif
+