- normal() respects the "level" parameter to limit the recursion depth

[ginac.git] / ginac / normal.cpp

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 1126d27c2fffb738e3c6cf57decf92af6dc430f0..074759df4db92413e902b5ec6182af426dde48b6 100644 (file)
--- 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)
@@ -516,7 +567,8 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
         return true;
     }
 #endif
-    if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)))
+    if (check_args && (!a.info(info_flags::rational_polynomial) ||
+                       !b.info(info_flags::rational_polynomial)))
         throw(std::invalid_argument("divide: arguments must be polynomials over the rationals"));
 
     // Find first symbol
@@ -1796,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);
@@ -1872,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(); 
@@ -1898,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);
 
@@ -2046,7 +2113,8 @@ ex basic::to_rational(lst &repl_lst) const
 }
 
 
-/** Implementation of ex::to_rational() for symbols. This returns the unmodified symbol.
+/** Implementation of ex::to_rational() for symbols. This returns the
+ *  unmodified symbol.
  *  @see ex::to_rational */
 ex symbol::to_rational(lst &repl_lst) const
 {
@@ -2054,14 +2122,14 @@ ex symbol::to_rational(lst &repl_lst) const
 }
 
 
-/** Implementation of ex::to_rational() for a numeric. It splits complex numbers
- *  into re+I*im and replaces I and non-rational real numbers with a temporary
- *  symbol.
+/** Implementation of ex::to_rational() for a numeric. It splits complex
+ *  numbers into re+I*im and replaces I and non-rational real numbers with a
+ *  temporary symbol.
  *  @see ex::to_rational */
 ex numeric::to_rational(lst &repl_lst) const
 {
     if (is_real()) {
-        if (!is_integer())
+        if (!is_rational())
             return replace_with_symbol(*this, repl_lst);
     } else { // complex
         numeric re = real();