X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=db35bb8ab805cddf6161dd0a72c89f669cdbb75f;hp=ee27868165e02e148960407369b8d18da282f180;hb=9eed27a1cf527abf1c5a81dea96081f13a86fcd8;hpb=9f410a18b08213185a6ae473295c9945ebdd9985 diff --git a/ginac/normal.cpp b/ginac/normal.cpp index ee278681..db35bb8a 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-2001 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2003 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 @@ -92,14 +92,14 @@ 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)) { + if (is_a(e)) { x = &ex_to(e); return true; - } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) { + } else if (is_exactly_a(e) || is_exactly_a(e)) { for (unsigned i=0; i(e)) { if (get_first_symbol(e.op(0), x)) return true; } @@ -169,12 +169,12 @@ static void add_symbol(const symbol *s, sym_desc_vec &v) // 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)) { + if (is_a(e)) { add_symbol(&ex_to(e), v); - } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) { + } else if (is_exactly_a(e) || is_exactly_a(e)) { for (unsigned i=0; i(e)) { collect_symbols(e.op(0), v); } } @@ -230,18 +230,18 @@ static numeric lcmcoeff(const ex &e, const numeric &l) { if (e.info(info_flags::rational)) return lcm(ex_to(e).denom(), l); - else if (is_ex_exactly_of_type(e, add)) { + else if (is_exactly_a(e)) { numeric c = _num1; for (unsigned i=0; i(e)) { numeric c = _num1; for (unsigned i=0; i(e)) { + if (is_a(e.op(0))) return l; else return pow(lcmcoeff(e.op(0), l), ex_to(e.op(1))); @@ -268,7 +268,7 @@ static numeric lcm_of_coefficients_denominators(const ex &e) * @param lcm LCM to multiply in */ static ex multiply_lcm(const ex &e, const numeric &lcm) { - if (is_ex_exactly_of_type(e, mul)) { + if (is_exactly_a(e)) { unsigned num = e.nops(); exvector v; v.reserve(num + 1); numeric lcm_accum = _num1; @@ -279,14 +279,14 @@ static ex multiply_lcm(const ex &e, const numeric &lcm) } v.push_back(lcm / lcm_accum); return (new mul(v))->setflag(status_flags::dynallocated); - } else if (is_ex_exactly_of_type(e, add)) { + } else if (is_exactly_a(e)) { unsigned num = e.nops(); exvector v; v.reserve(num); for (unsigned i=0; isetflag(status_flags::dynallocated); - } else if (is_ex_exactly_of_type(e, power)) { - if (is_ex_exactly_of_type(e.op(0), symbol)) + } else if (is_exactly_a(e)) { + if (is_a(e.op(0))) return e * lcm; else return pow(multiply_lcm(e.op(0), lcm.power(ex_to(e.op(1)).inverse())), e.op(1)); @@ -364,7 +364,7 @@ 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)) + if (is_exactly_a(a) && is_exactly_a(b)) return a / b; #if FAST_COMPARE if (a.is_equal(b)) @@ -380,8 +380,8 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args) int bdeg = b.degree(x); int rdeg = r.degree(x); ex blcoeff = b.expand().coeff(x, bdeg); - bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric); - exvector v; v.reserve(rdeg - bdeg + 1); + bool blcoeff_is_numeric = is_exactly_a(blcoeff); + exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0)); while (rdeg >= bdeg) { ex term, rcoeff = r.coeff(x, rdeg); if (blcoeff_is_numeric) @@ -414,8 +414,8 @@ 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)) + if (is_exactly_a(a)) { + if (is_exactly_a(b)) return _ex0; else return a; @@ -434,7 +434,7 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) 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); + bool blcoeff_is_numeric = is_exactly_a(blcoeff); while (rdeg >= bdeg) { ex term, rcoeff = r.coeff(x, rdeg); if (blcoeff_is_numeric) @@ -464,27 +464,27 @@ 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)) + if (is_exactly_a(q)) return a; else return q + rem(numer, denom, x) / denom; } -/** Pseudo-remainder of polynomials a(x) and b(x) in Z[x]. +/** Pseudo-remainder of polynomials a(x) and b(x) in Q[x]. * * @param a first polynomial in x (dividend) * @param b second polynomial in x (divisor) * @param x a and b are polynomials in x * @param check_args check whether a and b are polynomials with rational * coefficients (defaults to "true") - * @return pseudo-remainder of a(x) and b(x) in Z[x] */ + * @return pseudo-remainder of a(x) and b(x) in Q[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)) + if (is_exactly_a(a)) { + if (is_exactly_a(b)) return _ex0; else return b; @@ -523,20 +523,20 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) } -/** Sparse pseudo-remainder of polynomials a(x) and b(x) in Z[x]. +/** Sparse pseudo-remainder of polynomials a(x) and b(x) in Q[x]. * * @param a first polynomial in x (dividend) * @param b second polynomial in x (divisor) * @param x a and b are polynomials in x * @param check_args check whether a and b are polynomials with rational * coefficients (defaults to "true") - * @return sparse pseudo-remainder of a(x) and b(x) in Z[x] */ + * @return sparse pseudo-remainder of a(x) and b(x) in Q[x] */ ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args) { if (b.is_zero()) throw(std::overflow_error("prem: division by zero")); - if (is_ex_exactly_of_type(a, numeric)) { - if (is_ex_exactly_of_type(b, numeric)) + if (is_exactly_a(a)) { + if (is_exactly_a(b)) return _ex0; else return b; @@ -581,18 +581,19 @@ ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args) * @param check_args check whether a and b are polynomials with rational * coefficients (defaults to "true") * @return "true" when exact division succeeds (quotient returned in q), - * "false" otherwise */ + * "false" otherwise (q left untouched) */ bool divide(const ex &a, const ex &b, ex &q, bool check_args) { - q = _ex0; if (b.is_zero()) throw(std::overflow_error("divide: division by zero")); - if (a.is_zero()) + if (a.is_zero()) { + q = _ex0; return true; - if (is_ex_exactly_of_type(b, numeric)) { + } + if (is_exactly_a(b)) { q = a / b; return true; - } else if (is_ex_exactly_of_type(a, numeric)) + } else if (is_exactly_a(a)) return false; #if FAST_COMPARE if (a.is_equal(b)) { @@ -611,13 +612,15 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) // Polynomial long division (recursive) ex r = a.expand(); - if (r.is_zero()) + if (r.is_zero()) { + q = _ex0; return true; + } int bdeg = b.degree(*x); int rdeg = r.degree(*x); ex blcoeff = b.expand().coeff(*x, bdeg); - bool blcoeff_is_numeric = is_ex_exactly_of_type(blcoeff, numeric); - exvector v; v.reserve(rdeg - bdeg + 1); + bool blcoeff_is_numeric = is_exactly_a(blcoeff); + exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0)); while (rdeg >= bdeg) { ex term, rcoeff = r.coeff(*x, rdeg); if (blcoeff_is_numeric) @@ -683,8 +686,8 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite q = a; return true; } - if (is_ex_exactly_of_type(a, numeric)) { - if (is_ex_exactly_of_type(b, numeric)) { + if (is_exactly_a(a)) { + if (is_exactly_a(b)) { q = a / b; return q.info(info_flags::integer); } else @@ -778,7 +781,7 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite int rdeg = adeg; ex eb = b.expand(); ex blcoeff = eb.coeff(*x, bdeg); - exvector v; v.reserve(rdeg - bdeg + 1); + exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0)); while (rdeg >= bdeg) { ex term, rcoeff = r.coeff(*x, rdeg); if (!divide_in_z(rcoeff, blcoeff, term, var+1)) @@ -818,7 +821,7 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite ex ex::unit(const symbol &x) const { ex c = expand().lcoeff(x); - if (is_ex_exactly_of_type(c, numeric)) + if (is_exactly_a(c)) return c < _ex0 ? _ex_1 : _ex1; else { const symbol *y; @@ -841,7 +844,7 @@ ex ex::content(const symbol &x) const { if (is_zero()) return _ex0; - if (is_ex_exactly_of_type(*this, numeric)) + if (is_exactly_a(*this)) return info(info_flags::negative) ? -*this : *this; ex e = expand(); if (e.is_zero()) @@ -877,14 +880,14 @@ ex ex::primpart(const symbol &x) const { if (is_zero()) return _ex0; - if (is_ex_exactly_of_type(*this, numeric)) + if (is_exactly_a(*this)) return _ex1; ex c = content(x); if (c.is_zero()) return _ex0; ex u = unit(x); - if (is_ex_exactly_of_type(c, numeric)) + if (is_exactly_a(c)) return *this / (c * u); else return quo(*this, c * u, x, false); @@ -904,11 +907,11 @@ ex ex::primpart(const symbol &x, const ex &c) const return _ex0; if (c.is_zero()) return _ex0; - if (is_ex_exactly_of_type(*this, numeric)) + if (is_exactly_a(*this)) return _ex1; ex u = unit(x); - if (is_ex_exactly_of_type(c, numeric)) + if (is_exactly_a(c)) return *this / (c * u); else return quo(*this, c * u, x, false); @@ -1107,7 +1110,7 @@ static ex red_gcd(const ex &a, const ex &b, const symbol *x) 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)) + if (is_exactly_a(r)) return gamma; else return gamma * r.primpart(*x); @@ -1182,7 +1185,7 @@ static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator 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)) + if (is_exactly_a(r)) return gamma; else return gamma * r.primpart(x); @@ -1353,7 +1356,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const 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)) { + if (is_exactly_a(a) && is_exactly_a(b)) { numeric g = gcd(ex_to(a), ex_to(b)); if (ca) *ca = ex_to(a) / g; @@ -1391,7 +1394,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const // 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)) { + if (!is_exactly_a(gamma)) { // Reconstruct polynomial from GCD of mapped polynomials ex g = interpolate(gamma, xi, x, maxdeg); @@ -1404,7 +1407,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(lc).is_negative()) + if (is_exactly_a(lc) && ex_to(lc).is_negative()) return -g; else return g; @@ -1417,7 +1420,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(lc).is_negative()) + if (is_exactly_a(lc) && ex_to(lc).is_negative()) return -g; else return g; @@ -1430,7 +1433,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(lc).is_negative()) + if (is_exactly_a(lc) && ex_to(lc).is_negative()) return -g; else return g; @@ -1462,7 +1465,7 @@ 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)) { + if (is_exactly_a(a) && is_exactly_a(b)) { numeric g = gcd(ex_to(a), ex_to(b)); if (ca || cb) { if (g.is_zero()) { @@ -1486,8 +1489,8 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) } // Partially factored cases (to avoid expanding large expressions) - if (is_ex_exactly_of_type(a, mul)) { - if (is_ex_exactly_of_type(b, mul) && b.nops() > a.nops()) + if (is_exactly_a(a)) { + if (is_exactly_a(b) && b.nops() > a.nops()) goto factored_b; factored_a: unsigned num = a.nops(); @@ -1505,8 +1508,8 @@ factored_a: if (cb) *cb = part_b; return (new mul(g))->setflag(status_flags::dynallocated); - } else if (is_ex_exactly_of_type(b, mul)) { - if (is_ex_exactly_of_type(a, mul) && a.nops() > b.nops()) + } else if (is_exactly_a(b)) { + if (is_exactly_a(a) && a.nops() > b.nops()) goto factored_a; factored_b: unsigned num = b.nops(); @@ -1528,9 +1531,9 @@ factored_b: #if FAST_COMPARE // Input polynomials of the form poly^n are sometimes also trivial - if (is_ex_exactly_of_type(a, power)) { + if (is_exactly_a(a)) { ex p = a.op(0); - if (is_ex_exactly_of_type(b, power)) { + if (is_exactly_a(b)) { if (p.is_equal(b.op(0))) { // a = p^n, b = p^m, gcd = p^min(n, m) ex exp_a = a.op(1), exp_b = b.op(1); @@ -1558,7 +1561,7 @@ factored_b: return p; } } - } else if (is_ex_exactly_of_type(b, power)) { + } else if (is_exactly_a(b)) { ex p = b.op(0); if (p.is_equal(a)) { // a = p, b = p^n, gcd = p @@ -1647,7 +1650,7 @@ factored_b: } catch (gcdheu_failed) { g = fail(); } - if (is_ex_exactly_of_type(g, fail)) { + if (is_exactly_a(g)) { //std::clog << "heuristics failed" << std::endl; #if STATISTICS heur_gcd_failed++; @@ -1695,7 +1698,7 @@ 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)) + if (is_exactly_a(a) && is_exactly_a(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")); @@ -1738,14 +1741,45 @@ static exvector sqrfree_yun(const ex &a, const symbol &x) return res; } -/** Compute square-free factorization of multivariate polynomial in Q[X]. + +/** Compute a square-free factorization of a 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 polynomial a in square-free factored form. */ + * @return a square-free factorization of \p a. + * + * \note + * A polynomial \f$p(X) \in C[X]\f$ is said square-free + * if, whenever any two polynomials \f$q(X)\f$ and \f$r(X)\f$ + * are such that + * \f[ + * p(X) = q(X)^2 r(X), + * \f] + * we have \f$q(X) \in C\f$. + * This means that \f$p(X)\f$ has no repeated factors, apart + * eventually from constants. + * Given a polynomial \f$p(X) \in C[X]\f$, we say that the + * decomposition + * \f[ + * p(X) = b \cdot p_1(X)^{a_1} \cdot p_2(X)^{a_2} \cdots p_r(X)^{a_r} + * \f] + * is a square-free factorization of \f$p(X)\f$ if the + * following conditions hold: + * -# \f$b \in C\f$ and \f$b \neq 0\f$; + * -# \f$a_i\f$ is a positive integer for \f$i = 1, \ldots, r\f$; + * -# the degree of the polynomial \f$p_i\f$ is strictly positive + * for \f$i = 1, \ldots, r\f$; + * -# the polynomial \f$\Pi_{i=1}^r p_i(X)\f$ is square-free. + * + * Square-free factorizations need not be unique. For example, if + * \f$a_i\f$ is even, we could change the polynomial \f$p_i(X)\f$ + * into \f$-p_i(X)\f$. + * Observe also that the factors \f$p_i(X)\f$ need not be irreducible + * polynomials. + */ ex sqrfree(const ex &a, const lst &l) { - if (is_a(a) || // algorithm does not trap a==0 + if (is_exactly_a(a) || // algorithm does not trap a==0 is_a(a)) // shortcut return a; @@ -1809,6 +1843,7 @@ ex sqrfree(const ex &a, const lst &l) return result * lcm.inverse(); } + /** Compute square-free partial fraction decomposition of rational function * a(x). * @@ -2187,13 +2222,11 @@ ex power::normal(lst &sym_lst, lst &repl_lst, int level) const // (a/b)^-x -> {sym((b/a)^x), 1} 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 { // n_exponent not numeric - - // (a/b)^x -> {sym((a/b)^x, 1} - 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); } } + + // (a/b)^x -> {sym((a/b)^x, 1} + 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); } @@ -2380,4 +2413,106 @@ ex expairseq::to_rational(lst &repl_lst) const } +/** Remove the common factor in the terms of a sum 'e' by calculating the GCD, + * and multiply it into the expression 'factor' (which needs to be initialized + * to 1, unless you're accumulating factors). */ +static ex find_common_factor(const ex & e, ex & factor, lst & repl) +{ + if (is_exactly_a(e)) { + + unsigned num = e.nops(); + exvector terms; terms.reserve(num); + ex gc; + + // Find the common GCD + for (unsigned i=0; i(x) || is_exactly_a(x)) { + ex f = 1; + x = find_common_factor(x, f, repl); + x *= f; + } + + if (i == 0) + gc = x; + else + gc = gcd(gc, x); + + terms.push_back(x); + } + + if (gc.is_equal(_ex1)) + return e; + + // The GCD is the factor we pull out + factor *= gc; + + // Now divide all terms by the GCD + for (unsigned i=0; i(t)) { + for (unsigned j=0; jsetflag(status_flags::dynallocated); + goto term_done; + } + } + } + + divide(t, gc, x); + t = x; +term_done: ; + } + return (new add(terms))->setflag(status_flags::dynallocated); + + } else if (is_exactly_a(e)) { + + unsigned num = e.nops(); + exvector v; v.reserve(num); + + for (unsigned i=0; isetflag(status_flags::dynallocated); + + } else if (is_exactly_a(e)) { + + ex x = e.to_rational(repl); + if (is_exactly_a(x) && x.op(1).info(info_flags::negative)) + return replace_with_symbol(x, repl); + else + return x; + + } else + return e; +} + + +/** Collect common factors in sums. This converts expressions like + * 'a*(b*x+b*y)' to 'a*b*(x+y)'. */ +ex collect_common_factors(const ex & e) +{ + if (is_exactly_a(e) || is_exactly_a(e)) { + + lst repl; + ex factor = 1; + ex r = find_common_factor(e, factor, repl); + return factor.subs(repl) * r.subs(repl); + + } else + return e; +} + + } // namespace GiNaC