X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=c1698ff79af62e63dbd30abd36c3499c7a0a8705;hp=514423c97930f4eff87ad74205553735069f3b26;hb=dad107ff48f68d45e72469a8716df375ae145cf3;hpb=55d35dcf72dc411c8265628fcad2bd67d320a8c9 diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 514423c9..c1698ff7 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -6,7 +6,7 @@ * computation, square-free factorization and rational function normalization. */ /* - * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2001 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 @@ -23,7 +23,6 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ -#include #include #include @@ -34,7 +33,6 @@ #include "constant.h" #include "expairseq.h" #include "fail.h" -#include "indexed.h" #include "inifcns.h" #include "lst.h" #include "mul.h" @@ -46,9 +44,7 @@ #include "symbol.h" #include "utils.h" -#ifndef NO_NAMESPACE_GINAC namespace GiNaC { -#endif // ndef NO_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 @@ -96,18 +92,18 @@ static struct _stat_print { * @return "false" if no symbol was found, "true" otherwise */ static bool get_first_symbol(const ex &e, const symbol *&x) { - if (is_ex_exactly_of_type(e, symbol)) { - x = static_cast(e.bp); - return true; - } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) { - for (unsigned i=0; i(e.bp); + return true; + } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) { + for (unsigned i=0; i sym_desc_vec; // Add symbol the sym_desc_vec (used internally by get_symbol_stats()) static void add_symbol(const symbol *s, sym_desc_vec &v) { - sym_desc_vec::iterator it = v.begin(), itend = v.end(); - while (it != itend) { - if (it->sym->compare(*s) == 0) // If it's already in there, don't add it a second time - return; - it++; - } - sym_desc d; - d.sym = s; - v.push_back(d); + sym_desc_vec::iterator it = v.begin(), itend = v.end(); + while (it != itend) { + if (it->sym->compare(*s) == 0) // If it's already in there, don't add it a second time + return; + it++; + } + sym_desc d; + d.sym = s; + v.push_back(d); } // Collect all symbols of an expression (used internally by get_symbol_stats()) static void collect_symbols(const ex &e, sym_desc_vec &v) { - if (is_ex_exactly_of_type(e, symbol)) { - add_symbol(static_cast(e.bp), v); - } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) { - for (unsigned i=0; i(e.bp), v); + } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) { + for (unsigned i=0; isym)); - int deg_b = b.degree(*(it->sym)); - it->deg_a = deg_a; - it->deg_b = deg_b; - it->max_deg = max(deg_a, deg_b); - it->ldeg_a = a.ldegree(*(it->sym)); - it->ldeg_b = b.ldegree(*(it->sym)); - it++; - } - sort(v.begin(), v.end()); + collect_symbols(a.eval(), v); // eval() to expand assigned symbols + collect_symbols(b.eval(), v); + sym_desc_vec::iterator it = v.begin(), itend = v.end(); + while (it != itend) { + int deg_a = a.degree(*(it->sym)); + int deg_b = b.degree(*(it->sym)); + it->deg_a = deg_a; + it->deg_b = deg_b; + it->max_deg = std::max(deg_a, deg_b); + it->max_lcnops = std::max(a.lcoeff(*(it->sym)).nops(), b.lcoeff(*(it->sym)).nops()); + it->ldeg_a = a.ldegree(*(it->sym)); + it->ldeg_b = b.ldegree(*(it->sym)); + it++; + } + sort(v.begin(), v.end()); #if 0 std::clog << "Symbols:\n"; it = v.begin(); itend = v.end(); while (it != itend) { - std::clog << " " << *it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << endl; + std::clog << " " << *it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << endl; std::clog << " lcoeff_a=" << a.lcoeff(*(it->sym)) << ", lcoeff_b=" << b.lcoeff(*(it->sym)) << endl; it++; } @@ -222,21 +228,21 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) // expression recursively (used internally by lcm_of_coefficients_denominators()) static numeric lcmcoeff(const ex &e, const numeric &l) { - if (e.info(info_flags::rational)) - return lcm(ex_to_numeric(e).denom(), l); - else if (is_ex_exactly_of_type(e, add)) { - numeric c = _num1(); - for (unsigned i=0; iinteger_content(); + GINAC_ASSERT(bp!=0); + return bp->integer_content(); } numeric basic::integer_content(void) const { - return _num1(); + return _num1(); } numeric numeric::integer_content(void) const { - return abs(*this); + return abs(*this); } numeric add::integer_content(void) const { - epvector::const_iterator it = seq.begin(); - epvector::const_iterator itend = seq.end(); - numeric c = _num0(); - while (it != itend) { - GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); - GINAC_ASSERT(is_ex_exactly_of_type(it->coeff,numeric)); - c = gcd(ex_to_numeric(it->coeff), c); - it++; - } - GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - c = gcd(ex_to_numeric(overall_coeff),c); - return c; + epvector::const_iterator it = seq.begin(); + epvector::const_iterator itend = seq.end(); + numeric c = _num0(); + while (it != itend) { + GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); + GINAC_ASSERT(is_ex_exactly_of_type(it->coeff,numeric)); + c = gcd(ex_to_numeric(it->coeff), c); + it++; + } + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); + c = gcd(ex_to_numeric(overall_coeff),c); + return c; } numeric mul::integer_content(void) const { #ifdef DO_GINAC_ASSERT - epvector::const_iterator it = seq.begin(); - epvector::const_iterator itend = seq.end(); - while (it != itend) { - GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); - ++it; - } + epvector::const_iterator it = seq.begin(); + epvector::const_iterator itend = seq.end(); + while (it != itend) { + GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); + ++it; + } #endif // def DO_GINAC_ASSERT - GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - return abs(ex_to_numeric(overall_coeff)); + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); + return abs(ex_to_numeric(overall_coeff)); } @@ -347,42 +353,42 @@ numeric mul::integer_content(void) const * @return quotient of a and b in Q[x] */ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args) { - if (b.is_zero()) - throw(std::overflow_error("quo: division by zero")); - if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) - return a / b; + if (b.is_zero()) + throw(std::overflow_error("quo: division by zero")); + if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) + return a / b; #if FAST_COMPARE - if (a.is_equal(b)) - return _ex1(); + if (a.is_equal(b)) + return _ex1(); #endif - if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) - throw(std::invalid_argument("quo: arguments must be polynomials over the rationals")); - - // Polynomial long division - ex q = _ex0(); - ex r = a.expand(); - if (r.is_zero()) - return r; - 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); - while (rdeg >= bdeg) { - ex term, rcoeff = r.coeff(x, rdeg); - if (blcoeff_is_numeric) - term = rcoeff / blcoeff; - else { - if (!divide(rcoeff, blcoeff, term, false)) - return *new ex(fail()); - } - term *= power(x, rdeg - bdeg); - q += term; - r -= (term * b).expand(); - if (r.is_zero()) - break; - rdeg = r.degree(x); - } - return q; + if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) + throw(std::invalid_argument("quo: arguments must be polynomials over the rationals")); + + // Polynomial long division + ex q = _ex0(); + ex r = a.expand(); + if (r.is_zero()) + return r; + 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); + while (rdeg >= bdeg) { + ex term, rcoeff = r.coeff(x, rdeg); + if (blcoeff_is_numeric) + term = rcoeff / blcoeff; + else { + if (!divide(rcoeff, blcoeff, term, false)) + return *new ex(fail()); + } + term *= power(x, rdeg - bdeg); + q += term; + r -= (term * b).expand(); + if (r.is_zero()) + break; + rdeg = r.degree(x); + } + return q; } @@ -397,44 +403,44 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args) * @return remainder of a(x) and b(x) in Q[x] */ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) { - if (b.is_zero()) - throw(std::overflow_error("rem: 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 (b.is_zero()) + throw(std::overflow_error("rem: 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 FAST_COMPARE - if (a.is_equal(b)) - return _ex0(); + if (a.is_equal(b)) + return _ex0(); #endif - if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) - throw(std::invalid_argument("rem: arguments must be polynomials over the rationals")); - - // Polynomial long division - ex r = a.expand(); - if (r.is_zero()) - return r; - 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); - while (rdeg >= bdeg) { - ex term, rcoeff = r.coeff(x, rdeg); - if (blcoeff_is_numeric) - term = rcoeff / blcoeff; - else { - if (!divide(rcoeff, blcoeff, term, false)) - return *new ex(fail()); - } - term *= power(x, rdeg - bdeg); - r -= (term * b).expand(); - if (r.is_zero()) - break; - rdeg = r.degree(x); - } - return r; + if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) + throw(std::invalid_argument("rem: arguments must be polynomials over the rationals")); + + // Polynomial long division + ex r = a.expand(); + if (r.is_zero()) + return r; + 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); + while (rdeg >= bdeg) { + ex term, rcoeff = r.coeff(x, rdeg); + if (blcoeff_is_numeric) + term = rcoeff / blcoeff; + else { + if (!divide(rcoeff, blcoeff, term, false)) + return *new ex(fail()); + } + term *= power(x, rdeg - bdeg); + r -= (term * b).expand(); + if (r.is_zero()) + break; + rdeg = r.degree(x); + } + return r; } @@ -448,45 +454,45 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) * @return pseudo-remainder of a(x) and b(x) in Z[x] */ ex prem(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(); - - int delta = rdeg - bdeg + 1, i = 0; - 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); - i++; - } - return power(blcoeff, delta - i) * r; + 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(); + + int delta = rdeg - bdeg + 1, i = 0; + 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); + i++; + } + return power(blcoeff, delta - i) * r; } @@ -501,43 +507,43 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) 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; + 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; } @@ -552,54 +558,54 @@ ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args) * "false" otherwise */ 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()) - return true; - if (is_ex_exactly_of_type(b, numeric)) { - q = a / b; - return true; - } else if (is_ex_exactly_of_type(a, numeric)) - return false; + q = _ex0(); + if (b.is_zero()) + throw(std::overflow_error("divide: division by zero")); + if (a.is_zero()) + return true; + if (is_ex_exactly_of_type(b, numeric)) { + q = a / b; + return true; + } else if (is_ex_exactly_of_type(a, numeric)) + return false; #if FAST_COMPARE - if (a.is_equal(b)) { - q = _ex1(); - return true; - } + if (a.is_equal(b)) { + q = _ex1(); + return true; + } #endif - 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 - const symbol *x; - if (!get_first_symbol(a, x) && !get_first_symbol(b, x)) - throw(std::invalid_argument("invalid expression in divide()")); - - // Polynomial long division (recursive) - ex r = a.expand(); - if (r.is_zero()) - 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); - while (rdeg >= bdeg) { - ex term, rcoeff = r.coeff(*x, rdeg); - if (blcoeff_is_numeric) - term = rcoeff / blcoeff; - else - if (!divide(rcoeff, blcoeff, term, false)) - return false; - term *= power(*x, rdeg - bdeg); - q += term; - r -= (term * b).expand(); - if (r.is_zero()) - return true; - rdeg = r.degree(*x); - } - return false; + 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 + const symbol *x; + if (!get_first_symbol(a, x) && !get_first_symbol(b, x)) + throw(std::invalid_argument("invalid expression in divide()")); + + // Polynomial long division (recursive) + ex r = a.expand(); + if (r.is_zero()) + 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); + while (rdeg >= bdeg) { + ex term, rcoeff = r.coeff(*x, rdeg); + if (blcoeff_is_numeric) + term = rcoeff / blcoeff; + else + if (!divide(rcoeff, blcoeff, term, false)) + return false; + term *= power(*x, rdeg - bdeg); + q += term; + r -= (term * b).expand(); + if (r.is_zero()) + return true; + rdeg = r.degree(*x); + } + return false; } @@ -612,10 +618,11 @@ typedef std::pair ex2; typedef std::pair exbool; struct ex2_less { - bool operator() (const ex2 p, const ex2 q) const - { - return p.first.compare(q.first) < 0 || (!(q.first.compare(p.first) < 0) && p.second.compare(q.second) < 0); - } + bool operator() (const ex2 &p, const ex2 &q) const + { + int cmp = p.first.compare(q.first); + return ((cmp<0) || (!(cmp>0) && p.second.compare(q.second)<0)); + } }; typedef std::map ex2_exbool_remember; @@ -640,127 +647,127 @@ typedef std::map ex2_exbool_remember; * @see get_symbol_stats, heur_gcd */ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_iterator var) { - q = _ex0(); - if (b.is_zero()) - throw(std::overflow_error("divide_in_z: division by zero")); - if (b.is_equal(_ex1())) { - q = a; - return true; - } - if (is_ex_exactly_of_type(a, numeric)) { - if (is_ex_exactly_of_type(b, numeric)) { - q = a / b; - return q.info(info_flags::integer); - } else - return false; - } + q = _ex0(); + if (b.is_zero()) + throw(std::overflow_error("divide_in_z: division by zero")); + if (b.is_equal(_ex1())) { + q = a; + return true; + } + if (is_ex_exactly_of_type(a, numeric)) { + if (is_ex_exactly_of_type(b, numeric)) { + q = a / b; + return q.info(info_flags::integer); + } else + return false; + } #if FAST_COMPARE - if (a.is_equal(b)) { - q = _ex1(); - return true; - } + if (a.is_equal(b)) { + q = _ex1(); + return true; + } #endif #if USE_REMEMBER - // Remembering - static ex2_exbool_remember dr_remember; - ex2_exbool_remember::const_iterator remembered = dr_remember.find(ex2(a, b)); - if (remembered != dr_remember.end()) { - q = remembered->second.first; - return remembered->second.second; - } + // Remembering + static ex2_exbool_remember dr_remember; + ex2_exbool_remember::const_iterator remembered = dr_remember.find(ex2(a, b)); + if (remembered != dr_remember.end()) { + q = remembered->second.first; + return remembered->second.second; + } #endif - // Main symbol - const symbol *x = var->sym; + // Main symbol + const symbol *x = var->sym; - // Compare degrees - int adeg = a.degree(*x), bdeg = b.degree(*x); - if (bdeg > adeg) - return false; + // Compare degrees + int adeg = a.degree(*x), bdeg = b.degree(*x); + if (bdeg > adeg) + return false; #if USE_TRIAL_DIVISION - // Trial division with polynomial interpolation - int i, k; - - // Compute values at evaluation points 0..adeg - vector alpha; alpha.reserve(adeg + 1); - exvector u; u.reserve(adeg + 1); - numeric point = _num0(); - ex c; - for (i=0; i<=adeg; i++) { - ex bs = b.subs(*x == point); - while (bs.is_zero()) { - point += _num1(); - bs = b.subs(*x == point); - } - if (!divide_in_z(a.subs(*x == point), bs, c, var+1)) - return false; - alpha.push_back(point); - u.push_back(c); - point += _num1(); - } - - // Compute inverses - vector rcp; rcp.reserve(adeg + 1); - rcp.push_back(_num0()); - for (k=1; k<=adeg; k++) { - numeric product = alpha[k] - alpha[0]; - for (i=1; i=0; i--) - temp = temp * (alpha[k] - alpha[i]) + v[i]; - v.push_back((u[k] - temp) * rcp[k]); - } - - // Convert from Newton form to standard form - c = v[adeg]; - for (k=adeg-1; k>=0; k--) - c = c * (*x - alpha[k]) + v[k]; - - if (c.degree(*x) == (adeg - bdeg)) { - q = c.expand(); - return true; - } else - return false; + // Trial division with polynomial interpolation + int i, k; + + // Compute values at evaluation points 0..adeg + vector alpha; alpha.reserve(adeg + 1); + exvector u; u.reserve(adeg + 1); + numeric point = _num0(); + ex c; + for (i=0; i<=adeg; i++) { + ex bs = b.subs(*x == point); + while (bs.is_zero()) { + point += _num1(); + bs = b.subs(*x == point); + } + if (!divide_in_z(a.subs(*x == point), bs, c, var+1)) + return false; + alpha.push_back(point); + u.push_back(c); + point += _num1(); + } + + // Compute inverses + vector rcp; rcp.reserve(adeg + 1); + rcp.push_back(_num0()); + for (k=1; k<=adeg; k++) { + numeric product = alpha[k] - alpha[0]; + for (i=1; i=0; i--) + temp = temp * (alpha[k] - alpha[i]) + v[i]; + v.push_back((u[k] - temp) * rcp[k]); + } + + // Convert from Newton form to standard form + c = v[adeg]; + for (k=adeg-1; k>=0; k--) + c = c * (*x - alpha[k]) + v[k]; + + if (c.degree(*x) == (adeg - bdeg)) { + q = c.expand(); + return true; + } else + return false; #else - // Polynomial long division (recursive) - ex r = a.expand(); - if (r.is_zero()) - return true; - int rdeg = adeg; - ex eb = b.expand(); - ex blcoeff = eb.coeff(*x, bdeg); - while (rdeg >= bdeg) { - ex term, rcoeff = r.coeff(*x, rdeg); - if (!divide_in_z(rcoeff, blcoeff, term, var+1)) - break; - term = (term * power(*x, rdeg - bdeg)).expand(); - q += term; - r -= (term * eb).expand(); - if (r.is_zero()) { + // Polynomial long division (recursive) + ex r = a.expand(); + if (r.is_zero()) + return true; + int rdeg = adeg; + ex eb = b.expand(); + ex blcoeff = eb.coeff(*x, bdeg); + while (rdeg >= bdeg) { + ex term, rcoeff = r.coeff(*x, rdeg); + if (!divide_in_z(rcoeff, blcoeff, term, var+1)) + break; + term = (term * power(*x, rdeg - bdeg)).expand(); + q += term; + r -= (term * eb).expand(); + if (r.is_zero()) { #if USE_REMEMBER - dr_remember[ex2(a, b)] = exbool(q, true); + dr_remember[ex2(a, b)] = exbool(q, true); #endif - return true; - } - rdeg = r.degree(*x); - } + return true; + } + rdeg = r.degree(*x); + } #if USE_REMEMBER - dr_remember[ex2(a, b)] = exbool(q, false); + dr_remember[ex2(a, b)] = exbool(q, false); #endif - return false; + return false; #endif } @@ -779,16 +786,16 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite * @see ex::content, ex::primpart */ ex ex::unit(const symbol &x) const { - ex c = expand().lcoeff(x); - if (is_ex_exactly_of_type(c, numeric)) - return c < _ex0() ? _ex_1() : _ex1(); - else { - const symbol *y; - if (get_first_symbol(c, y)) - return c.unit(*y); - else - throw(std::invalid_argument("invalid expression in unit()")); - } + ex c = expand().lcoeff(x); + if (is_ex_exactly_of_type(c, numeric)) + return c < _ex0() ? _ex_1() : _ex1(); + else { + const symbol *y; + if (get_first_symbol(c, y)) + return c.unit(*y); + else + throw(std::invalid_argument("invalid expression in unit()")); + } } @@ -801,30 +808,30 @@ ex ex::unit(const symbol &x) const * @see ex::unit, ex::primpart */ ex ex::content(const symbol &x) const { - if (is_zero()) - return _ex0(); - if (is_ex_exactly_of_type(*this, numeric)) - return info(info_flags::negative) ? -*this : *this; - ex e = expand(); - if (e.is_zero()) - return _ex0(); - - // First, try the integer content - ex c = e.integer_content(); - ex r = e / c; - ex lcoeff = r.lcoeff(x); - if (lcoeff.info(info_flags::integer)) - return c; - - // GCD of all coefficients - int deg = e.degree(x); - int ldeg = e.ldegree(x); - if (deg == ldeg) - return e.lcoeff(x) / e.unit(x); - c = _ex0(); - for (int i=ldeg; i<=deg; i++) - c = gcd(e.coeff(x, i), c, NULL, NULL, false); - return c; + if (is_zero()) + return _ex0(); + if (is_ex_exactly_of_type(*this, numeric)) + return info(info_flags::negative) ? -*this : *this; + ex e = expand(); + if (e.is_zero()) + return _ex0(); + + // First, try the integer content + ex c = e.integer_content(); + ex r = e / c; + ex lcoeff = r.lcoeff(x); + if (lcoeff.info(info_flags::integer)) + return c; + + // GCD of all coefficients + int deg = e.degree(x); + int ldeg = e.ldegree(x); + if (deg == ldeg) + return e.lcoeff(x) / e.unit(x); + c = _ex0(); + for (int i=ldeg; i<=deg; i++) + c = gcd(e.coeff(x, i), c, NULL, NULL, false); + return c; } @@ -837,19 +844,19 @@ ex ex::content(const symbol &x) const * @see ex::unit, ex::content */ ex ex::primpart(const symbol &x) const { - if (is_zero()) - return _ex0(); - if (is_ex_exactly_of_type(*this, numeric)) - return _ex1(); - - ex c = content(x); - if (c.is_zero()) - return _ex0(); - ex u = unit(x); - if (is_ex_exactly_of_type(c, numeric)) - return *this / (c * u); - else - return quo(*this, c * u, x, false); + if (is_zero()) + return _ex0(); + if (is_ex_exactly_of_type(*this, numeric)) + return _ex1(); + + ex c = content(x); + if (c.is_zero()) + return _ex0(); + ex u = unit(x); + if (is_ex_exactly_of_type(c, numeric)) + return *this / (c * u); + else + return quo(*this, c * u, x, false); } @@ -862,18 +869,18 @@ ex ex::primpart(const symbol &x) const * @return primitive part */ ex ex::primpart(const symbol &x, const ex &c) const { - if (is_zero()) - return _ex0(); - if (c.is_zero()) - return _ex0(); - if (is_ex_exactly_of_type(*this, numeric)) - return _ex1(); - - ex u = unit(x); - if (is_ex_exactly_of_type(c, numeric)) - return *this / (c * u); - else - return quo(*this, c * u, x, false); + if (is_zero()) + return _ex0(); + if (c.is_zero()) + return _ex0(); + if (is_ex_exactly_of_type(*this, numeric)) + return _ex1(); + + ex u = unit(x); + if (is_ex_exactly_of_type(c, numeric)) + return *this / (c * u); + else + return quo(*this, c * u, x, false); } @@ -895,31 +902,31 @@ static ex eu_gcd(const ex &a, const ex &b, const symbol *x) { //std::clog << "eu_gcd(" << a << "," << b << ")\n"; - // Sort c and d so that c has higher degree - ex c, d; - int adeg = a.degree(*x), bdeg = b.degree(*x); - if (adeg >= bdeg) { - c = a; - d = b; - } else { - c = b; - d = a; - } + // Sort c and d so that c has higher degree + ex c, d; + int adeg = a.degree(*x), bdeg = b.degree(*x); + if (adeg >= bdeg) { + c = a; + d = b; + } else { + c = b; + d = a; + } // Normalize in Q[x] c = c / c.lcoeff(*x); d = d / d.lcoeff(*x); // Euclidean algorithm - ex r; - for (;;) { + ex r; + for (;;) { //std::clog << " d = " << d << endl; - r = rem(c, d, *x, false); - if (r.is_zero()) - return d / d.lcoeff(*x); - c = d; + r = rem(c, d, *x, false); + if (r.is_zero()) + return d / d.lcoeff(*x); + c = d; d = r; - } + } } @@ -937,30 +944,30 @@ static ex euprem_gcd(const ex &a, const ex &b, const symbol *x) { //std::clog << "euprem_gcd(" << a << "," << b << ")\n"; - // Sort c and d so that c has higher degree - ex c, d; - int adeg = a.degree(*x), bdeg = b.degree(*x); - if (adeg >= bdeg) { - c = a; - d = b; - } else { - c = b; - d = a; - } + // Sort c and d so that c has higher degree + ex c, d; + int adeg = a.degree(*x), bdeg = b.degree(*x); + if (adeg >= bdeg) { + c = a; + d = b; + } else { + c = b; + d = a; + } // Calculate GCD of contents ex gamma = gcd(c.content(*x), d.content(*x), NULL, NULL, false); // Euclidean algorithm with pseudo-remainders - ex r; - for (;;) { + ex r; + for (;;) { //std::clog << " d = " << d << endl; - r = prem(c, d, *x, false); - if (r.is_zero()) - return d.primpart(*x) * gamma; - c = d; + r = prem(c, d, *x, false); + if (r.is_zero()) + return d.primpart(*x) * gamma; + c = d; d = r; - } + } } @@ -978,39 +985,39 @@ static ex peu_gcd(const ex &a, const ex &b, const symbol *x) { //std::clog << "peu_gcd(" << a << "," << b << ")\n"; - // Sort c and d so that c has higher degree - ex c, d; - int adeg = a.degree(*x), bdeg = b.degree(*x); - int ddeg; - if (adeg >= bdeg) { - c = a; - d = b; - ddeg = bdeg; - } else { - c = b; - d = a; - ddeg = adeg; - } - - // Remove content from c and d, to be attached to GCD later - ex cont_c = c.content(*x); - ex cont_d = d.content(*x); - ex gamma = gcd(cont_c, cont_d, NULL, NULL, false); - if (ddeg == 0) - return gamma; - c = c.primpart(*x, cont_c); - d = d.primpart(*x, cont_d); - - // Euclidean algorithm with content removal + // Sort c and d so that c has higher degree + ex c, d; + int adeg = a.degree(*x), bdeg = b.degree(*x); + int ddeg; + if (adeg >= bdeg) { + c = a; + d = b; + ddeg = bdeg; + } else { + c = b; + d = a; + ddeg = adeg; + } + + // Remove content from c and d, to be attached to GCD later + ex cont_c = c.content(*x); + ex cont_d = d.content(*x); + ex gamma = gcd(cont_c, cont_d, NULL, NULL, false); + if (ddeg == 0) + return gamma; + c = c.primpart(*x, cont_c); + d = d.primpart(*x, cont_d); + + // Euclidean algorithm with content removal ex r; - for (;;) { + for (;;) { //std::clog << " d = " << d << endl; - r = prem(c, d, *x, false); - if (r.is_zero()) - return gamma * d; - c = d; + r = prem(c, d, *x, false); + if (r.is_zero()) + return gamma * d; + c = d; d = r.primpart(*x); - } + } } @@ -1027,57 +1034,57 @@ static ex red_gcd(const ex &a, const ex &b, const symbol *x) { //std::clog << "red_gcd(" << a << "," << b << ")\n"; - // Sort c and d so that c has higher degree - ex c, d; - int adeg = a.degree(*x), bdeg = b.degree(*x); - int cdeg, ddeg; - if (adeg >= bdeg) { - c = a; - d = b; - cdeg = adeg; - ddeg = bdeg; - } else { - c = b; - d = a; - cdeg = bdeg; - ddeg = adeg; - } - - // Remove content from c and d, to be attached to GCD later - ex cont_c = c.content(*x); - ex cont_d = d.content(*x); - ex gamma = gcd(cont_c, cont_d, NULL, NULL, false); - if (ddeg == 0) - return gamma; - c = c.primpart(*x, cont_c); - d = d.primpart(*x, cont_d); - - // First element of divisor sequence - ex r, ri = _ex1(); - int delta = cdeg - ddeg; - - for (;;) { - // Calculate polynomial pseudo-remainder + // Sort c and d so that c has higher degree + ex c, d; + int adeg = a.degree(*x), bdeg = b.degree(*x); + int cdeg, ddeg; + if (adeg >= bdeg) { + c = a; + d = b; + cdeg = adeg; + ddeg = bdeg; + } else { + c = b; + d = a; + cdeg = bdeg; + ddeg = adeg; + } + + // Remove content from c and d, to be attached to GCD later + ex cont_c = c.content(*x); + ex cont_d = d.content(*x); + ex gamma = gcd(cont_c, cont_d, NULL, NULL, false); + if (ddeg == 0) + return gamma; + c = c.primpart(*x, cont_c); + d = d.primpart(*x, cont_d); + + // First element of divisor sequence + ex r, ri = _ex1(); + int delta = cdeg - ddeg; + + for (;;) { + // Calculate polynomial pseudo-remainder //std::clog << " d = " << d << endl; - r = prem(c, d, *x, false); - if (r.is_zero()) - return gamma * d.primpart(*x); - c = d; - cdeg = ddeg; - - if (!divide(r, pow(ri, delta), d, false)) - throw(std::runtime_error("invalid expression in red_gcd(), division failed")); - ddeg = d.degree(*x); - if (ddeg == 0) { - if (is_ex_exactly_of_type(r, numeric)) - return gamma; - else - return gamma * r.primpart(*x); - } - - ri = c.expand().lcoeff(*x); - delta = cdeg - ddeg; - } + r = prem(c, d, *x, false); + if (r.is_zero()) + return gamma * d.primpart(*x); + c = d; + cdeg = ddeg; + + if (!divide(r, pow(ri, delta), d, false)) + throw(std::runtime_error("invalid expression in red_gcd(), division failed")); + ddeg = d.degree(*x); + if (ddeg == 0) { + if (is_ex_exactly_of_type(r, numeric)) + return gamma; + else + return gamma * r.primpart(*x); + } + + ri = c.expand().lcoeff(*x); + delta = cdeg - ddeg; + } } @@ -1097,68 +1104,68 @@ static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var) sr_gcd_called++; #endif - // The first symbol is our main variable - const symbol &x = *(var->sym); - - // Sort c and d so that c has higher degree - ex c, d; - int adeg = a.degree(x), bdeg = b.degree(x); - int cdeg, ddeg; - if (adeg >= bdeg) { - c = a; - d = b; - cdeg = adeg; - ddeg = bdeg; - } else { - c = b; - d = a; - cdeg = bdeg; - ddeg = adeg; - } - - // Remove content from c and d, to be attached to GCD later - ex cont_c = c.content(x); - ex cont_d = d.content(x); - ex gamma = gcd(cont_c, cont_d, NULL, NULL, false); - if (ddeg == 0) - return gamma; - c = c.primpart(x, cont_c); - d = d.primpart(x, cont_d); + // The first symbol is our main variable + const symbol &x = *(var->sym); + + // Sort c and d so that c has higher degree + ex c, d; + int adeg = a.degree(x), bdeg = b.degree(x); + int cdeg, ddeg; + if (adeg >= bdeg) { + c = a; + d = b; + cdeg = adeg; + ddeg = bdeg; + } else { + c = b; + d = a; + cdeg = bdeg; + ddeg = adeg; + } + + // Remove content from c and d, to be attached to GCD later + ex cont_c = c.content(x); + ex cont_d = d.content(x); + ex gamma = gcd(cont_c, cont_d, NULL, NULL, false); + if (ddeg == 0) + return gamma; + c = c.primpart(x, cont_c); + d = d.primpart(x, cont_d); //std::clog << " content " << gamma << " removed, continuing with sr_gcd(" << c << "," << d << ")\n"; - // First element of subresultant sequence - ex r = _ex0(), ri = _ex1(), psi = _ex1(); - int delta = cdeg - ddeg; + // First element of subresultant sequence + ex r = _ex0(), ri = _ex1(), psi = _ex1(); + int delta = cdeg - ddeg; - for (;;) { - // Calculate polynomial pseudo-remainder + for (;;) { + // Calculate polynomial pseudo-remainder //std::clog << " start of loop, psi = " << psi << ", calculating pseudo-remainder...\n"; //std::clog << " d = " << d << endl; - r = prem(c, d, x, false); - if (r.is_zero()) - return gamma * d.primpart(x); - c = d; - cdeg = ddeg; + r = prem(c, d, x, false); + if (r.is_zero()) + return gamma * d.primpart(x); + c = d; + cdeg = ddeg; //std::clog << " dividing...\n"; - if (!divide_in_z(r, ri * pow(psi, delta), d, var)) - throw(std::runtime_error("invalid expression in sr_gcd(), division failed")); - ddeg = d.degree(x); - if (ddeg == 0) { - if (is_ex_exactly_of_type(r, numeric)) - return gamma; - else - return gamma * r.primpart(x); - } - - // Next element of subresultant sequence + if (!divide_in_z(r, ri * pow(psi, delta), d, var)) + throw(std::runtime_error("invalid expression in sr_gcd(), division failed")); + ddeg = d.degree(x); + if (ddeg == 0) { + if (is_ex_exactly_of_type(r, numeric)) + return gamma; + else + return gamma * r.primpart(x); + } + + // Next element of subresultant sequence //std::clog << " calculating next subresultant...\n"; - ri = c.expand().lcoeff(x); - if (delta == 1) - psi = ri; - else if (delta) - divide_in_z(pow(ri, delta), pow(psi, delta-1), psi, var+1); - delta = cdeg - ddeg; - } + ri = c.expand().lcoeff(x); + if (delta == 1) + psi = ri; + else if (delta) + divide_in_z(pow(ri, delta), pow(psi, delta-1), psi, var+1); + delta = cdeg - ddeg; + } } @@ -1170,49 +1177,51 @@ static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var) * @see heur_gcd */ numeric ex::max_coefficient(void) const { - GINAC_ASSERT(bp!=0); - return bp->max_coefficient(); + GINAC_ASSERT(bp!=0); + return bp->max_coefficient(); } +/** Implementation ex::max_coefficient(). + * @see heur_gcd */ numeric basic::max_coefficient(void) const { - return _num1(); + return _num1(); } numeric numeric::max_coefficient(void) const { - return abs(*this); + return abs(*this); } numeric add::max_coefficient(void) const { - epvector::const_iterator it = seq.begin(); - epvector::const_iterator itend = seq.end(); - GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - numeric cur_max = abs(ex_to_numeric(overall_coeff)); - while (it != itend) { - numeric a; - GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); - a = abs(ex_to_numeric(it->coeff)); - if (a > cur_max) - cur_max = a; - it++; - } - return cur_max; + epvector::const_iterator it = seq.begin(); + epvector::const_iterator itend = seq.end(); + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); + numeric cur_max = abs(ex_to_numeric(overall_coeff)); + while (it != itend) { + numeric a; + GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); + a = abs(ex_to_numeric(it->coeff)); + if (a > cur_max) + cur_max = a; + it++; + } + return cur_max; } numeric mul::max_coefficient(void) const { #ifdef DO_GINAC_ASSERT - epvector::const_iterator it = seq.begin(); - epvector::const_iterator itend = seq.end(); - while (it != itend) { - GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); - it++; - } + epvector::const_iterator it = seq.begin(); + epvector::const_iterator itend = seq.end(); + while (it != itend) { + GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); + it++; + } #endif // def DO_GINAC_ASSERT - GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - return abs(ex_to_numeric(overall_coeff)); + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); + return abs(ex_to_numeric(overall_coeff)); } @@ -1225,70 +1234,54 @@ numeric mul::max_coefficient(void) const * @see heur_gcd */ ex ex::smod(const numeric &xi) const { - GINAC_ASSERT(bp!=0); - return bp->smod(xi); + GINAC_ASSERT(bp!=0); + return bp->smod(xi); } ex basic::smod(const numeric &xi) const { - return *this; + return *this; } ex numeric::smod(const numeric &xi) const { -#ifndef NO_NAMESPACE_GINAC - return GiNaC::smod(*this, xi); -#else // ndef NO_NAMESPACE_GINAC - return ::smod(*this, xi); -#endif // ndef NO_NAMESPACE_GINAC + return GiNaC::smod(*this, xi); } ex add::smod(const numeric &xi) const { - epvector newseq; - newseq.reserve(seq.size()+1); - epvector::const_iterator it = seq.begin(); - epvector::const_iterator itend = seq.end(); - while (it != itend) { - GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); -#ifndef NO_NAMESPACE_GINAC - numeric coeff = GiNaC::smod(ex_to_numeric(it->coeff), xi); -#else // ndef NO_NAMESPACE_GINAC - numeric coeff = ::smod(ex_to_numeric(it->coeff), xi); -#endif // ndef NO_NAMESPACE_GINAC - if (!coeff.is_zero()) - newseq.push_back(expair(it->rest, coeff)); - it++; - } - GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); -#ifndef NO_NAMESPACE_GINAC - numeric coeff = GiNaC::smod(ex_to_numeric(overall_coeff), xi); -#else // ndef NO_NAMESPACE_GINAC - numeric coeff = ::smod(ex_to_numeric(overall_coeff), xi); -#endif // ndef NO_NAMESPACE_GINAC - return (new add(newseq,coeff))->setflag(status_flags::dynallocated); + epvector newseq; + newseq.reserve(seq.size()+1); + epvector::const_iterator it = seq.begin(); + epvector::const_iterator itend = seq.end(); + while (it != itend) { + GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); + numeric coeff = GiNaC::smod(ex_to_numeric(it->coeff), xi); + if (!coeff.is_zero()) + newseq.push_back(expair(it->rest, coeff)); + it++; + } + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); + numeric coeff = GiNaC::smod(ex_to_numeric(overall_coeff), xi); + return (new add(newseq,coeff))->setflag(status_flags::dynallocated); } ex mul::smod(const numeric &xi) const { #ifdef DO_GINAC_ASSERT - epvector::const_iterator it = seq.begin(); - epvector::const_iterator itend = seq.end(); - while (it != itend) { - GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); - it++; - } + epvector::const_iterator it = seq.begin(); + epvector::const_iterator itend = seq.end(); + while (it != itend) { + GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); + it++; + } #endif // def DO_GINAC_ASSERT - mul * mulcopyp=new mul(*this); - GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); -#ifndef NO_NAMESPACE_GINAC - mulcopyp->overall_coeff = GiNaC::smod(ex_to_numeric(overall_coeff),xi); -#else // ndef NO_NAMESPACE_GINAC - mulcopyp->overall_coeff = ::smod(ex_to_numeric(overall_coeff),xi); -#endif // ndef NO_NAMESPACE_GINAC - mulcopyp->clearflag(status_flags::evaluated); - mulcopyp->clearflag(status_flags::hash_calculated); - return mulcopyp->setflag(status_flags::dynallocated); + mul * mulcopyp=new mul(*this); + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); + 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); } @@ -1336,61 +1329,62 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const return *new ex(fail()); // GCD of two numeric values -> CLN - 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 (ca) - *ca = ex_to_numeric(a) / g; - if (cb) - *cb = ex_to_numeric(b) / g; - return g; - } - - // The first symbol is our main variable - const symbol &x = *(var->sym); - - // Remove integer content - numeric gc = gcd(a.integer_content(), b.integer_content()); - numeric rgc = gc.inverse(); - ex p = a * rgc; - ex q = b * rgc; - int maxdeg = max(p.degree(x), q.degree(x)); - - // Find evaluation point - numeric mp = p.max_coefficient(), mq = q.max_coefficient(); - numeric xi; - if (mp > mq) - xi = mq * _num2() + _num2(); - else - xi = mp * _num2() + _num2(); - - // 6 tries maximum - for (int t=0; t<6; t++) { - if (xi.int_length() * maxdeg > 100000) { + 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 (ca) + *ca = ex_to_numeric(a) / g; + if (cb) + *cb = ex_to_numeric(b) / g; + return g; + } + + // The first symbol is our main variable + const symbol &x = *(var->sym); + + // Remove integer content + numeric gc = gcd(a.integer_content(), b.integer_content()); + numeric rgc = gc.inverse(); + ex p = a * rgc; + ex q = b * rgc; + int maxdeg = std::max(p.degree(x),q.degree(x)); + + // Find evaluation point + numeric mp = p.max_coefficient(); + numeric mq = q.max_coefficient(); + numeric xi; + if (mp > mq) + xi = mq * _num2() + _num2(); + else + xi = mp * _num2() + _num2(); + + // 6 tries maximum + for (int t=0; t<6; t++) { + if (xi.int_length() * maxdeg > 100000) { //std::clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << endl; - throw gcdheu_failed(); + throw gcdheu_failed(); } - // Apply evaluation homomorphism and calculate GCD + // Apply evaluation homomorphism and calculate GCD ex cp, cq; - ex gamma = heur_gcd(p.subs(x == xi), q.subs(x == xi), &cp, &cq, var+1).expand(); - if (!is_ex_exactly_of_type(gamma, fail)) { + ex gamma = heur_gcd(p.subs(x == xi), q.subs(x == xi), &cp, &cq, var+1).expand(); + if (!is_ex_exactly_of_type(gamma, fail)) { - // Reconstruct polynomial from GCD of mapped polynomials + // Reconstruct polynomial from GCD of mapped polynomials ex g = interpolate(gamma, xi, x); - // Remove integer content - g /= g.integer_content(); - - // If the calculated polynomial divides both p and q, this is the GCD - ex dummy; - if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) { - g *= gc; - ex lc = g.lcoeff(x); - if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) - return -g; - else - return g; - } + // Remove integer content + g /= g.integer_content(); + + // If the calculated polynomial divides both p and q, this is the GCD + ex dummy; + if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) { + g *= gc; + ex lc = g.lcoeff(x); + if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) + return -g; + else + return g; + } #if 0 cp = interpolate(cp, xi, x); if (divide_in_z(cp, p, g, var)) { @@ -1398,11 +1392,11 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const g *= gc; if (ca) *ca = cp; - ex lc = g.lcoeff(x); - if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) - return -g; - else - return g; + ex lc = g.lcoeff(x); + if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) + return -g; + else + return g; } } cq = interpolate(cq, xi, x); @@ -1411,20 +1405,20 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const g *= gc; if (cb) *cb = cq; - ex lc = g.lcoeff(x); - if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) - return -g; - else - return g; + ex lc = g.lcoeff(x); + if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) + return -g; + else + return g; } } #endif - } + } - // Next evaluation point - xi = iquo(xi * isqrt(isqrt(xi)) * numeric(73794), numeric(27011)); - } - return *new ex(fail()); + // Next evaluation point + xi = iquo(xi * isqrt(isqrt(xi)) * numeric(73794), numeric(27011)); + } + return *new ex(fail()); } @@ -1444,8 +1438,8 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) #endif // GCD of numerics -> CLN - 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(a, numeric) && is_ex_exactly_of_type(b, numeric)) { + numeric g = gcd(ex_to_numeric(a), ex_to_numeric(b)); if (ca || cb) { if (g.is_zero()) { if (ca) @@ -1453,19 +1447,19 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) if (cb) *cb = _ex0(); } else { - if (ca) - *ca = ex_to_numeric(a) / g; - if (cb) - *cb = ex_to_numeric(b) / g; + if (ca) + *ca = ex_to_numeric(a) / g; + if (cb) + *cb = ex_to_numeric(b) / g; } } - return g; - } + return g; + } // Check arguments - if (check_args && !a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) { - throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals")); - } + if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) { + throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals")); + } // Partially factored cases (to avoid expanding large expressions) if (is_ex_exactly_of_type(a, mul)) { @@ -1551,83 +1545,83 @@ factored_b: } #endif - // Some trivial cases + // Some trivial cases ex aex = a.expand(), bex = b.expand(); - if (aex.is_zero()) { - if (ca) - *ca = _ex0(); - if (cb) - *cb = _ex1(); - return b; - } - if (bex.is_zero()) { - if (ca) - *ca = _ex1(); - if (cb) - *cb = _ex0(); - return a; - } - if (aex.is_equal(_ex1()) || bex.is_equal(_ex1())) { - if (ca) - *ca = a; - if (cb) - *cb = b; - return _ex1(); - } + if (aex.is_zero()) { + if (ca) + *ca = _ex0(); + if (cb) + *cb = _ex1(); + return b; + } + if (bex.is_zero()) { + if (ca) + *ca = _ex1(); + if (cb) + *cb = _ex0(); + return a; + } + if (aex.is_equal(_ex1()) || bex.is_equal(_ex1())) { + if (ca) + *ca = a; + if (cb) + *cb = b; + return _ex1(); + } #if FAST_COMPARE - if (a.is_equal(b)) { - if (ca) - *ca = _ex1(); - if (cb) - *cb = _ex1(); - return a; - } + if (a.is_equal(b)) { + if (ca) + *ca = _ex1(); + if (cb) + *cb = _ex1(); + return a; + } #endif - // Gather symbol statistics - sym_desc_vec sym_stats; - get_symbol_stats(a, b, sym_stats); + // Gather symbol statistics + sym_desc_vec sym_stats; + get_symbol_stats(a, b, sym_stats); - // The symbol with least degree is our main variable - sym_desc_vec::const_iterator var = sym_stats.begin(); - const symbol &x = *(var->sym); + // The symbol with least degree is our main variable + sym_desc_vec::const_iterator var = sym_stats.begin(); + const symbol &x = *(var->sym); - // Cancel trivial common factor - int ldeg_a = var->ldeg_a; - int ldeg_b = var->ldeg_b; - int min_ldeg = min(ldeg_a, ldeg_b); - if (min_ldeg > 0) { - ex common = power(x, min_ldeg); + // Cancel trivial common factor + int ldeg_a = var->ldeg_a; + int ldeg_b = var->ldeg_b; + int min_ldeg = std::min(ldeg_a,ldeg_b); + if (min_ldeg > 0) { + ex common = power(x, min_ldeg); //std::clog << "trivial common factor " << common << endl; - return gcd((aex / common).expand(), (bex / 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) { + // Try to eliminate variables + if (var->deg_a == 0) { //std::clog << "eliminating variable " << x << " from b" << endl; - ex c = bex.content(x); - ex g = gcd(aex, c, ca, cb, false); - if (cb) - *cb *= bex.unit(x) * bex.primpart(x, c); - return g; - } else if (var->deg_b == 0) { + ex c = bex.content(x); + ex g = gcd(aex, c, ca, cb, false); + if (cb) + *cb *= bex.unit(x) * bex.primpart(x, c); + return g; + } else if (var->deg_b == 0) { //std::clog << "eliminating variable " << x << " from a" << endl; - ex c = aex.content(x); - ex g = gcd(c, bex, ca, cb, false); - if (ca) - *ca *= aex.unit(x) * aex.primpart(x, c); - return g; - } - - ex g; + ex c = aex.content(x); + ex g = gcd(c, bex, ca, cb, false); + if (ca) + *ca *= aex.unit(x) * aex.primpart(x, c); + return g; + } + + ex g; #if 1 - // Try heuristic algorithm first, fall back to PRS if that failed - try { - g = heur_gcd(aex, bex, ca, cb, var); - } catch (gcdheu_failed) { - g = *new ex(fail()); - } - if (is_ex_exactly_of_type(g, fail)) { + // Try heuristic algorithm first, fall back to PRS if that failed + try { + g = heur_gcd(aex, bex, ca, cb, var); + } catch (gcdheu_failed) { + g = *new ex(fail()); + } + if (is_ex_exactly_of_type(g, fail)) { //std::clog << "heuristics failed" << endl; #if STATISTICS heur_gcd_failed++; @@ -1646,13 +1640,13 @@ factored_b: if (cb) *cb = b; } else { - if (ca) - divide(aex, g, *ca, false); - if (cb) - divide(bex, g, *cb, false); + if (ca) + divide(aex, g, *ca, false); + if (cb) + divide(bex, g, *cb, false); } #if 1 - } else { + } else { if (g.is_equal(_ex1())) { // Keep cofactors factored if possible if (ca) @@ -1662,7 +1656,7 @@ factored_b: } } #endif - return g; + return g; } @@ -1675,14 +1669,14 @@ factored_b: * @return the LCM as a new expression */ ex lcm(const ex &a, const ex &b, bool check_args) { - if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) - return lcm(ex_to_numeric(a), ex_to_numeric(b)); - if (check_args && !a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) - throw(std::invalid_argument("lcm: arguments must be polynomials over the rationals")); - - ex ca, cb; - ex g = gcd(a, b, &ca, &cb, false); - return ca * cb * g; + if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) + return lcm(ex_to_numeric(a), ex_to_numeric(b)); + if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) + throw(std::invalid_argument("lcm: arguments must be polynomials over the rationals")); + + ex ca, cb; + ex g = gcd(a, b, &ca, &cb, false); + return ca * cb * g; } @@ -1694,34 +1688,34 @@ ex lcm(const ex &a, const ex &b, bool check_args) // a and b can be multivariate polynomials but they are treated as univariate polynomials in x. static ex univariate_gcd(const ex &a, const ex &b, const symbol &x) { - if (a.is_zero()) - return b; - if (b.is_zero()) - return a; - if (a.is_equal(_ex1()) || b.is_equal(_ex1())) - return _ex1(); - if (is_ex_of_type(a, numeric) && is_ex_of_type(b, numeric)) - return gcd(ex_to_numeric(a), ex_to_numeric(b)); - if (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) - throw(std::invalid_argument("univariate_gcd: arguments must be polynomials over the rationals")); - - // Euclidean algorithm - ex c, d, r; - if (a.degree(x) >= b.degree(x)) { - c = a; - d = b; - } else { - c = b; - d = a; - } - for (;;) { - r = rem(c, d, x, false); - if (r.is_zero()) - break; - c = d; - d = r; - } - return d / d.lcoeff(x); + if (a.is_zero()) + return b; + if (b.is_zero()) + return a; + if (a.is_equal(_ex1()) || b.is_equal(_ex1())) + return _ex1(); + if (is_ex_of_type(a, numeric) && is_ex_of_type(b, numeric)) + return gcd(ex_to_numeric(a), ex_to_numeric(b)); + if (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) + throw(std::invalid_argument("univariate_gcd: arguments must be polynomials over the rationals")); + + // Euclidean algorithm + ex c, d, r; + if (a.degree(x) >= b.degree(x)) { + c = a; + d = b; + } else { + c = b; + d = a; + } + for (;;) { + r = rem(c, d, x, false); + if (r.is_zero()) + break; + c = d; + d = r; + } + return d / d.lcoeff(x); } @@ -1733,27 +1727,27 @@ static ex univariate_gcd(const ex &a, const ex &b, const symbol &x) * @return factored polynomial */ ex sqrfree(const ex &a, const symbol &x) { - int i = 1; - ex res = _ex1(); - ex b = a.diff(x); - ex c = univariate_gcd(a, b, x); - ex w; - if (c.is_equal(_ex1())) { - w = a; - } else { - w = quo(a, c, x); - ex y = quo(b, c, x); - ex z = y - w.diff(x); - while (!z.is_zero()) { - ex g = univariate_gcd(w, z, x); - res *= power(g, i); - w = quo(w, g, x); - y = quo(z, g, x); - z = y - w.diff(x); - i++; - } - } - return res * power(w, i); + int i = 1; + ex res = _ex1(); + ex b = a.diff(x); + ex c = univariate_gcd(a, b, x); + ex w; + if (c.is_equal(_ex1())) { + w = a; + } else { + w = quo(a, c, x); + ex y = quo(b, c, x); + ex z = y - w.diff(x); + while (!z.is_zero()) { + ex g = univariate_gcd(w, z, x); + res *= power(g, i); + w = quo(w, g, x); + y = quo(z, g, x); + z = y - w.diff(x); + i++; + } + } + return res * power(w, i); } @@ -1775,20 +1769,20 @@ ex sqrfree(const ex &a, const symbol &x) * @see ex::normal */ static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst) { - // Expression already in repl_lst? Then return the assigned symbol - for (unsigned i=0; isetflag(status_flags::dynallocated); + return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); } @@ -1825,7 +1819,7 @@ ex basic::normal(lst &sym_lst, lst &repl_lst, int level) const * @see ex::normal */ ex symbol::normal(lst &sym_lst, lst &repl_lst, int level) const { - return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated); + return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated); } @@ -1838,15 +1832,15 @@ ex numeric::normal(lst &sym_lst, lst &repl_lst, int level) const numeric num = numer(); ex numex = num; - if (num.is_real()) { - if (!num.is_integer()) - numex = replace_with_symbol(numex, sym_lst, repl_lst); - } else { // complex - numeric re = num.real(), im = num.imag(); - ex re_ex = re.is_rational() ? re : replace_with_symbol(re, sym_lst, repl_lst); - ex im_ex = im.is_rational() ? im : replace_with_symbol(im, sym_lst, repl_lst); - numex = re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst); - } + if (num.is_real()) { + if (!num.is_integer()) + numex = replace_with_symbol(numex, sym_lst, repl_lst); + } else { // complex + numeric re = num.real(), im = num.imag(); + ex re_ex = re.is_rational() ? re : replace_with_symbol(re, sym_lst, repl_lst); + ex im_ex = im.is_rational() ? im : replace_with_symbol(im, sym_lst, repl_lst); + numex = re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst); + } // Denominator is always a real integer (see numeric::denom()) return (new lst(numex, denom()))->setflag(status_flags::dynallocated); @@ -1859,29 +1853,33 @@ ex numeric::normal(lst &sym_lst, lst &repl_lst, int level) const * @return cancelled fraction {n, d} as a list */ static ex frac_cancel(const ex &n, const ex &d) { - ex num = n; - ex den = d; - numeric pre_factor = _num1(); + ex num = n; + ex den = d; + numeric pre_factor = _num1(); //std::clog << "frac_cancel num = " << num << ", den = " << den << endl; - // Handle special cases where numerator or denominator is 0 - if (num.is_zero()) - return (new lst(_ex0(), _ex1()))->setflag(status_flags::dynallocated); - if (den.expand().is_zero()) - throw(std::overflow_error("frac_cancel: division by zero in frac_cancel")); + // Handle trivial case where denominator is 1 + if (den.is_equal(_ex1())) + return (new lst(num, den))->setflag(status_flags::dynallocated); - // Bring numerator and denominator to Z[X] by multiplying with - // LCM of all coefficients' denominators - numeric num_lcm = lcm_of_coefficients_denominators(num); - numeric den_lcm = lcm_of_coefficients_denominators(den); + // Handle special cases where numerator or denominator is 0 + if (num.is_zero()) + return (new lst(num, _ex1()))->setflag(status_flags::dynallocated); + if (den.expand().is_zero()) + throw(std::overflow_error("frac_cancel: division by zero in frac_cancel")); + + // Bring numerator and denominator to Z[X] by multiplying with + // LCM of all coefficients' denominators + numeric num_lcm = lcm_of_coefficients_denominators(num); + numeric den_lcm = lcm_of_coefficients_denominators(den); num = multiply_lcm(num, num_lcm); den = multiply_lcm(den, den_lcm); - pre_factor = den_lcm / num_lcm; + pre_factor = den_lcm / num_lcm; - // Cancel GCD from numerator and denominator - ex cnum, cden; - if (gcd(num, den, &cnum, &cden, false) != _ex1()) { + // Cancel GCD from numerator and denominator + ex cnum, cden; + if (gcd(num, den, &cnum, &cden, false) != _ex1()) { num = cnum; den = cden; } @@ -1890,7 +1888,7 @@ static ex frac_cancel(const ex &n, const ex &d) // as defined by get_first_symbol() is made positive) const symbol *x; if (get_first_symbol(den, x)) { - GINAC_ASSERT(is_ex_exactly_of_type(den.unit(*x),numeric)); + GINAC_ASSERT(is_ex_exactly_of_type(den.unit(*x),numeric)); if (ex_to_numeric(den.unit(*x)).is_negative()) { num *= _ex_1(); den *= _ex_1(); @@ -1899,7 +1897,7 @@ static ex frac_cancel(const ex &n, const ex &d) // Return result as list //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << endl; - return (new lst(num * pre_factor.numer(), den * pre_factor.denom()))->setflag(status_flags::dynallocated); + return (new lst(num * pre_factor.numer(), den * pre_factor.denom()))->setflag(status_flags::dynallocated); } @@ -1909,78 +1907,56 @@ static ex frac_cancel(const ex &n, const ex &d) ex add::normal(lst &sym_lst, lst &repl_lst, int level) const { if (level == 1) - return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated); + return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _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); - epvector::const_iterator it = seq.begin(), itend = seq.end(); - while (it != itend) { - - // Normalize and expand child - ex n = recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1).expand(); - - // If numerator is a sum, chop into summands - if (is_ex_exactly_of_type(n.op(0), add)) { - epvector::const_iterator bit = ex_to_add(n.op(0)).seq.begin(), bitend = ex_to_add(n.op(0)).seq.end(); - while (bit != bitend) { - o.push_back((new lst(recombine_pair_to_ex(*bit), n.op(1)))->setflag(status_flags::dynallocated)); - bit++; - } - - // The overall_coeff is already normalized (== rational), we just - // split it into numerator and denominator - GINAC_ASSERT(ex_to_numeric(ex_to_add(n.op(0)).overall_coeff).is_rational()); - numeric overall = ex_to_numeric(ex_to_add(n.op(0)).overall_coeff); - o.push_back((new lst(overall.numer(), overall.denom() * n.op(1)))->setflag(status_flags::dynallocated)); - } else - o.push_back(n); - it++; - } - o.push_back(overall_coeff.bp->normal(sym_lst, repl_lst, level-1)); - - // o is now a vector of {numerator, denominator} lists - - // Determine common denominator - ex den = _ex1(); - exvector::const_iterator ait = o.begin(), aitend = o.end(); -//std::clog << "add::normal uses the following summands:\n"; - while (ait != aitend) { -//std::clog << " num = " << ait->op(0) << ", den = " << ait->op(1) << endl; - den = lcm(ait->op(1), den, false); - ait++; - } -//std::clog << " common denominator = " << den << endl; + throw(std::runtime_error("max recursion level reached")); - // Add fractions - if (den.is_equal(_ex1())) { - - // Common denominator is 1, simply add all fractions - exvector num_seq; - for (ait=o.begin(); ait!=aitend; ait++) { - num_seq.push_back(ait->op(0) / ait->op(1)); + // Normalize children and split each one into numerator and denominator + exvector nums, dens; + nums.reserve(seq.size()+1); + dens.reserve(seq.size()+1); + epvector::const_iterator it = seq.begin(), itend = seq.end(); + while (it != itend) { + ex n = recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1); + nums.push_back(n.op(0)); + dens.push_back(n.op(1)); + it++; + } + ex n = overall_coeff.bp->normal(sym_lst, repl_lst, level-1); + nums.push_back(n.op(0)); + dens.push_back(n.op(1)); + GINAC_ASSERT(nums.size() == dens.size()); + + // Now, nums is a vector of all numerators and dens is a vector of + // all denominators +//std::clog << "add::normal uses " << nums.size() << " summands:\n"; + + // Add fractions sequentially + exvector::const_iterator num_it = nums.begin(), num_itend = nums.end(); + exvector::const_iterator den_it = dens.begin(), den_itend = dens.end(); +//std::clog << " num = " << *num_it << ", den = " << *den_it << endl; + ex num = *num_it++, den = *den_it++; + while (num_it != num_itend) { +//std::clog << " num = " << *num_it << ", den = " << *den_it << endl; + ex next_num = *num_it++, next_den = *den_it++; + + // Trivially add sequences of fractions with identical denominators + while ((den_it != den_itend) && next_den.is_equal(*den_it)) { + next_num += *num_it; + num_it++; den_it++; } - return (new lst((new add(num_seq))->setflag(status_flags::dynallocated), den))->setflag(status_flags::dynallocated); - } else { + // Additiion of two fractions, taking advantage of the fact that + // the heuristic GCD algorithm computes the cofactors at no extra cost + ex co_den1, co_den2; + ex g = gcd(den, next_den, &co_den1, &co_den2, false); + num = ((num * co_den2) + (next_num * co_den1)).expand(); + den *= co_den2; // this is the lcm(den, next_den) + } +//std::clog << " common denominator = " << den << endl; - // Perform fractional addition - exvector num_seq; - for (ait=o.begin(); ait!=aitend; ait++) { - ex q; - if (!divide(den, ait->op(1), q, false)) { - // should not happen - throw(std::runtime_error("invalid expression in add::normal, division failed")); - } - num_seq.push_back((ait->op(0) * q).expand()); - } - ex num = (new add(num_seq))->setflag(status_flags::dynallocated); - - // Cancel common factors from num/den - return frac_cancel(num, den); - } + // Cancel common factors from num/den + return frac_cancel(num, den); } @@ -1990,27 +1966,27 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const { if (level == 1) - return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated); + return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); else if (level == -max_recursion_level) - throw(std::runtime_error("max recursion level reached")); + throw(std::runtime_error("max recursion level reached")); - // Normalize children, separate into numerator and denominator + // Normalize children, separate into numerator and denominator ex num = _ex1(); ex den = _ex1(); ex n; - epvector::const_iterator it = seq.begin(), itend = seq.end(); - while (it != itend) { + epvector::const_iterator it = seq.begin(), itend = seq.end(); + while (it != itend) { n = recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1); num *= n.op(0); den *= n.op(1); - it++; - } + it++; + } n = overall_coeff.bp->normal(sym_lst, repl_lst, level-1); num *= n.op(0); den *= n.op(1); // Perform fraction cancellation - return frac_cancel(num, den); + return frac_cancel(num, den); } @@ -2021,69 +1997,69 @@ ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const ex power::normal(lst &sym_lst, lst &repl_lst, int level) const { if (level == 1) - return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated); + return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); else if (level == -max_recursion_level) - throw(std::runtime_error("max recursion level reached")); + throw(std::runtime_error("max recursion level reached")); - // Normalize basis - ex n = basis.bp->normal(sym_lst, repl_lst, level-1); + // Normalize basis and exponent (exponent gets reassembled) + ex n_basis = basis.bp->normal(sym_lst, repl_lst, level-1); + ex n_exponent = exponent.bp->normal(sym_lst, repl_lst, level-1); + n_exponent = n_exponent.op(0) / n_exponent.op(1); - if (exponent.info(info_flags::integer)) { + if (n_exponent.info(info_flags::integer)) { - if (exponent.info(info_flags::positive)) { + if (n_exponent.info(info_flags::positive)) { // (a/b)^n -> {a^n, b^n} - return (new lst(power(n.op(0), exponent), power(n.op(1), exponent)))->setflag(status_flags::dynallocated); + return (new lst(power(n_basis.op(0), n_exponent), power(n_basis.op(1), n_exponent)))->setflag(status_flags::dynallocated); - } else if (exponent.info(info_flags::negative)) { + } else if (n_exponent.info(info_flags::negative)) { // (a/b)^-n -> {b^n, a^n} - return (new lst(power(n.op(1), -exponent), power(n.op(0), -exponent)))->setflag(status_flags::dynallocated); + return (new lst(power(n_basis.op(1), -n_exponent), power(n_basis.op(0), -n_exponent)))->setflag(status_flags::dynallocated); } } else { - if (exponent.info(info_flags::positive)) { + if (n_exponent.info(info_flags::positive)) { // (a/b)^x -> {sym((a/b)^x), 1} - return (new lst(replace_with_symbol(power(n.op(0) / n.op(1), exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); + return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); - } else if (exponent.info(info_flags::negative)) { + } else if (n_exponent.info(info_flags::negative)) { - if (n.op(1).is_equal(_ex1())) { + if (n_basis.op(1).is_equal(_ex1())) { // a^-x -> {1, sym(a^x)} - return (new lst(_ex1(), replace_with_symbol(power(n.op(0), -exponent), sym_lst, repl_lst)))->setflag(status_flags::dynallocated); + return (new lst(_ex1(), replace_with_symbol(power(n_basis.op(0), -n_exponent), sym_lst, repl_lst)))->setflag(status_flags::dynallocated); } else { // (a/b)^-x -> {sym((b/a)^x), 1} - return (new lst(replace_with_symbol(power(n.op(1) / n.op(0), -exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); + return (new lst(replace_with_symbol(power(n_basis.op(1) / n_basis.op(0), -n_exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); } - } else { // exponent not numeric + } else { // n_exponent not numeric // (a/b)^x -> {sym((a/b)^x, 1} - return (new lst(replace_with_symbol(power(n.op(0) / n.op(1), exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); + return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); } - } + } } -/** Implementation of ex::normal() for pseries. It normalizes each coefficient and - * replaces the series by a temporary symbol. +/** Implementation of ex::normal() for pseries. It normalizes each coefficient + * and replaces the series by a temporary symbol. * @see ex::normal */ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const { - epvector new_seq; - new_seq.reserve(seq.size()); - - epvector::const_iterator it = seq.begin(), itend = seq.end(); - while (it != itend) { - new_seq.push_back(expair(it->rest.normal(), it->coeff)); - it++; - } - ex n = pseries(relational(var,point), new_seq); + epvector newseq; + for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) { + ex restexp = i->rest.normal(); + if (!restexp.is_zero()) + newseq.push_back(expair(restexp, i->coeff)); + } + ex n = pseries(relational(var,point), newseq); return (new lst(replace_with_symbol(n, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); } @@ -2110,17 +2086,17 @@ ex relational::normal(lst &sym_lst, lst &repl_lst, int level) const * @return normalized expression */ ex ex::normal(int level) const { - lst sym_lst, repl_lst; + lst sym_lst, repl_lst; - ex e = bp->normal(sym_lst, repl_lst, level); + ex e = bp->normal(sym_lst, repl_lst, level); GINAC_ASSERT(is_ex_of_type(e, lst)); // Re-insert replaced symbols - if (sym_lst.nops() > 0) - e = e.subs(sym_lst, repl_lst); + if (sym_lst.nops() > 0) + e = e.subs(sym_lst, repl_lst); // Convert {numerator, denominator} form back to fraction - return e.op(0) / e.op(1); + return e.op(0) / e.op(1); } /** Numerator of an expression. If the expression is not of the normal form @@ -2131,14 +2107,14 @@ ex ex::normal(int level) const * @return numerator */ ex ex::numer(void) const { - lst sym_lst, repl_lst; + lst sym_lst, repl_lst; - ex e = bp->normal(sym_lst, repl_lst, 0); + ex e = bp->normal(sym_lst, repl_lst, 0); GINAC_ASSERT(is_ex_of_type(e, lst)); // Re-insert replaced symbols - if (sym_lst.nops() > 0) - return e.op(0).subs(sym_lst, repl_lst); + if (sym_lst.nops() > 0) + return e.op(0).subs(sym_lst, repl_lst); else return e.op(0); } @@ -2151,14 +2127,14 @@ ex ex::numer(void) const * @return denominator */ ex ex::denom(void) const { - lst sym_lst, repl_lst; + lst sym_lst, repl_lst; - ex e = bp->normal(sym_lst, repl_lst, 0); + ex e = bp->normal(sym_lst, repl_lst, 0); GINAC_ASSERT(is_ex_of_type(e, lst)); // Re-insert replaced symbols - if (sym_lst.nops() > 0) - return e.op(1).subs(sym_lst, repl_lst); + if (sym_lst.nops() > 0) + return e.op(1).subs(sym_lst, repl_lst); else return e.op(1); } @@ -2178,7 +2154,7 @@ ex basic::to_rational(lst &repl_lst) const * @see ex::to_rational */ ex symbol::to_rational(lst &repl_lst) const { - return *this; + return *this; } @@ -2188,16 +2164,16 @@ ex symbol::to_rational(lst &repl_lst) const * @see ex::to_rational */ ex numeric::to_rational(lst &repl_lst) const { - if (is_real()) { - if (!is_rational()) - return replace_with_symbol(*this, repl_lst); - } else { // complex - numeric re = real(); - numeric im = imag(); - ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst); - ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst); - return re_ex + im_ex * replace_with_symbol(I, repl_lst); - } + if (is_real()) { + if (!is_rational()) + return replace_with_symbol(*this, repl_lst); + } else { // complex + numeric re = real(); + numeric im = imag(); + ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst); + ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst); + return re_ex + im_ex * replace_with_symbol(I, repl_lst); + } return *this; } @@ -2218,17 +2194,17 @@ ex power::to_rational(lst &repl_lst) const * @see ex::to_rational */ ex expairseq::to_rational(lst &repl_lst) const { - epvector s; - s.reserve(seq.size()); - for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) { - s.push_back(split_ex_to_pair(recombine_pair_to_ex(*it).to_rational(repl_lst))); - // s.push_back(combine_ex_with_coeff_to_pair((*it).rest.to_rational(repl_lst), - } - ex oc = overall_coeff.to_rational(repl_lst); - if (oc.info(info_flags::numeric)) - return thisexpairseq(s, overall_coeff); - else s.push_back(combine_ex_with_coeff_to_pair(oc,_ex1())); - return thisexpairseq(s, default_overall_coeff()); + epvector s; + s.reserve(seq.size()); + for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) { + s.push_back(split_ex_to_pair(recombine_pair_to_ex(*it).to_rational(repl_lst))); + // s.push_back(combine_ex_with_coeff_to_pair((*it).rest.to_rational(repl_lst), + } + ex oc = overall_coeff.to_rational(repl_lst); + if (oc.info(info_flags::numeric)) + return thisexpairseq(s, overall_coeff); + else s.push_back(combine_ex_with_coeff_to_pair(oc,_ex1())); + return thisexpairseq(s, default_overall_coeff()); } @@ -2252,6 +2228,4 @@ ex ex::to_rational(lst &repl_lst) const } -#ifndef NO_NAMESPACE_GINAC } // namespace GiNaC -#endif // ndef NO_NAMESPACE_GINAC