Improve normalization with nested functions.

[ginac.git] / ginac / normal.cpp

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 09773d37489693a4854964f129002abb186b2c5b..f22bf27b131ec1952b6f0f4192da37b842e51ba9 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -6,7 +6,7 @@
  *  computation, square-free factorization and rational function normalization. */
 
 /*
- *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2015 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
@@ -23,9 +23,6 @@
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#include <algorithm>
-#include <map>
-
 #include "normal.h"
 #include "basic.h"
 #include "ex.h"
@@ -44,6 +41,10 @@
 #include "pseries.h"
 #include "symbol.h"
 #include "utils.h"
+#include "polynomial/chinrem_gcd.h"
+
+#include <algorithm>
+#include <map>
 
 namespace GiNaC {
 
@@ -672,12 +673,13 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
                        q = rem_i*power(ab, a_exp - 1);
                        return true;
                }
-               for (int i=2; i < a_exp; i++) {
-                       if (divide(power(ab, i), b, rem_i, false)) {
-                               q = rem_i*power(ab, a_exp - i);
-                               return true;
-                       }
-               } // ... so we *really* need to expand expression.
+// code below is commented-out because it leads to a significant slowdown
+//             for (int i=2; i < a_exp; i++) {
+//                     if (divide(power(ab, i), b, rem_i, false)) {
+//                             q = rem_i*power(ab, a_exp - i);
+//                             return true;
+//                     }
+//             } // ... so we *really* need to expand expression.
        }
        
        // Polynomial long division (recursive)
@@ -1622,8 +1624,15 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio
                }
 #endif
        }
+       if (options & gcd_options::use_sr_gcd) {
+               g = sr_gcd(aex, bex, var);
+       } else {
+               exvector vars;
+               for (std::size_t n = sym_stats.size(); n-- != 0; )
+                       vars.push_back(sym_stats[n].sym);
+               g = chinrem_gcd(aex, bex, vars);
+       }
 
-       g = sr_gcd(aex, bex, var);
        if (g.is_equal(_ex1)) {
                // Keep cofactors factored if possible
                if (ca)
@@ -1997,16 +2006,18 @@ ex sqrfree_parfrac(const ex & a, const symbol & x)
  *  @see ex::normal */
 static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup)
 {
+       // Since the repl contains replaced expressions we should search for them
+       ex e_replaced = e.subs(repl, subs_options::no_pattern);
+
        // Expression already replaced? Then return the assigned symbol
-       exmap::const_iterator it = rev_lookup.find(e);
+       exmap::const_iterator it = rev_lookup.find(e_replaced);
        if (it != rev_lookup.end())
                return it->second;
-       
+
        // Otherwise create new symbol and add to list, taking care that the
        // replacement expression doesn't itself contain symbols from repl,
        // because subs() is not recursive
        ex es = (new symbol)->setflag(status_flags::dynallocated);
-       ex e_replaced = e.subs(repl, subs_options::no_pattern);
        repl.insert(std::make_pair(es, e_replaced));
        rev_lookup.insert(std::make_pair(e_replaced, es));
        return es;
@@ -2019,16 +2030,18 @@ static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup)
  *  @see basic::to_polynomial */
 static ex replace_with_symbol(const ex & e, exmap & repl)
 {
+       // Since the repl contains replaced expressions we should search for them
+       ex e_replaced = e.subs(repl, subs_options::no_pattern);
+
        // Expression already replaced? Then return the assigned symbol
        for (exmap::const_iterator it = repl.begin(); it != repl.end(); ++it)
-               if (it->second.is_equal(e))
+               if (it->second.is_equal(e_replaced))
                        return it->first;
-       
+
        // Otherwise create new symbol and add to list, taking care that the
        // replacement expression doesn't itself contain symbols from repl,
        // because subs() is not recursive
        ex es = (new symbol)->setflag(status_flags::dynallocated);
-       ex e_replaced = e.subs(repl, subs_options::no_pattern);
        repl.insert(std::make_pair(es, e_replaced));
        return es;
 }