Fixed bug in multivariate factorization. Fractional numerical content was not
authorJens Vollinga <jensv@balin.nikhef.nl>
Thu, 18 Dec 2008 10:10:50 +0000 (11:10 +0100)
committerJens Vollinga <jensv@balin.nikhef.nl>
Thu, 18 Dec 2008 10:10:50 +0000 (11:10 +0100)
extracted and led to a not so finite running behavior. Added new checks for that
bug.

Added static keywords to hide debugging output operators.

check/exam_factor.cpp
ginac/factor.cpp

index 8851b7c..71a63f2 100644 (file)
@@ -41,7 +41,7 @@ static unsigned check_factor(const ex& e)
 static unsigned exam_factor1()
 {
        unsigned result = 0;
-       ex e, d;
+       ex e;
        symbol x("x");
        lst syms;
        syms.append(x);
@@ -110,7 +110,7 @@ static unsigned exam_factor1()
 static unsigned exam_factor2()
 {
        unsigned result = 0;
-       ex e, d;
+       ex e;
        symbol x("x"), y("y"), z("z");
        lst syms;
        syms = x, y, z;
@@ -166,6 +166,23 @@ static unsigned exam_factor2()
        return result;
 }
 
+static unsigned exam_factor3()
+{
+       unsigned result = 0;
+       ex e;
+       symbol k("k"), n("n");
+       lst syms;
+       syms = k, n;
+       
+       e = ex("1/2*(-3+3*k-n)*(-2+3*k-n)*(-1+3*k-n)", syms);
+       result += check_factor(e);
+
+       e = ex("1/4*(2*k-n)*(-1+2*k-n)", syms);
+       result += check_factor(e);
+
+       return result;
+}
+
 unsigned exam_factor()
 {
        unsigned result = 0;
@@ -174,6 +191,7 @@ unsigned exam_factor()
 
        result += exam_factor1(); cout << '.' << flush;
        result += exam_factor2(); cout << '.' << flush;
+       result += exam_factor3(); cout << '.' << flush;
 
        return result;
 }
index b712913..da870d3 100644 (file)
@@ -92,7 +92,7 @@ ostream& operator<<(ostream& o, const vector<int>& v)
        }
        return o;
 }
-ostream& operator<<(ostream& o, const vector<cl_I>& v)
+static ostream& operator<<(ostream& o, const vector<cl_I>& v)
 {
        vector<cl_I>::const_iterator i = v.begin(), end = v.end();
        while ( i != end ) {
@@ -101,7 +101,7 @@ ostream& operator<<(ostream& o, const vector<cl_I>& v)
        }
        return o;
 }
-ostream& operator<<(ostream& o, const vector<cl_MI>& v)
+static ostream& operator<<(ostream& o, const vector<cl_MI>& v)
 {
        vector<cl_MI>::const_iterator i = v.begin(), end = v.end();
        while ( i != end ) {
@@ -2237,13 +2237,8 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
        const ex& x = *syms.begin();
 
        // make polynomial primitive
-       ex p = poly.collect(x);
-       ex cont = p.lcoeff(x);
-       for ( int i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
-               cont = gcd(cont, p.coeff(x,i));
-               if ( cont == 1 ) break;
-       }
-       ex pp = expand(normal(p / cont));
+       ex unit, cont, pp;
+       poly.unitcontprim(x, unit, cont, pp);
        if ( !is_a<numeric>(cont) ) {
                return factor_sqrfree(cont) * factor_sqrfree(pp);
        }
@@ -2427,7 +2422,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
                // try Hensel lifting
                ex res = hensel_multivar(pp, x, epv, prime, l, modfactors, C);
                if ( res != lst() ) {
-                       ex result = cont;
+                       ex result = cont * unit;
                        for ( size_t i=0; i<res.nops(); ++i ) {
                                result *= res.op(i).content(x) * res.op(i).unit(x);
                                result *= res.op(i).primpart(x);