From: Christian Bauer <Christian.Bauer@uni-mainz.de>
Date: Fri, 5 May 2000 17:09:04 +0000 (+0000)
Subject: - normal() respects the "level" parameter to limit the recursion depth
X-Git-Tag: release_0-6-0~9
X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=e7a89b57b362bf02a13c0e72292ff490924a550c;hp=e009cb7984c4971df3b9e036eaead640095f46d5

- normal() respects the "level" parameter to limit the recursion depth
- added sprem() function (sparse pseudo-remainder), as in Maple
---

diff --git a/ginac/inifcns.cpp b/ginac/inifcns.cpp
index f7d2864e..85a7365f 100644
--- a/ginac/inifcns.cpp
+++ b/ginac/inifcns.cpp
@@ -84,7 +84,7 @@ static ex csgn_eval(const ex & x)
     if (is_ex_exactly_of_type(x, numeric))
         return csgn(ex_to_numeric(x));
     
-    if (is_ex_exactly_of_type(x, mul)) {
+    else if (is_ex_exactly_of_type(x, mul)) {
         numeric oc = ex_to_numeric(x.op(x.nops()-1));
         if (oc.is_real()) {
             if (oc > 0)
@@ -102,8 +102,8 @@ static ex csgn_eval(const ex & x)
                 // csgn(-42*I*x) -> -csgn(I*x)
                 return -csgn(I*x/oc).hold();
         }
-    }
-    
+	}
+   
     return csgn(x).hold();
 }
 
diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index f93ecc86..074759df 100644
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -489,6 +489,57 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args)
 }
 
 
+/** Sparse pseudo-remainder of polynomials a(x) and b(x) in Z[x].
+ *
+ *  @param a  first polynomial in x (dividend)
+ *  @param b  second polynomial in x (divisor)
+ *  @param x  a and b are polynomials in x
+ *  @param check_args  check whether a and b are polynomials with rational
+ *         coefficients (defaults to "true")
+ *  @return sparse pseudo-remainder of a(x) and b(x) in Z[x] */
+
+ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args)
+{
+    if (b.is_zero())
+        throw(std::overflow_error("prem: division by zero"));
+    if (is_ex_exactly_of_type(a, numeric)) {
+        if (is_ex_exactly_of_type(b, numeric))
+            return _ex0();
+        else
+            return b;
+    }
+    if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)))
+        throw(std::invalid_argument("prem: arguments must be polynomials over the rationals"));
+
+    // Polynomial long division
+    ex r = a.expand();
+    ex eb = b.expand();
+    int rdeg = r.degree(x);
+    int bdeg = eb.degree(x);
+    ex blcoeff;
+    if (bdeg <= rdeg) {
+        blcoeff = eb.coeff(x, bdeg);
+        if (bdeg == 0)
+            eb = _ex0();
+        else
+            eb -= blcoeff * power(x, bdeg);
+    } else
+        blcoeff = _ex1();
+
+    while (rdeg >= bdeg && !r.is_zero()) {
+        ex rlcoeff = r.coeff(x, rdeg);
+        ex term = (power(x, rdeg - bdeg) * eb * rlcoeff).expand();
+        if (rdeg == 0)
+            r = _ex0();
+        else
+            r -= rlcoeff * power(x, rdeg);
+        r = (blcoeff * r).expand() - term;
+        rdeg = r.degree(x);
+    }
+    return r;
+}
+
+
 /** Exact polynomial division of a(X) by b(X) in Q[X].
  *  
  *  @param a  first multivariate polynomial (dividend)
@@ -1797,6 +1848,11 @@ static ex frac_cancel(const ex &n, const ex &d)
  *  @see ex::normal */
 ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
 {
+	if (level == 1)
+		return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated);
+	else if (level == -max_recursion_level)
+        throw(std::runtime_error("max recursion level reached"));
+
     // Normalize and expand children, chop into summands
     exvector o;
     o.reserve(seq.size()+1);
@@ -1873,6 +1929,11 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const
  *  @see ex::normal() */
 ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
 {
+	if (level == 1)
+		return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated);
+	else if (level == -max_recursion_level)
+        throw(std::runtime_error("max recursion level reached"));
+
     // Normalize children, separate into numerator and denominator
 	ex num = _ex1();
 	ex den = _ex1(); 
@@ -1899,6 +1960,11 @@ ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const
  *  @see ex::normal */
 ex power::normal(lst &sym_lst, lst &repl_lst, int level) const
 {
+	if (level == 1)
+		return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated);
+	else if (level == -max_recursion_level)
+        throw(std::runtime_error("max recursion level reached"));
+
 	// Normalize basis
     ex n = basis.bp->normal(sym_lst, repl_lst, level-1);