[bugfix] chinrem_gcd: handle polynomials over rationals properly.

[ginac.git] / ginac / polynomial / mgcd.cpp

diff --git a/ginac/polynomial/mgcd.cpp b/ginac/polynomial/mgcd.cpp
index b39475cc5fa05e88b7ea23a5911966d391c64e62..5b9da128223d8a6d9024b24e775398b2e3ed8b53 100644 (file)
--- a/ginac/polynomial/mgcd.cpp
+++ b/ginac/polynomial/mgcd.cpp
@@ -1,18 +1,52 @@
+/** @file mgcd.cpp
+ *
+ *  Chinese remainder algorithm. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2010 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#include "operators.h"
 #include "chinrem_gcd.h"
-#include <cln/integer.h>
 #include "pgcd.h"
 #include "collect_vargs.h"
 #include "primes_factory.h"
 #include "divide_in_z_p.h"
 #include "poly_cra.h"
+#include <numeric> // std::accumulate
 
-namespace GiNaC
-{
+#include <cln/integer.h>
+#include <cln/integer_ring.h>
+#include <cln/rational.h>
+#include <cln/rational_ring.h>
+
+namespace GiNaC {
 
 static cln::cl_I extract_integer_content(ex& Apr, const ex& A)
 {
        static const cln::cl_I n1(1);
        const numeric icont_ = A.integer_content();
+       if (cln::instanceof(icont_.to_cl_N(), cln::cl_RA_ring)) {
+               Apr = (A/icont_).expand();
+               // A is a polynomail over rationals, so GCD is defined
+               // up to arbitrary rational number.
+               return n1;
+       }
+       GINAC_ASSERT(cln::instanceof(icont_.to_cl_N(), cln::cl_I_ring));
        const cln::cl_I icont = cln::the<cln::cl_I>(icont_.to_cl_N());
        if (icont != 1) {
                Apr = (A/icont_).expand();
@@ -35,12 +69,15 @@ ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
        const cln::cl_I g_lc = cln::gcd(a_lc, b_lc);
 
        const ex& x(vars.back());
-       int n = std::min(A.degree(x), B.degree(x));
+       exp_vector_t n = std::min(degree_vector(A, vars), degree_vector(B, vars));
+       const int nTot = std::accumulate(n.begin(), n.end(), 0);
        const cln::cl_I A_max_coeff = to_cl_I(A.max_coefficient()); 
        const cln::cl_I B_max_coeff = to_cl_I(B.max_coefficient());
-       const cln::cl_I lcoeff_limit = (cln::cl_I(1) << n)*cln::abs(g_lc)*
+
+       const cln::cl_I lcoeff_limit = (cln::cl_I(1) << nTot)*cln::abs(g_lc)*
                std::min(A_max_coeff, B_max_coeff);
 
+
        cln::cl_I q = 0;
        ex H = 0;
 
@@ -60,8 +97,8 @@ ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
                const cln::cl_I Cp_lc = integer_lcoeff(Cp, vars);
                const cln::cl_I nlc = smod(recip(Cp_lc, p)*g_lcp, p);
                Cp = (Cp*numeric(nlc)).expand().smod(pnum);
-               int cp_deg = Cp.degree(x);
-               if (cp_deg == 0)
+               exp_vector_t cp_deg = degree_vector(Cp, vars);
+               if (zerop(cp_deg))
                        return numeric(g_lc);
                if (zerop(q)) {
                        H = Cp;
@@ -93,4 +130,3 @@ ex chinrem_gcd(const ex& A_, const ex& B_, const exvector& vars)
 }
 
 } // namespace GiNaC
-