Happy New Year!

[ginac.git] / ginac / factor.cpp

diff --git a/ginac/factor.cpp b/ginac/factor.cpp
index 2c0ab8a5c07b04ab37aea5a78e4e22996fab9821..2b77d1b0110f3e327216a11a20510b901a6b93ee 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-2019 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>
@@ -1874,7 +1873,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 +1947,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);
                        }
                }
@@ -2500,13 +2497,41 @@ struct apply_factor_map : public map_function {
        }
 };
 
-} // anonymous namespace
+/** Iterate through explicit factors of e, call yield(f, k) for
+ *  each factor of the form f^k.
+ *
+ *  Note that this function doesn't factor e itself, it only
+ *  iterates through the factors already explicitly present.
+ */
+template <typename F> void
+factor_iter(const ex &e, F yield)
+{
+       if (is_a<mul>(e)) {
+               for (const auto &f : e) {
+                       if (is_a<power>(f)) {
+                               yield(f.op(0), f.op(1));
+                       } else {
+                               yield(f, ex(1));
+                       }
+               }
+       } else {
+               if (is_a<power>(e)) {
+                       yield(e.op(0), e.op(1));
+               } else {
+                       yield(e, ex(1));
+               }
+       }
+}
 
-/** Interface function to the outside world. It checks the arguments, tries a
- *  square free factorization, and then calls factor_sqrfree to do the hard
- *  work.
+/** This function factorizes a polynomial. It checks the arguments,
+ *  tries a square free factorization, and then calls factor_sqrfree
+ *  to do the hard work.
+ *
+ *  This function expands its argument, so for polynomials with
+ *  explicit factors it's better to call it on each one separately
+ *  (or use factor() which does just that).
  */
-ex factor(const ex& poly, unsigned options)
+static ex factor1(const ex& poly, unsigned options)
 {
        // check arguments
        if ( !poly.info(info_flags::polynomial) ) {
@@ -2533,44 +2558,35 @@ ex factor(const ex& poly, unsigned options)
        ex sfpoly = sqrfree(poly.expand(), syms);
 
        // factorize the square free components
-       if ( is_a<power>(sfpoly) ) {
-               // case: (polynomial)^exponent
-               const ex& base = sfpoly.op(0);
-               if ( !is_a<add>(base) ) {
-                       // simple case: (monomial)^exponent
-                       return sfpoly;
-               }
-               ex f = factor_sqrfree(base);
-               return pow(f, sfpoly.op(1));
-       }
-       if ( is_a<mul>(sfpoly) ) {
-               // case: multiple factors
-               ex res = 1;
-               for ( size_t i=0; i<sfpoly.nops(); ++i ) {
-                       const ex& t = sfpoly.op(i);
-                       if ( is_a<power>(t) ) {
-                               const ex& base = t.op(0);
-                               if ( !is_a<add>(base) ) {
-                                       res *= t;
-                               } else {
-                                       ex f = factor_sqrfree(base);
-                                       res *= pow(f, t.op(1));
-                               }
-                       } else if ( is_a<add>(t) ) {
-                               ex f = factor_sqrfree(t);
-                               res *= f;
+       ex res = 1;
+       factor_iter(sfpoly,
+               [&](const ex &f, const ex &k) {
+                       if ( is_a<add>(f) ) {
+                               res *= pow(factor_sqrfree(f), k);
                        } else {
-                               res *= t;
+                               // simple case: (monomial)^exponent
+                               res *= pow(f, k);
                        }
-               }
-               return res;
-       }
-       if ( is_a<symbol>(sfpoly) ) {
-               return poly;
-       }
-       // case: (polynomial)
-       ex f = factor_sqrfree(sfpoly);
-       return f;
+               });
+       return res;
+}
+
+} // anonymous namespace
+
+/** Interface function to the outside world. It uses factor1()
+ *  on each of the explicitly present factors of poly.
+ */
+ex factor(const ex& poly, unsigned options)
+{
+       ex result = 1;
+       factor_iter(poly,
+               [&](const ex &f1, const ex &k1) {
+                       factor_iter(factor1(f1, options),
+                               [&](const ex &f2, const ex &k2) {
+                                       result *= pow(f2, k1*k2);
+                               });
+               });
+       return result;
 }
 
 } // namespace GiNaC