X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=619a2325c03b95543a0c1391f7131c5453fa4ce6;hb=2c9eca6dcf983bbca109ed386d548504f3cdfff4;hp=f0e24aa2c5680179b6d53b4c44152c5bb52ca6bb;hpb=003197fb7b25e822db436a77b70ce9d5ccb82714;p=ginac.git diff --git a/ginac/normal.cpp b/ginac/normal.cpp index f0e24aa2..619a2325 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -23,7 +23,6 @@ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ -#include #include #include @@ -34,21 +33,18 @@ #include "constant.h" #include "expairseq.h" #include "fail.h" -#include "indexed.h" #include "inifcns.h" #include "lst.h" #include "mul.h" -#include "ncmul.h" #include "numeric.h" #include "power.h" #include "relational.h" +#include "matrix.h" #include "pseries.h" #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 @@ -140,8 +136,17 @@ struct sym_desc { /** Maximum of deg_a and deg_b (Used for sorting) */ int max_deg; + /** Maximum number of terms of leading coefficient of symbol in both polynomials */ + int max_lcnops; + /** Commparison operator for sorting */ - bool operator<(const sym_desc &x) const {return max_deg < x.max_deg;} + bool operator<(const sym_desc &x) const + { + if (max_deg == x.max_deg) + return max_lcnops < x.max_lcnops; + else + return max_deg < x.max_deg; + } }; // Vector of sym_desc structures @@ -196,7 +201,8 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) 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_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++; @@ -206,7 +212,7 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) 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++; } @@ -223,7 +229,7 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) static numeric lcmcoeff(const ex &e, const numeric &l) { if (e.info(info_flags::rational)) - return lcm(ex_to_numeric(e).denom(), l); + return lcm(ex_to(e).denom(), l); else if (is_ex_exactly_of_type(e, add)) { numeric c = _num1(); for (unsigned i=0; i(e.op(1))); + } return l; } @@ -274,7 +284,10 @@ static ex multiply_lcm(const ex &e, const numeric &lcm) c += multiply_lcm(e.op(i), lcm); return c; } else if (is_ex_exactly_of_type(e, power)) { - return pow(multiply_lcm(e.op(0), lcm.power(ex_to_numeric(e.op(1)).inverse())), e.op(1)); + if (is_ex_exactly_of_type(e.op(0), symbol)) + return e * lcm; + else + return pow(multiply_lcm(e.op(0), lcm.power(ex_to(e.op(1)).inverse())), e.op(1)); } else return e * lcm; } @@ -309,11 +322,11 @@ numeric add::integer_content(void) const 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); + c = gcd(ex_to(it->coeff), c); it++; } GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - c = gcd(ex_to_numeric(overall_coeff),c); + c = gcd(ex_to(overall_coeff),c); return c; } @@ -328,7 +341,7 @@ numeric mul::integer_content(void) const } #endif // def DO_GINAC_ASSERT GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - return abs(ex_to_numeric(overall_coeff)); + return abs(ex_to(overall_coeff)); } @@ -373,7 +386,7 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args) term = rcoeff / blcoeff; else { if (!divide(rcoeff, blcoeff, term, false)) - return *new ex(fail()); + return (new fail())->setflag(status_flags::dynallocated); } term *= power(x, rdeg - bdeg); q += term; @@ -403,7 +416,7 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) if (is_ex_exactly_of_type(b, numeric)) return _ex0(); else - return b; + return a; } #if FAST_COMPARE if (a.is_equal(b)) @@ -426,7 +439,7 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) term = rcoeff / blcoeff; else { if (!divide(rcoeff, blcoeff, term, false)) - return *new ex(fail()); + return (new fail())->setflag(status_flags::dynallocated); } term *= power(x, rdeg - bdeg); r -= (term * b).expand(); @@ -438,6 +451,24 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) } +/** Decompose rational function a(x)=N(x)/D(x) into P(x)+n(x)/D(x) + * with degree(n, x) < degree(D, x). + * + * @param a rational function in x + * @param x a is a function of x + * @return decomposed function. */ +ex decomp_rational(const ex &a, const symbol &x) +{ + ex nd = numer_denom(a); + ex numer = nd.op(0), denom = nd.op(1); + ex q = quo(numer, denom, x); + if (is_ex_exactly_of_type(q, fail)) + return a; + else + return q + rem(numer, denom, x) / denom; +} + + /** Pseudo-remainder of polynomials a(x) and b(x) in Z[x]. * * @param a first polynomial in x (dividend) @@ -612,9 +643,10 @@ typedef std::pair ex2; typedef std::pair exbool; struct ex2_less { - bool operator() (const ex2 p, const ex2 q) const + 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); + int cmp = p.first.compare(q.first); + return ((cmp<0) || (!(cmp>0) && p.second.compare(q.second)<0)); } }; @@ -1174,6 +1206,8 @@ numeric ex::max_coefficient(void) const return bp->max_coefficient(); } +/** Implementation ex::max_coefficient(). + * @see heur_gcd */ numeric basic::max_coefficient(void) const { return _num1(); @@ -1189,11 +1223,11 @@ 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)); + numeric cur_max = abs(ex_to(overall_coeff)); while (it != itend) { numeric a; GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); - a = abs(ex_to_numeric(it->coeff)); + a = abs(ex_to(it->coeff)); if (a > cur_max) cur_max = a; it++; @@ -1212,7 +1246,7 @@ numeric mul::max_coefficient(void) const } #endif // def DO_GINAC_ASSERT GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - return abs(ex_to_numeric(overall_coeff)); + return abs(ex_to(overall_coeff)); } @@ -1236,11 +1270,7 @@ ex basic::smod(const numeric &xi) const 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 } ex add::smod(const numeric &xi) const @@ -1251,21 +1281,13 @@ ex add::smod(const numeric &xi) const 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 + numeric coeff = GiNaC::smod(ex_to(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)); -#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 + numeric coeff = GiNaC::smod(ex_to(overall_coeff), xi); return (new add(newseq,coeff))->setflag(status_flags::dynallocated); } @@ -1279,13 +1301,9 @@ ex mul::smod(const numeric &xi) const it++; } #endif // def DO_GINAC_ASSERT - mul * mulcopyp=new mul(*this); + 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->overall_coeff = GiNaC::smod(ex_to(overall_coeff),xi); mulcopyp->clearflag(status_flags::evaluated); mulcopyp->clearflag(status_flags::hash_calculated); return mulcopyp->setflag(status_flags::dynallocated); @@ -1333,15 +1351,15 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const // Algorithms only works for non-vanishing input polynomials if (a.is_zero() || b.is_zero()) - return *new ex(fail()); + return (new fail())->setflag(status_flags::dynallocated); // 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)); + numeric g = gcd(ex_to(a), ex_to(b)); if (ca) - *ca = ex_to_numeric(a) / g; + *ca = ex_to(a) / g; if (cb) - *cb = ex_to_numeric(b) / g; + *cb = ex_to(b) / g; return g; } @@ -1367,7 +1385,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const // 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; +//std::clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << std::endl; throw gcdheu_failed(); } @@ -1387,7 +1405,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const 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()) + if (is_ex_exactly_of_type(lc, numeric) && ex_to(lc).is_negative()) return -g; else return g; @@ -1400,7 +1418,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const if (ca) *ca = cp; ex lc = g.lcoeff(x); - if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) + if (is_ex_exactly_of_type(lc, numeric) && ex_to(lc).is_negative()) return -g; else return g; @@ -1413,7 +1431,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const if (cb) *cb = cq; ex lc = g.lcoeff(x); - if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) + if (is_ex_exactly_of_type(lc, numeric) && ex_to(lc).is_negative()) return -g; else return g; @@ -1425,7 +1443,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const // Next evaluation point xi = iquo(xi * isqrt(isqrt(xi)) * numeric(73794), numeric(27011)); } - return *new ex(fail()); + return (new fail())->setflag(status_flags::dynallocated); } @@ -1446,7 +1464,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) // 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)); + numeric g = gcd(ex_to(a), ex_to(b)); if (ca || cb) { if (g.is_zero()) { if (ca) @@ -1455,9 +1473,9 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) *cb = _ex0(); } else { if (ca) - *ca = ex_to_numeric(a) / g; + *ca = ex_to(a) / g; if (cb) - *cb = ex_to_numeric(b) / g; + *cb = ex_to(b) / g; } } return g; @@ -1599,20 +1617,20 @@ factored_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; +//std::clog << "trivial common factor " << common << std::endl; return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common; } // Try to eliminate variables if (var->deg_a == 0) { -//std::clog << "eliminating variable " << x << " from b" << endl; +//std::clog << "eliminating variable " << x << " from b" << std::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) { -//std::clog << "eliminating variable " << x << " from a" << endl; +//std::clog << "eliminating variable " << x << " from a" << std::endl; ex c = aex.content(x); ex g = gcd(c, bex, ca, cb, false); if (ca) @@ -1626,10 +1644,10 @@ factored_b: try { g = heur_gcd(aex, bex, ca, cb, var); } catch (gcdheu_failed) { - g = *new ex(fail()); + g = fail(); } if (is_ex_exactly_of_type(g, fail)) { -//std::clog << "heuristics failed" << endl; +//std::clog << "heuristics failed" << std::endl; #if STATISTICS heur_gcd_failed++; #endif @@ -1677,7 +1695,7 @@ factored_b: 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)); + return lcm(ex_to(a), ex_to(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")); @@ -1691,70 +1709,153 @@ ex lcm(const ex &a, const ex &b, bool check_args) * Square-free factorization */ -// Univariate GCD of polynomials in Q[x] (used internally by sqrfree()). -// 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) +/** Compute square-free factorization of multivariate polynomial a(x) using + * Yun´s algorithm. Used internally by sqrfree(). + * + * @param a multivariate polynomial over Z[X], treated here as univariate + * polynomial in x. + * @param x variable to factor in + * @return vector of factors sorted in ascending degree */ +static exvector sqrfree_yun(const ex &a, 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")); + exvector res; + ex w = a; + ex z = w.diff(x); + ex g = gcd(w, z); + if (g.is_equal(_ex1())) { + res.push_back(a); + return res; + } + ex y; + do { + w = quo(w, g, x); + y = quo(z, g, x); + z = y - w.diff(x); + g = gcd(w, z); + res.push_back(g); + } while (!z.is_zero()); + return res; +} - // Euclidean algorithm - ex c, d, r; - if (a.degree(x) >= b.degree(x)) { - c = a; - d = b; +/** Compute square-free factorization of multivariate polynomial in Q[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. */ +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 + return a; + // If no lst of variables to factorize in was specified we have to + // invent one now. Maybe one can optimize here by reversing the order + // or so, I don't know. + lst args; + if (l.nops()==0) { + sym_desc_vec sdv; + get_symbol_stats(a, _ex0(), sdv); + for (sym_desc_vec::iterator it=sdv.begin(); it!=sdv.end(); ++it) + args.append(*it->sym); } 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); + args = 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(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); + // find the factors + exvector factors = sqrfree_yun(tmp,x); + // construct the next list of symbols with the first element popped + lst newargs; + for (int i=1; i0) { + for (exvector::iterator i=factors.begin(); i!=factors.end(); ++i) + *i = sqrfree(*i, newargs); + } + // Done with recursion, now construct the final result + ex result = _ex1(); + exvector::iterator it = factors.begin(); + for (int p = 1; it!=factors.end(); ++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); + return result * lcm.inverse(); } - -/** Compute square-free factorization of multivariate polynomial a(x) using - * Yun´s algorithm. +/** Compute square-free partial fraction decomposition of rational function + * a(x). * - * @param a multivariate polynomial - * @param x variable to factor in - * @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++; + * @param a rational function over Z[x], treated as univariate polynomial + * in x + * @param x variable to factor in + * @return decomposed rational function */ +ex sqrfree_parfrac(const ex & a, const symbol & x) +{ + // Find numerator and denominator + ex nd = numer_denom(a); + ex numer = nd.op(0), denom = nd.op(1); +//clog << "numer = " << numer << ", denom = " << denom << endl; + + // Convert N(x)/D(x) -> Q(x) + R(x)/D(x), so degree(R) < degree(D) + ex red_poly = quo(numer, denom, x), red_numer = rem(numer, denom, x).expand(); +//clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << endl; + + // Factorize denominator and compute cofactors + exvector yun = sqrfree_yun(denom, x); +//clog << "yun factors: " << exprseq(yun) << endl; + int num_yun = yun.size(); + exvector factor; factor.reserve(num_yun); + exvector cofac; cofac.reserve(num_yun); + for (unsigned i=0; isetflag(status_flags::dynallocated); + if (nops() == 0) + return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); + else { + if (level == 1) + 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")); + else { + normal_map_function map_normal(level - 1); + return (new lst(replace_with_symbol(map(map_normal), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); + } + } } @@ -1864,7 +1985,7 @@ static ex frac_cancel(const ex &n, const ex &d) ex den = d; numeric pre_factor = _num1(); -//std::clog << "frac_cancel num = " << num << ", den = " << den << endl; +//std::clog << "frac_cancel num = " << num << ", den = " << den << std::endl; // Handle trivial case where denominator is 1 if (den.is_equal(_ex1())) @@ -1896,14 +2017,14 @@ static ex frac_cancel(const ex &n, const ex &d) const symbol *x; if (get_first_symbol(den, x)) { GINAC_ASSERT(is_ex_exactly_of_type(den.unit(*x),numeric)); - if (ex_to_numeric(den.unit(*x)).is_negative()) { + if (ex_to(den.unit(*x)).is_negative()) { num *= _ex_1(); den *= _ex_1(); } } // Return result as list -//std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << endl; +//std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << std::endl; return (new lst(num * pre_factor.numer(), den * pre_factor.denom()))->setflag(status_flags::dynallocated); } @@ -1941,10 +2062,10 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const // 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; +//std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl; ex num = *num_it++, den = *den_it++; while (num_it != num_itend) { -//std::clog << " num = " << *num_it << ", den = " << *den_it << endl; +//std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl; ex next_num = *num_it++, next_den = *den_it++; // Trivially add sequences of fractions with identical denominators @@ -1960,7 +2081,7 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const 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; +//std::clog << " common denominator = " << den << std::endl; // Cancel common factors from num/den return frac_cancel(num, den); @@ -2008,77 +2129,69 @@ ex power::normal(lst &sym_lst, lst &repl_lst, int level) const 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 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++; + 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), new_seq); + ex n = pseries(relational(var,point), newseq); return (new lst(replace_with_symbol(n, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); } -/** Implementation of ex::normal() for relationals. It normalizes both sides. - * @see ex::normal */ -ex relational::normal(lst &sym_lst, lst &repl_lst, int level) const -{ - return (new lst(relational(lh.normal(), rh.normal(), o), _ex1()))->setflag(status_flags::dynallocated); -} - - /** Normalization of rational functions. * This function converts an expression to its normal form * "numerator/denominator", where numerator and denominator are (relatively @@ -2106,9 +2219,9 @@ ex ex::normal(int level) const return e.op(0) / e.op(1); } -/** Numerator of an expression. If the expression is not of the normal form - * "numerator/denominator", it is first converted to this form and then the - * numerator is returned. +/** Get numerator of an expression. If the expression is not of the normal + * form "numerator/denominator", it is first converted to this form and + * then the numerator is returned. * * @see ex::normal * @return numerator */ @@ -2126,9 +2239,9 @@ ex ex::numer(void) const return e.op(0); } -/** Denominator of an expression. If the expression is not of the normal form - * "numerator/denominator", it is first converted to this form and then the - * denominator is returned. +/** Get denominator of an expression. If the expression is not of the normal + * form "numerator/denominator", it is first converted to this form and + * then the denominator is returned. * * @see ex::normal * @return denominator */ @@ -2146,6 +2259,26 @@ ex ex::denom(void) const return e.op(1); } +/** Get numerator and denominator of an expression. If the expresison is not + * of the normal form "numerator/denominator", it is first converted to this + * form and then a list [numerator, denominator] is returned. + * + * @see ex::normal + * @return a list [numerator, denominator] */ +ex ex::numer_denom(void) const +{ + lst sym_lst, repl_lst; + + 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.subs(sym_lst, repl_lst); + else + return e; +} + /** Default implementation of ex::to_rational(). It replaces the object with a * temporary symbol. @@ -2235,6 +2368,4 @@ ex ex::to_rational(lst &repl_lst) const } -#ifndef NO_NAMESPACE_GINAC } // namespace GiNaC -#endif // ndef NO_NAMESPACE_GINAC