- Complete revamp of methods in class matrix. Some redundant (and poor)

[ginac.git] / ginac / normal.cpp

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index a15ec62d4ef1434ef1c38294ab8be41b732cb184..514423c97930f4eff87ad74205553735069f3b26 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -1140,7 +1140,7 @@ static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
         c = d;
         cdeg = ddeg;
 //std::clog << " dividing...\n";
-        if (!divide_in_z(r, ri * pow(psi, delta), d, var+1))
+        if (!divide_in_z(r, ri * pow(psi, delta), d, var))
             throw(std::runtime_error("invalid expression in sr_gcd(), division failed"));
         ddeg = d.degree(x);
         if (ddeg == 0) {
@@ -1331,16 +1331,17 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
        heur_gcd_called++;
 #endif
 
+       // Algorithms only works for non-vanishing input polynomials
+       if (a.is_zero() || b.is_zero())
+               return *new ex(fail());
+
        // GCD of two numeric values -> 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));
-        numeric rg;
-        if (ca || cb)
-            rg = g.inverse();
         if (ca)
-            *ca = ex_to_numeric(a).mul(rg);
+            *ca = ex_to_numeric(a) / g;
         if (cb)
-            *cb = ex_to_numeric(b).mul(rg);
+            *cb = ex_to_numeric(b) / g;
         return g;
     }
 
@@ -1445,10 +1446,19 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
        // 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;
+               if (ca || cb) {
+                       if (g.is_zero()) {
+                               if (ca)
+                                       *ca = _ex0();
+                               if (cb)
+                                       *cb = _ex0();
+                       } else {
+                       if (ca)
+                           *ca = ex_to_numeric(a) / g;
+                       if (cb)
+                           *cb = ex_to_numeric(b) / g;
+                       }
+               }
         return g;
     }
 
@@ -1947,10 +1957,10 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
     // Add fractions
     if (den.is_equal(_ex1())) {
 
-               // Common denominator is 1, simply add all numerators
+               // Common denominator is 1, simply add all fractions
         exvector num_seq;
                for (ait=o.begin(); ait!=aitend; ait++) {
-                       num_seq.push_back(ait->op(0));
+                       num_seq.push_back(ait->op(0) / ait->op(1));
                }
                return (new lst((new add(num_seq))->setflag(status_flags::dynallocated), den))->setflag(status_flags::dynallocated);