]> www.ginac.de Git - ginac.git/blobdiff - ginac/normal.cpp
- is_zero() is now called on expanded expressions in gcd()
[ginac.git] / ginac / normal.cpp
index 614b1db1bfc4a00d5630ecfeeadcbee69fdc5956..70447acfc5cb94c793c579af7b39abb5e0913d01 100644 (file)
@@ -6,9 +6,47 @@
  *  computation, square-free factorization and rational function normalization.
  */
 
+/*
+ *  GiNaC Copyright (C) 1999 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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ */
+
 #include <stdexcept>
+#include <algorithm>
+#include <map>
 
-#include "ginac.h"
+#include "normal.h"
+#include "basic.h"
+#include "ex.h"
+#include "add.h"
+#include "constant.h"
+#include "expairseq.h"
+#include "fail.h"
+#include "indexed.h"
+#include "inifcns.h"
+#include "lst.h"
+#include "mul.h"
+#include "ncmul.h"
+#include "numeric.h"
+#include "power.h"
+#include "relational.h"
+#include "series.h"
+#include "symbol.h"
+
+namespace GiNaC {
 
 // If comparing expressions (ex::compare()) is fast, you can set this to 1.
 // Some routines like quo(), rem() and gcd() will then return a quick answer
@@ -48,8 +86,6 @@ static bool get_first_symbol(const ex &e, const symbol *&x)
  *  Statistical information about symbols in polynomials
  */
 
-#include <algorithm>
-
 /** This structure holds information about the highest and lowest degrees
  *  in which a symbol appears in two multivariate polynomials "a" and "b".
  *  A vector of these structures with information about all symbols in
@@ -455,8 +491,6 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
  *  Remembering
  */
 
-#include <map>
-
 typedef pair<ex, ex> ex2;
 typedef pair<ex, bool> exbool;
 
@@ -874,7 +908,7 @@ ex basic::smod(const numeric &xi) const
 
 ex numeric::smod(const numeric &xi) const
 {
-    return ::smod(*this, xi);
+    return GiNaC::smod(*this, xi);
 }
 
 ex add::smod(const numeric &xi) const
@@ -885,13 +919,13 @@ ex add::smod(const numeric &xi) const
     epvector::const_iterator itend = seq.end();
     while (it != itend) {
         ASSERT(!is_ex_exactly_of_type(it->rest,numeric));
-        numeric coeff = ::smod(ex_to_numeric(it->coeff), xi);
+        numeric coeff = GiNaC::smod(ex_to_numeric(it->coeff), xi);
         if (!coeff.is_zero())
             newseq.push_back(expair(it->rest, coeff));
         it++;
     }
     ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
-    numeric coeff = ::smod(ex_to_numeric(overall_coeff), xi);
+    numeric coeff = GiNaC::smod(ex_to_numeric(overall_coeff), xi);
     return (new add(newseq,coeff))->setflag(status_flags::dynallocated);
 }
 
@@ -907,7 +941,7 @@ ex mul::smod(const numeric &xi) const
 #endif // def DOASSERT
     mul * mulcopyp=new mul(*this);
     ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
-    mulcopyp->overall_coeff=::smod(ex_to_numeric(overall_coeff),xi);
+    mulcopyp->overall_coeff = GiNaC::smod(ex_to_numeric(overall_coeff),xi);
     mulcopyp->clearflag(status_flags::evaluated);
     mulcopyp->clearflag(status_flags::hash_calculated);
     return mulcopyp->setflag(status_flags::dynallocated);
@@ -1016,21 +1050,22 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const
 ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
 {
     // Some trivial cases
-    if (a.is_zero()) {
+       ex aex = a.expand(), bex = b.expand();
+    if (aex.is_zero()) {
         if (ca)
             *ca = exZERO();
         if (cb)
             *cb = exONE();
         return b;
     }
-    if (b.is_zero()) {
+    if (bex.is_zero()) {
         if (ca)
             *ca = exONE();
         if (cb)
             *cb = exZERO();
         return a;
     }
-    if (a.is_equal(exONE()) || b.is_equal(exONE())) {
+    if (aex.is_equal(exONE()) || bex.is_equal(exONE())) {
         if (ca)
             *ca = a;
         if (cb)
@@ -1046,17 +1081,15 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
         return a;
     }
 #endif
-    if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) {
-        numeric g = gcd(ex_to_numeric(a), ex_to_numeric(b));
+    if (is_ex_exactly_of_type(aex, numeric) && is_ex_exactly_of_type(bex, numeric)) {
+        numeric g = gcd(ex_to_numeric(aex), ex_to_numeric(bex));
         if (ca)
-            *ca = ex_to_numeric(a) / g;
+            *ca = ex_to_numeric(aex) / g;
         if (cb)
-            *cb = ex_to_numeric(b) / g;
+            *cb = ex_to_numeric(bex) / g;
         return g;
     }
     if (check_args && !a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) {
-        cerr << "a=" << a << endl;
-        cerr << "b=" << b << endl;
         throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals"));
     }
 
@@ -1075,40 +1108,40 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
     if (min_ldeg > 0) {
         ex common = power(*x, min_ldeg);
 //clog << "trivial common factor " << common << endl;
-        return gcd((a / common).expand(), (b / common).expand(), ca, cb, false) * common;
+        return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common;
     }
 
     // Try to eliminate variables
     if (var->deg_a == 0) {
 //clog << "eliminating variable " << *x << " from b" << endl;
-        ex c = b.content(*x);
-        ex g = gcd(a, c, ca, cb, false);
+        ex c = bex.content(*x);
+        ex g = gcd(aex, c, ca, cb, false);
         if (cb)
-            *cb *= b.unit(*x) * b.primpart(*x, c);
+            *cb *= bex.unit(*x) * bex.primpart(*x, c);
         return g;
     } else if (var->deg_b == 0) {
 //clog << "eliminating variable " << *x << " from a" << endl;
-        ex c = a.content(*x);
-        ex g = gcd(c, b, ca, cb, false);
+        ex c = aex.content(*x);
+        ex g = gcd(c, bex, ca, cb, false);
         if (ca)
-            *ca *= a.unit(*x) * a.primpart(*x, c);
+            *ca *= aex.unit(*x) * aex.primpart(*x, c);
         return g;
     }
 
     // Try heuristic algorithm first, fall back to PRS if that failed
     ex g;
     try {
-        g = heur_gcd(a.expand(), b.expand(), ca, cb, var);
+        g = heur_gcd(aex, bex, ca, cb, var);
     } catch (gcdheu_failed) {
         g = *new ex(fail());
     }
     if (is_ex_exactly_of_type(g, fail)) {
 //clog << "heuristics failed\n";
-        g = sr_gcd(a, b, x);
+        g = sr_gcd(aex, bex, x);
         if (ca)
-            divide(a, g, *ca, false);
+            divide(aex, g, *ca, false);
         if (cb)
-            divide(b, g, *cb, false);
+            divide(bex, g, *cb, false);
     }
     return g;
 }
@@ -1450,3 +1483,5 @@ ex ex::normal(int level) const
     else
         return e;
 }
+
+} // namespace GiNaC