- 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 f7d2864ef9ae64035eae269d140b6c36b6d41169..85a7365fa980ece03a596b9b9cff70724c8da529 100644 (file)
--- 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, 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)
         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();
         }
                 // csgn(-42*I*x) -> -csgn(I*x)
                 return -csgn(I*x/oc).hold();
         }

-    }
-    
+       }
+   

     return csgn(x).hold();
 }
 
     return csgn(x).hold();
 }
 

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index f93ecc8621b9e85e1082e9a100c366f9f8fa7ce8..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)
 /** 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
 {
  *  @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);
     // 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
 {
  *  @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(); 
     // 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
 {
  *  @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);
 
        // Normalize basis
     ex n = basis.bp->normal(sym_lst, repl_lst, level-1);