G_do_hoelder: fix case with real x values which are not of type cl_R.

[ginac.git] / ginac / normal.cpp

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 9f8b7b4a406fdb5fe5ee39454e86eb6340f0d86a..ef1af87705563a79441929f91f919a9e52e2366a 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -212,10 +212,10 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v)
 
 #if 0
        std::clog << "Symbols:\n";
-       it = v.begin(); itend = v.end();
+       auto it = v.begin(), itend = v.end();
        while (it != itend) {
-               std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << endl;
-               std::clog << "  lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << endl;
+               std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << std::endl;
+               std::clog << "  lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << std::endl;
                ++it;
        }
 #endif
@@ -1470,7 +1470,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
        }
 
        // Some trivial cases
-       ex aex = a.expand(), bex = b.expand();
+       ex aex = a.expand();
        if (aex.is_zero()) {
                if (ca)
                        *ca = _ex0;
@@ -1478,6 +1478,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
                        *cb = _ex1;
                return b;
        }
+       ex bex = b.expand();
        if (bex.is_zero()) {
                if (ca)
                        *ca = _ex1;
@@ -1563,8 +1564,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
                        *cb = b;
                return _ex1;
        }
-       // move symbols which contained only in one of the polynomials
-       // to the end:
+       // move symbol contained only in one of the polynomials to the end:
        rotate(sym_stats.begin(), vari, sym_stats.end());
 
        sym_desc_vec::const_iterator var = sym_stats.begin();
@@ -1676,7 +1676,6 @@ static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb)
                        if (cb)
                                *cb = b;
                        return _ex1;
-                       // XXX: do I need to check for p_gcd = -1?
        }
 
        // there are common factors:
@@ -1706,17 +1705,27 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb)
        if (p.is_equal(b)) {
                // a = p^n, b = p, gcd = p
                if (ca)
-                       *ca = pow(p, a.op(1) - 1);
+                       *ca = pow(p, exp_a - 1);
                if (cb)
                        *cb = _ex1;
                return p;
-       } 
+       }
+       if (is_a<symbol>(p)) {
+               // Cancel trivial common factor
+               int ldeg_a = ex_to<numeric>(exp_a).to_int();
+               int ldeg_b = b.ldegree(p);
+               int min_ldeg = std::min(ldeg_a, ldeg_b);
+               if (min_ldeg > 0) {
+                       ex common = pow(p, min_ldeg);
+                       return gcd(pow(p, ldeg_a - min_ldeg), (b / common).expand(), ca, cb, false) * common;
+               }
+       }
 
        ex p_co, bpart_co;
        ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
 
-       // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
        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)
@@ -1802,7 +1811,7 @@ static epvector sqrfree_yun(const ex &a, const symbol &x)
        do {
                w = quo(w, g, x);
                if (w.is_zero()) {
-                       return res;
+                       return results;
                }
                z = quo(z, g, x) - w.diff(x);
                exponent = exponent + 1;