Remove obsolete debug #include file from factor.cpp.

[ginac.git] / ginac / factor.cpp

diff --git a/ginac/factor.cpp b/ginac/factor.cpp
index c978bab5574308015ec0af5e9b63feb267b1e847..dbe3ed3ea3bf2c494650cd13e3b0a190c6f77dc0 100644 (file)
--- a/ginac/factor.cpp
+++ b/ginac/factor.cpp
@@ -33,7 +33,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2017 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2020 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
@@ -66,7 +66,6 @@
 #include "add.h"
 
 #include <algorithm>
-#include <cmath>
 #include <limits>
 #include <list>
 #include <vector>
@@ -1492,6 +1491,9 @@ static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
        poly.unitcontprim(x, unit, cont, prim_ex);
        upoly prim;
        upoly_from_ex(prim, prim_ex, x);
+       if (prim_ex.is_equal(1)) {
+               return poly;
+       }
 
        // determine proper prime and minimize number of modular factors
        prime = 3;
@@ -1874,7 +1876,8 @@ static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex&
 {
        vector<ex> a = a_;
 
-       const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
+       const cl_I modulus = expt_pos(cl_I(p),k);
+       const cl_modint_ring R = find_modint_ring(modulus);
        const size_t r = a.size();
        const size_t nu = I.size() + 1;
 
@@ -1947,16 +1950,13 @@ static vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex&
                        cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
                        upvec delta_s = univar_diophant(amod, x, m, p, k);
                        cl_MI modcm;
-                       cl_I poscm = cm;
-                       while ( poscm < 0 ) {
-                               poscm = poscm + expt_pos(cl_I(p),k);
-                       }
+                       cl_I poscm = plusp(cm) ? cm : mod(cm, modulus);
                        modcm = cl_MI(R, poscm);
                        for ( size_t j=0; j<delta_s.size(); ++j ) {
                                delta_s[j] = delta_s[j] * modcm;
                                sigma[j] = sigma[j] + umodpoly_to_ex(delta_s[j], x);
                        }
-                       if ( nterms > 1 ) {
+                       if ( nterms > 1 && i+1 != nterms ) {
                                z = c.op(i+1);
                        }
                }
@@ -2232,14 +2232,13 @@ static ex factor_sqrfree(const ex& poly);
  */
 static ex factor_multivariate(const ex& poly, const exset& syms)
 {
-       exset::const_iterator s;
        const ex& x = *syms.begin();
 
        // make polynomial primitive
        ex unit, cont, pp;
        poly.unitcontprim(x, unit, cont, pp);
        if ( !is_a<numeric>(cont) ) {
-               return factor_sqrfree(cont) * factor_sqrfree(pp);
+               return unit * factor_sqrfree(cont) * factor_sqrfree(pp);
        }
 
        // factor leading coefficient
@@ -2312,7 +2311,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
                        vector<numeric> ftilde(vnlst.nops()-1);
                        for ( size_t i=0; i<ftilde.size(); ++i ) {
                                ex ft = vnlst.op(i+1);
-                               s = syms.begin();
+                               auto s = syms.begin();
                                ++s;
                                for ( size_t j=0; j<a.size(); ++j ) {
                                        ft = ft.subs(*s == a[j]);
@@ -2386,7 +2385,7 @@ static ex factor_multivariate(const ex& poly, const exset& syms)
                // set up evaluation points
                EvalPoint ep;
                vector<EvalPoint> epv;
-               s = syms.begin();
+               auto s = syms.begin();
                ++s;
                for ( size_t i=0; i<a.size(); ++i ) {
                        ep.x = *s++;
@@ -2458,9 +2457,9 @@ static ex factor_sqrfree(const ex& poly)
        if ( findsymbols.syms.size() == 1 ) {
                // univariate case
                const ex& x = *(findsymbols.syms.begin());
-               if ( poly.ldegree(x) > 0 ) {
+               int ld = poly.ldegree(x);
+               if ( ld > 0 ) {
                        // pull out direct factors
-                       int ld = poly.ldegree(x);
                        ex res = factor_univariate(expand(poly/pow(x, ld)), x);
                        return res * pow(x,ld);
                } else {
@@ -2593,7 +2592,3 @@ ex factor(const ex& poly, unsigned options)
 }
 
 } // namespace GiNaC
-
-#ifdef DEBUGFACTOR
-#include "test.h"
-#endif