-static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb, bool check_args)
-{
- if (is_exactly_a<power>(a)) {
- ex p = a.op(0);
- const ex& exp_a = a.op(1);
- if (is_exactly_a<power>(b)) {
- ex pb = b.op(0);
- const ex& exp_b = b.op(1);
- if (p.is_equal(pb)) {
- // a = p^n, b = p^m, gcd = p^min(n, m)
- 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 {
- ex p_co, pb_co;
- ex p_gcd = gcd(p, pb, &p_co, &pb_co, check_args);
- if (p_gcd.is_equal(_ex1)) {
- // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==>
- // gcd(a,b) = 1
- if (ca)
- *ca = a;
- if (cb)
- *cb = b;
- return _ex1;
- // XXX: do I need to check for p_gcd = -1?
- } else {
- // there are common factors:
- // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
- // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
- if (exp_a < exp_b) {
- return power(p_gcd, exp_a)*
- gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false);
- } else {
- return power(p_gcd, exp_b)*
- gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false);
- }
- } // p_gcd.is_equal(_ex1)
- } // p.is_equal(pb)
-
- } 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;
- }
-
- ex p_co, bpart_co;
- ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
-
- if (p_gcd.is_equal(_ex1)) {
- // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
- if (ca)
- *ca = a;
- if (cb)
- *cb = b;
- return _ex1;
- } else {
- // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
- return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false);
- }
- } // is_exactly_a<power>(b)