fixed a bug where quo() would call vector::reserve() with a negative argument

[ginac.git] / ginac / normal.cpp

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 4e2feb39a546957438891578486e0b07a2479b7a..4b490ed0e9824d301cb979fc830433b98aa22019 100644 (file)
--- a/ginac/normal.cpp
+++ b/ginac/normal.cpp
@@ -381,7 +381,7 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args)
        int rdeg = r.degree(x);
        ex blcoeff = b.expand().coeff(x, bdeg);
        bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric);
-       exvector v; v.reserve(rdeg - bdeg + 1);
+       exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
        while (rdeg >= bdeg) {
                ex term, rcoeff = r.coeff(x, rdeg);
                if (blcoeff_is_numeric)
@@ -581,14 +581,15 @@ ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args)
  *  @param check_args  check whether a and b are polynomials with rational
  *         coefficients (defaults to "true")
  *  @return "true" when exact division succeeds (quotient returned in q),
- *          "false" otherwise */
+ *          "false" otherwise (q left untouched) */
 bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 {
-       q = _ex0;
        if (b.is_zero())
                throw(std::overflow_error("divide: division by zero"));
-       if (a.is_zero())
+       if (a.is_zero()) {
+               q = _ex0;
                return true;
+       }
        if (is_ex_exactly_of_type(b, numeric)) {
                q = a / b;
                return true;
@@ -611,13 +612,15 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args)
 
        // Polynomial long division (recursive)
        ex r = a.expand();
-       if (r.is_zero())
+       if (r.is_zero()) {
+               q = _ex0;
                return true;
+       }
        int bdeg = b.degree(*x);
        int rdeg = r.degree(*x);
        ex blcoeff = b.expand().coeff(*x, bdeg);
        bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric);
-       exvector v; v.reserve(rdeg - bdeg + 1);
+       exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
        while (rdeg >= bdeg) {
                ex term, rcoeff = r.coeff(*x, rdeg);
                if (blcoeff_is_numeric)
@@ -778,7 +781,7 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
        int rdeg = adeg;
        ex eb = b.expand();
        ex blcoeff = eb.coeff(*x, bdeg);
-       exvector v; v.reserve(rdeg - bdeg + 1);
+       exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
        while (rdeg >= bdeg) {
                ex term, rcoeff = r.coeff(*x, rdeg);
                if (!divide_in_z(rcoeff, blcoeff, term, var+1))
@@ -1742,11 +1745,11 @@ static exvector sqrfree_yun(const ex &a, const symbol &x)
  *
  *  @param a  multivariate polynomial over Q[X]
  *  @param x  lst of variables to factor in, may be left empty for autodetection
- *  @return   polynomail a in square-free factored form. */
+ *  @return   polynomial a in square-free factored form. */
 ex sqrfree(const ex &a, const lst &l)
 {
-       if (is_ex_of_type(a,numeric) ||     // algorithm does not trap a==0
-           is_ex_of_type(a,symbol))        // shortcut
+       if (is_a<numeric>(a) ||     // algorithm does not trap a==0
+           is_a<symbol>(a))        // shortcut
                return a;
 
        // If no lst of variables to factorize in was specified we have to
@@ -1768,11 +1771,11 @@ ex sqrfree(const ex &a, const lst &l)
        // Find the symbol to factor in at this stage
        if (!is_ex_of_type(args.op(0), symbol))
                throw (std::runtime_error("sqrfree(): invalid factorization variable"));
-       const symbol x = ex_to<symbol>(args.op(0));
+       const symbol &x = ex_to<symbol>(args.op(0));
 
        // convert the argument from something in Q[X] to something in Z[X]
-       numeric lcm = lcm_of_coefficients_denominators(a);
-       ex tmp = multiply_lcm(a,lcm);
+       const numeric lcm = lcm_of_coefficients_denominators(a);
+       const ex tmp = multiply_lcm(a,lcm);
 
        // find the factors
        exvector factors = sqrfree_yun(tmp,x);
@@ -1781,10 +1784,10 @@ ex sqrfree(const ex &a, const lst &l)
        lst newargs = args;
        newargs.remove_first();
 
-       // recurse down the factors in remaining vars
+       // recurse down the factors in remaining variables
        if (newargs.nops()>0) {
-               exvector::iterator i = factors.begin(), end = factors.end();
-               while (i != end) {
+               exvector::iterator i = factors.begin();
+               while (i != factors.end()) {
                        *i = sqrfree(*i, newargs);
                        ++i;
                }
@@ -1796,10 +1799,16 @@ ex sqrfree(const ex &a, const lst &l)
        for (int p = 1; it!=itend; ++it, ++p)
                result *= power(*it, p);
 
-       // Yun's algorithm does not account for constant factors.  (For
-       // univariate polynomials it works only in the monic case.)  We can
-       // correct this by inserting what has been lost back into the result:
-       result = result * quo(tmp, result, x);
+       // Yun's algorithm does not account for constant factors.  (For univariate
+       // polynomials it works only in the monic case.)  We can correct this by
+       // inserting what has been lost back into the result.  For completeness
+       // we'll also have to recurse down that factor in the remaining variables.
+       if (newargs.nops()>0)
+               result *= sqrfree(quo(tmp, result, x), newargs);
+       else
+               result *= quo(tmp, result, x);
+
+       // Put in the reational overall factor again and return
        return result * lcm.inverse();
 }