* computation, square-free factorization and rational function normalization. */
/*
- * GiNaC Copyright (C) 1999-2002 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
* @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<symbol>(e)) {
x = &ex_to<symbol>(e);
return true;
- } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) {
+ } else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
for (unsigned i=0; i<e.nops(); i++)
if (get_first_symbol(e.op(i), x))
return true;
- } else if (is_ex_exactly_of_type(e, power)) {
+ } else if (is_exactly_a<power>(e)) {
if (get_first_symbol(e.op(0), x))
return true;
}
// 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<symbol>(e)) {
add_symbol(&ex_to<symbol>(e), v);
- } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) {
+ } else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
for (unsigned i=0; i<e.nops(); i++)
collect_symbols(e.op(i), v);
- } else if (is_ex_exactly_of_type(e, power)) {
+ } else if (is_exactly_a<power>(e)) {
collect_symbols(e.op(0), v);
}
}
{
if (e.info(info_flags::rational))
return lcm(ex_to<numeric>(e).denom(), l);
- else if (is_ex_exactly_of_type(e, add)) {
+ else if (is_exactly_a<add>(e)) {
numeric c = _num1;
for (unsigned i=0; i<e.nops(); i++)
c = lcmcoeff(e.op(i), c);
return lcm(c, l);
- } else if (is_ex_exactly_of_type(e, mul)) {
+ } else if (is_exactly_a<mul>(e)) {
numeric c = _num1;
for (unsigned i=0; i<e.nops(); i++)
c *= lcmcoeff(e.op(i), _num1);
return lcm(c, l);
- } else if (is_ex_exactly_of_type(e, power)) {
- if (is_ex_exactly_of_type(e.op(0), symbol))
+ } else if (is_exactly_a<power>(e)) {
+ if (is_a<symbol>(e.op(0)))
return l;
else
return pow(lcmcoeff(e.op(0), l), ex_to<numeric>(e.op(1)));
* @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<mul>(e)) {
unsigned num = e.nops();
exvector v; v.reserve(num + 1);
numeric lcm_accum = _num1;
}
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<add>(e)) {
unsigned num = e.nops();
exvector v; v.reserve(num);
for (unsigned i=0; i<num; i++)
v.push_back(multiply_lcm(e.op(i), lcm));
return (new add(v))->setflag(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<power>(e)) {
+ if (is_a<symbol>(e.op(0)))
return e * lcm;
else
return pow(multiply_lcm(e.op(0), lcm.power(ex_to<numeric>(e.op(1)).inverse())), e.op(1));
{
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<numeric>(a) && is_exactly_a<numeric>(b))
return a / b;
#if FAST_COMPARE
if (a.is_equal(b))
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<numeric>(blcoeff);
exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
while (rdeg >= bdeg) {
ex term, rcoeff = r.coeff(x, rdeg);
{
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<numeric>(a)) {
+ if (is_exactly_a<numeric>(b))
return _ex0;
else
return a;
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<numeric>(blcoeff);
while (rdeg >= bdeg) {
ex term, rcoeff = r.coeff(x, rdeg);
if (blcoeff_is_numeric)
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<fail>(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<numeric>(a)) {
+ if (is_exactly_a<numeric>(b))
return _ex0;
else
return b;
}
-/** 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<numeric>(a)) {
+ if (is_exactly_a<numeric>(b))
return _ex0;
else
return b;
q = _ex0;
return true;
}
- if (is_ex_exactly_of_type(b, numeric)) {
+ if (is_exactly_a<numeric>(b)) {
q = a / b;
return true;
- } else if (is_ex_exactly_of_type(a, numeric))
+ } else if (is_exactly_a<numeric>(a))
return false;
#if FAST_COMPARE
if (a.is_equal(b)) {
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<numeric>(blcoeff);
exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
while (rdeg >= bdeg) {
ex term, rcoeff = r.coeff(*x, rdeg);
q = a;
return true;
}
- if (is_ex_exactly_of_type(a, numeric)) {
- if (is_ex_exactly_of_type(b, numeric)) {
+ if (is_exactly_a<numeric>(a)) {
+ if (is_exactly_a<numeric>(b)) {
q = a / b;
return q.info(info_flags::integer);
} else
ex ex::unit(const symbol &x) const
{
ex c = expand().lcoeff(x);
- if (is_ex_exactly_of_type(c, numeric))
+ if (is_exactly_a<numeric>(c))
return c < _ex0 ? _ex_1 : _ex1;
else {
const symbol *y;
{
if (is_zero())
return _ex0;
- if (is_ex_exactly_of_type(*this, numeric))
+ if (is_exactly_a<numeric>(*this))
return info(info_flags::negative) ? -*this : *this;
ex e = expand();
if (e.is_zero())
{
if (is_zero())
return _ex0;
- if (is_ex_exactly_of_type(*this, numeric))
+ if (is_exactly_a<numeric>(*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<numeric>(c))
return *this / (c * u);
else
return quo(*this, c * u, x, false);
return _ex0;
if (c.is_zero())
return _ex0;
- if (is_ex_exactly_of_type(*this, numeric))
+ if (is_exactly_a<numeric>(*this))
return _ex1;
ex u = unit(x);
- if (is_ex_exactly_of_type(c, numeric))
+ if (is_exactly_a<numeric>(c))
return *this / (c * u);
else
return quo(*this, c * u, x, 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))
+ if (is_exactly_a<numeric>(r))
return gamma;
else
return gamma * r.primpart(*x);
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<numeric>(r))
return gamma;
else
return gamma * r.primpart(x);
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<numeric>(a) && is_exactly_a<numeric>(b)) {
numeric g = gcd(ex_to<numeric>(a), ex_to<numeric>(b));
if (ca)
*ca = ex_to<numeric>(a) / g;
// 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<fail>(gamma)) {
// Reconstruct polynomial from GCD of mapped polynomials
ex g = interpolate(gamma, xi, x, maxdeg);
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_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
return -g;
else
return g;
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_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
return -g;
else
return g;
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_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
return -g;
else
return g;
#endif
// GCD of numerics -> CLN
- if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) {
+ if (is_exactly_a<numeric>(a) && is_exactly_a<numeric>(b)) {
numeric g = gcd(ex_to<numeric>(a), ex_to<numeric>(b));
if (ca || cb) {
if (g.is_zero()) {
}
// 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<mul>(a)) {
+ if (is_exactly_a<mul>(b) && b.nops() > a.nops())
goto factored_b;
factored_a:
unsigned num = a.nops();
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<mul>(b)) {
+ if (is_exactly_a<mul>(a) && a.nops() > b.nops())
goto factored_a;
factored_b:
unsigned num = b.nops();
#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<power>(a)) {
ex p = a.op(0);
- if (is_ex_exactly_of_type(b, power)) {
+ if (is_exactly_a<power>(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);
return p;
}
}
- } else if (is_ex_exactly_of_type(b, power)) {
+ } else if (is_exactly_a<power>(b)) {
ex p = b.op(0);
if (p.is_equal(a)) {
// a = p, b = p^n, gcd = p
} catch (gcdheu_failed) {
g = fail();
}
- if (is_ex_exactly_of_type(g, fail)) {
+ if (is_exactly_a<fail>(g)) {
//std::clog << "heuristics failed" << std::endl;
#if STATISTICS
heur_gcd_failed++;
* @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<numeric>(a) && is_exactly_a<numeric>(b))
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"));
return res;
}
+
/** Compute a square-free factorization of a multivariate polynomial in Q[X].
*
* @param a multivariate polynomial over Q[X]
*/
ex sqrfree(const ex &a, const lst &l)
{
- if (is_a<numeric>(a) || // algorithm does not trap a==0
+ if (is_exactly_a<numeric>(a) || // algorithm does not trap a==0
is_a<symbol>(a)) // shortcut
return a;
return result * lcm.inverse();
}
+
/** Compute square-free partial fraction decomposition of rational function
* a(x).
*
// Make denominator unit normal (i.e. coefficient of first symbol
// as defined by get_first_symbol() is made positive)
- const symbol *x;
- if (get_first_symbol(den, x)) {
- GINAC_ASSERT(is_exactly_a<numeric>(den.unit(*x)));
- if (ex_to<numeric>(den.unit(*x)).is_negative()) {
+ if (is_exactly_a<numeric>(den)) {
+ if (ex_to<numeric>(den).is_negative()) {
num *= _ex_1;
den *= _ex_1;
}
+ } else {
+ const symbol *x;
+ if (get_first_symbol(den, x)) {
+ GINAC_ASSERT(is_exactly_a<numeric>(den.unit(*x)));
+ if (ex_to<numeric>(den.unit(*x)).is_negative()) {
+ num *= _ex_1;
+ den *= _ex_1;
+ }
+ }
}
// Return result as list
// (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);
}
}
+/** 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<add>(e)) {
+
+ unsigned num = e.nops();
+ exvector terms; terms.reserve(num);
+ ex gc;
+
+ // Find the common GCD
+ for (unsigned i=0; i<num; i++) {
+ ex x = e.op(i).to_rational(repl);
+
+ if (is_exactly_a<add>(x) || is_exactly_a<mul>(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<num; i++) {
+ ex x;
+
+ // Try to avoid divide() because it expands the polynomial
+ ex &t = terms[i];
+ if (is_exactly_a<mul>(t)) {
+ for (unsigned j=0; j<t.nops(); j++) {
+ if (t.op(j).is_equal(gc)) {
+ exvector v; v.reserve(t.nops());
+ for (unsigned k=0; k<t.nops(); k++) {
+ if (k == j)
+ v.push_back(_ex1);
+ else
+ v.push_back(t.op(k));
+ }
+ t = (new mul(v))->setflag(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<mul>(e)) {
+
+ unsigned num = e.nops();
+ exvector v; v.reserve(num);
+
+ for (unsigned i=0; i<num; i++)
+ v.push_back(find_common_factor(e.op(i), factor, repl));
+
+ return (new mul(v))->setflag(status_flags::dynallocated);
+
+ } else if (is_exactly_a<power>(e)) {
+
+ ex x = e.to_rational(repl);
+ if (is_exactly_a<power>(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<add>(e) || is_exactly_a<mul>(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
git://www.ginac.de / ginac.git / blobdiff