X-Git-Url: https://ginac.de/ginac.git/static/gitweb.css/ginac.git?a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=d8252abb02742fd4376fea26eb0988e140118acb;hb=44f3fac94db304f8a39745f425c7b831b7eec6ec;hp=7b90030de6b0fb595b781b9f3c7a532ad4e6c015;hpb=1c85bdb25fe8fed52caf39e5ec35d6d14e6d837c;p=ginac.git diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 7b90030d..d8252abb 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-2005 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2006 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 @@ -614,6 +614,72 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) if (!get_first_symbol(a, x) && !get_first_symbol(b, x)) throw(std::invalid_argument("invalid expression in divide()")); + // Try to avoid expanding partially factored expressions. + if (is_exactly_a(b)) { + // Divide sequentially by each term + ex rem_new, rem_old = a; + for (size_t i=0; i < b.nops(); i++) { + if (! divide(rem_old, b.op(i), rem_new, false)) + return false; + rem_old = rem_new; + } + q = rem_new; + return true; + } else if (is_exactly_a(b)) { + const ex& bb(b.op(0)); + int exp_b = ex_to(b.op(1)).to_int(); + ex rem_new, rem_old = a; + for (int i=exp_b; i>0; i--) { + if (! divide(rem_old, bb, rem_new, false)) + return false; + rem_old = rem_new; + } + q = rem_new; + return true; + } + + if (is_exactly_a(a)) { + // Divide sequentially each term. If some term in a is divisible + // by b we are done... and if not, we can't really say anything. + size_t i; + ex rem_i; + bool divisible_p = false; + for (i=0; i < a.nops(); ++i) { + if (divide(a.op(i), b, rem_i, false)) { + divisible_p = true; + break; + } + } + if (divisible_p) { + exvector resv; + resv.reserve(a.nops()); + for (size_t j=0; j < a.nops(); j++) { + if (j==i) + resv.push_back(rem_i); + else + resv.push_back(a.op(j)); + } + q = (new mul(resv))->setflag(status_flags::dynallocated); + return true; + } + } else if (is_exactly_a(a)) { + // The base itself might be divisible by b, in that case we don't + // need to expand a + const ex& ab(a.op(0)); + int a_exp = ex_to(a.op(1)).to_int(); + ex rem_i; + if (divide(ab, b, rem_i, false)) { + q = rem_i*power(ab, a_exp - 1); + return true; + } + for (int i=2; i < a_exp; i++) { + if (divide(power(ab, i), b, rem_i, false)) { + q = rem_i*power(ab, a_exp - i); + return true; + } + } // ... so we *really* need to expand expression. + } + // Polynomial long division (recursive) ex r = a.expand(); if (r.is_zero()) { @@ -1523,11 +1589,51 @@ factored_b: } } + if (is_exactly_a(aex)) { + numeric bcont = bex.integer_content(); + numeric g = gcd(ex_to(aex), bcont); + if (ca) + *ca = ex_to(aex)/g; + if (cb) + *cb = bex/g; + return g; + } + + if (is_exactly_a(bex)) { + numeric acont = aex.integer_content(); + numeric g = gcd(ex_to(bex), acont); + if (ca) + *ca = aex/g; + if (cb) + *cb = ex_to(bex)/g; + return g; + } + // Gather symbol statistics sym_desc_vec sym_stats; get_symbol_stats(a, b, sym_stats); - // The symbol with least degree is our main variable + // The symbol with least degree which is contained in both polynomials + // is our main variable + sym_desc_vec::iterator vari = sym_stats.begin(); + while ((vari != sym_stats.end()) && + (((vari->ldeg_b == 0) && (vari->deg_b == 0)) || + ((vari->ldeg_a == 0) && (vari->deg_a == 0)))) + vari++; + + // No common symbols at all, just return 1: + if (vari == sym_stats.end()) { + // N.B: keep cofactors factored + if (ca) + *ca = a; + if (cb) + *cb = b; + return _ex1; + } + // move symbols which contained only in one of the polynomials + // to the end: + rotate(sym_stats.begin(), vari, sym_stats.end()); + sym_desc_vec::const_iterator var = sym_stats.begin(); const ex &x = var->sym; @@ -1541,14 +1647,14 @@ factored_b: } // Try to eliminate variables - if (var->deg_a == 0) { + if (var->deg_a == 0 && var->deg_b != 0 ) { ex bex_u, bex_c, bex_p; bex.unitcontprim(x, bex_u, bex_c, bex_p); ex g = gcd(aex, bex_c, ca, cb, false); if (cb) *cb *= bex_u * bex_p; return g; - } else if (var->deg_b == 0) { + } else if (var->deg_b == 0 && var->deg_a != 0) { ex aex_u, aex_c, aex_p; aex.unitcontprim(x, aex_u, aex_c, aex_p); ex g = gcd(aex_c, bex, ca, cb, false); @@ -2389,7 +2495,16 @@ ex power::to_polynomial(exmap & repl) const if (exponent.info(info_flags::posint)) return power(basis.to_rational(repl), exponent); else if (exponent.info(info_flags::negint)) - return power(replace_with_symbol(power(basis, _ex_1), repl), -exponent); + { + ex basis_pref = collect_common_factors(basis); + if (is_exactly_a(basis_pref) || is_exactly_a(basis_pref)) { + // (A*B)^n will be automagically transformed to A^n*B^n + ex t = power(basis_pref, exponent); + return t.to_polynomial(repl); + } + else + return power(replace_with_symbol(power(basis, _ex_1), repl), -exponent); + } else return replace_with_symbol(*this, repl); } @@ -2447,7 +2562,7 @@ static ex find_common_factor(const ex & e, ex & factor, exmap & repl) for (size_t i=0; i(x) || is_exactly_a(x)) { + if (is_exactly_a(x) || is_exactly_a(x) || is_a(x)) { ex f = 1; x = find_common_factor(x, f, repl); x *= f; @@ -2507,7 +2622,7 @@ term_done: ; } else if (is_exactly_a(e)) { const ex e_exp(e.op(1)); - if (e_exp.info(info_flags::posint)) { + if (e_exp.info(info_flags::integer)) { ex eb = e.op(0).to_polynomial(repl); ex factor_local(_ex1); ex pre_res = find_common_factor(eb, factor_local, repl);