Fixed bug in the Q matrix calculation and the univariate factorization function.

diff --git a/check/exam_factor.cpp b/check/exam_factor.cpp
index ed064588e11bbabb42ff14829a946bfa42f89ec5..bf1751a572c02e7bf1cb13e0cad5d12514bbb78c 100644 (file)
--- a/check/exam_factor.cpp
+++ b/check/exam_factor.cpp
@@ -58,6 +58,9 @@ static unsigned exam_factor1()
        e = ex("(1+x)^3", syms);
        result += check_factor(e);
 
        e = ex("(1+x)^3", syms);
        result += check_factor(e);
 

+       e = ex("(x+1)*(x+4)", syms);
+       result += check_factor(e);
+

        e = ex("x^6-3*x^5+x^4-3*x^3-x^2-3*x+1", syms);
        result += check_factor(e);
 
        e = ex("x^6-3*x^5+x^4-3*x^3-x^2-3*x+1", syms);
        result += check_factor(e);
 


@@ -79,7 +82,11 @@ static unsigned exam_factor1()
        e = ex("x^8-40*x^6+352*x^4-960*x^2+576", syms);
        result += check_factor(e);
 
        e = ex("x^8-40*x^6+352*x^4-960*x^2+576", syms);
        result += check_factor(e);
 

-       e = ex("x*(2+x^2)*(1+x+x^3+x^2+x^6+x^5+x^4)*(1+x^3)^2*(-1+x)", syms);
+       e = ex("x*(2+x^2)*(1+x+x^3+x^2+x^6+x^5+x^4)*(1+x)^2*(1-x+x^2)^2*(-1+x)", syms);
+
+       result += check_factor(e);
+
+       e = ex("(x+4+x^2-x^3+43*x^4)*(x+1-x^2-3*x^3+4*x^4)", syms);

        result += check_factor(e);
 
        return result;
        result += check_factor(e);
 
        return result;

diff --git a/ginac/factor.cpp b/ginac/factor.cpp
index 18e7c04c78bc3ce86dab49e377f7a25de44d0aa3..5489b4ceb79eb87a99ba36ef9e17e5dac96b3005 100644 (file)
--- a/ginac/factor.cpp
+++ b/ginac/factor.cpp
@@ -700,19 +700,32 @@ static void q_matrix(const UniPoly& a, Matrix& Q)
 {
        unsigned int n = a.degree();
        unsigned int q = cl_I_to_uint(a.R->modulus);
 {
        unsigned int n = a.degree();
        unsigned int q = cl_I_to_uint(a.R->modulus);

-       vector<cl_MI> r(n, a.R->zero());
-       r[0] = a.R->one();
-       Q.set_row(0, r);
-       unsigned int max = (n-1) * q;
-       for ( size_t m=1; m<=max; ++m ) {
-               cl_MI rn_1 = r.back();
-               for ( size_t i=n-1; i>0; --i ) {
-                       r[i] = r[i-1] - rn_1 * a[i];
-               }
-               r[0] = -rn_1 * a[0];
-               if ( (m % q) == 0 ) {
-                       Q.set_row(m/q, r);
+// fast and buggy
+//     vector<cl_MI> r(n, a.R->zero());
+//     r[0] = a.R->one();
+//     Q.set_row(0, r);
+//     unsigned int max = (n-1) * q;
+//     for ( size_t m=1; m<=max; ++m ) {
+//             cl_MI rn_1 = r.back();
+//             for ( size_t i=n-1; i>0; --i ) {
+//                     r[i] = r[i-1] - rn_1 * a[i];
+//             }
+//             r[0] = -rn_1 * a[0];
+//             if ( (m % q) == 0 ) {
+//                     Q.set_row(m/q, r);
+//             }
+//     }
+// slow and (hopefully) correct
+       for ( size_t i=0; i<n; ++i ) {
+               UniPoly qk(a.R);
+               qk.set(i*q, a.R->one());
+               UniPoly r(a.R);
+               rem(qk, a, r);
+               Vec rvec;
+               for ( size_t j=0; j<n; ++j ) {
+                       rvec.push_back(r[j]);

                }
                }

+               Q.set_row(i, rvec);

        }
 }
 
        }
 }
 


@@ -795,7 +808,12 @@ static void berlekamp(const UniPoly& a, UniPolyVec& upv)
                                        *u = uo;
                                }
                                factors.push_back(g);
                                        *u = uo;
                                }
                                factors.push_back(g);

-                               ++size;
+                               size = 0;
+                               list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
+                               while ( i != end ) {
+                                       if ( i->degree() ) ++size; 
+                                       ++i;
+                               }

                                if ( size == k ) {
                                        list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
                                        while ( i != end ) {
                                if ( size == k ) {
                                        list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
                                        while ( i != end ) {


@@ -803,9 +821,9 @@ static void berlekamp(const UniPoly& a, UniPolyVec& upv)
                                        }
                                        return;
                                }
                                        }
                                        return;
                                }

-                               if ( u->degree() < nur.degree() ) {
-                                       break;
-                               }
+//                             if ( u->degree() < nur.degree() ) {
+//                                     break;
+//                             }

                        }
                }
                if ( ++r == k ) {
                        }
                }
                if ( ++r == k ) {


@@ -1107,6 +1125,7 @@ static ex factor_univariate(const ex& poly, const ex& x)
                                        mf.factors = newfactors2;
                                        mf.poly = answer.op(1);
                                        tocheck.push(mf);
                                        mf.factors = newfactors2;
                                        mf.poly = answer.op(1);
                                        tocheck.push(mf);

+                                       break;

                                }
                        }
                        else {
                                }
                        }
                        else {