* 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
#include "numeric.h"
#include "power.h"
#include "relational.h"
+#include "operators.h"
#include "matrix.h"
#include "pseries.h"
#include "symbol.h"
* @return "false" if no symbol was found, "true" otherwise */
static bool get_first_symbol(const ex &e, const symbol *&x)
{
- if (is_exactly_a<symbol>(e)) {
+ if (is_a<symbol>(e)) {
x = &ex_to<symbol>(e);
return true;
} else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
// 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_exactly_a<symbol>(e)) {
+ if (is_a<symbol>(e)) {
add_symbol(&ex_to<symbol>(e), v);
} else if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
for (unsigned i=0; i<e.nops(); i++)
++it;
}
std::sort(v.begin(), v.end());
+
#if 0
std::clog << "Symbols:\n";
it = v.begin(); itend = v.end();
c *= lcmcoeff(e.op(i), _num1);
return lcm(c, l);
} else if (is_exactly_a<power>(e)) {
- if (is_exactly_a<symbol>(e.op(0)))
+ if (is_a<symbol>(e.op(0)))
return l;
else
return pow(lcmcoeff(e.op(0), l), ex_to<numeric>(e.op(1)));
v.push_back(multiply_lcm(e.op(i), lcm));
return (new add(v))->setflag(status_flags::dynallocated);
} else if (is_exactly_a<power>(e)) {
- if (is_exactly_a<symbol>(e.op(0)))
+ 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));
}
-/** 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())
}
-/** 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())
* GCD of multivariate polynomials
*/
-/** Compute GCD of polynomials in Q[X] using the Euclidean algorithm (not
- * really suited for multivariate GCDs). This function is only provided for
- * testing purposes.
- *
- * @param a first multivariate polynomial
- * @param b second multivariate polynomial
- * @param x pointer to symbol (main variable) in which to compute the GCD in
- * @return the GCD as a new expression
- * @see gcd */
-
-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;
- }
-
- // Normalize in Q[x]
- c = c / c.lcoeff(*x);
- d = d / d.lcoeff(*x);
-
- // Euclidean algorithm
- 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;
- d = r;
- }
-}
-
-
-/** Compute GCD of multivariate polynomials using the Euclidean PRS algorithm
- * with pseudo-remainders ("World's Worst GCD Algorithm", staying in Z[X]).
- * This function is only provided for testing purposes.
- *
- * @param a first multivariate polynomial
- * @param b second multivariate polynomial
- * @param x pointer to symbol (main variable) in which to compute the GCD in
- * @return the GCD as a new expression
- * @see gcd */
-
-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;
- }
-
- // Calculate GCD of contents
- ex gamma = gcd(c.content(*x), d.content(*x), NULL, NULL, false);
-
- // Euclidean algorithm with pseudo-remainders
- 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;
- d = r;
- }
-}
-
-
-/** Compute GCD of multivariate polynomials using the primitive Euclidean
- * PRS algorithm (complete content removal at each step). This function is
- * only provided for testing purposes.
- *
- * @param a first multivariate polynomial
- * @param b second multivariate polynomial
- * @param x pointer to symbol (main variable) in which to compute the GCD in
- * @return the GCD as a new expression
- * @see gcd */
-
-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
- ex r;
- for (;;) {
-//std::clog << " d = " << d << endl;
- r = prem(c, d, *x, false);
- if (r.is_zero())
- return gamma * d;
- c = d;
- d = r.primpart(*x);
- }
-}
-
-
-/** Compute GCD of multivariate polynomials using the reduced PRS algorithm.
- * This function is only provided for testing purposes.
- *
- * @param a first multivariate polynomial
- * @param b second multivariate polynomial
- * @param x pointer to symbol (main variable) in which to compute the GCD in
- * @return the GCD as a new expression
- * @see gcd */
-
-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
-//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_exactly_a<numeric>(r))
- return gamma;
- else
- return gamma * r.primpart(*x);
- }
-
- ri = c.expand().lcoeff(*x);
- delta = cdeg - ddeg;
- }
-}
-
-
/** Compute GCD of multivariate polynomials using the subresultant PRS
* algorithm. This function is used internally by gcd().
*
static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var)
{
-//std::clog << "sr_gcd(" << a << "," << b << ")\n";
#if STATISTICS
sr_gcd_called++;
#endif
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;
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;
-//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);
}
// Next element of subresultant sequence
-//std::clog << " calculating next subresultant...\n";
ri = c.expand().lcoeff(x);
if (delta == 1)
psi = ri;
* @exception gcdheu_failed() */
static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var)
{
-//std::clog << "heur_gcd(" << a << "," << b << ")\n";
#if STATISTICS
heur_gcd_called++;
#endif
// 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 << std::endl;
throw gcdheu_failed();
}
else
return g;
}
-#if 0
- cp = interpolate(cp, xi, x);
- if (divide_in_z(cp, p, g, var)) {
- if (divide_in_z(g, q, cb ? *cb : dummy, var)) {
- g *= gc;
- if (ca)
- *ca = cp;
- ex lc = g.lcoeff(x);
- if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
- return -g;
- else
- return g;
- }
- }
- cq = interpolate(cq, xi, x);
- if (divide_in_z(cq, q, g, var)) {
- if (divide_in_z(g, p, ca ? *ca : dummy, var)) {
- g *= gc;
- if (cb)
- *cb = cq;
- ex lc = g.lcoeff(x);
- if (is_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
- return -g;
- else
- return g;
- }
- }
-#endif
}
// Next evaluation point
* @return the GCD as a new expression */
ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args)
{
-//std::clog << "gcd(" << a << "," << b << ")\n";
#if STATISTICS
gcd_called++;
#endif
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 << 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" << 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" << std::endl;
ex c = aex.content(x);
ex g = gcd(c, bex, ca, cb, false);
if (ca)
return g;
}
- ex g;
-#if 1
// Try heuristic algorithm first, fall back to PRS if that failed
+ ex g;
try {
g = heur_gcd(aex, bex, ca, cb, var);
} catch (gcdheu_failed) {
g = fail();
}
if (is_exactly_a<fail>(g)) {
-//std::clog << "heuristics failed" << std::endl;
#if STATISTICS
heur_gcd_failed++;
#endif
-#endif
-// g = heur_gcd(aex, bex, ca, cb, var);
-// g = eu_gcd(aex, bex, &x);
-// g = euprem_gcd(aex, bex, &x);
-// g = peu_gcd(aex, bex, &x);
-// g = red_gcd(aex, bex, &x);
g = sr_gcd(aex, bex, var);
if (g.is_equal(_ex1)) {
// Keep cofactors factored if possible
if (cb)
divide(bex, g, *cb, false);
}
-#if 1
} else {
if (g.is_equal(_ex1)) {
// Keep cofactors factored if possible
*cb = b;
}
}
-#endif
+
return g;
}
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).
*
/** Create a symbol for replacing the expression "e" (or return a previously
* assigned symbol). An expression of the form "symbol == expression" is added
* to repl_lst and the symbol is returned.
- * @see basic::to_rational */
+ * @see basic::to_rational
+ * @see basic::to_polynomial */
static ex replace_with_symbol(const ex &e, lst &repl_lst)
{
// Expression already in repl_lst? Then return the assigned symbol
// (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);
}
/** Rationalization of non-rational functions.
- * This function converts a general expression to a rational polynomial
+ * This function converts a general expression to a rational function
* by replacing all non-rational subexpressions (like non-rational numbers,
* non-integer powers or functions like sin(), cos() etc.) to temporary
* symbols. This makes it possible to use functions like gcd() and divide()
*
* @param repl_lst collects a list of all temporary symbols and their replacements
* @return rationalized expression */
+ex ex::to_rational(lst &repl_lst) const
+{
+ return bp->to_rational(repl_lst);
+}
+
+ex ex::to_polynomial(lst &repl_lst) const
+{
+ return bp->to_polynomial(repl_lst);
+}
+
+
+/** Default implementation of ex::to_rational(). This replaces the object with
+ * a temporary symbol. */
ex basic::to_rational(lst &repl_lst) const
{
return replace_with_symbol(*this, repl_lst);
}
+ex basic::to_polynomial(lst &repl_lst) const
+{
+ return replace_with_symbol(*this, repl_lst);
+}
+
/** Implementation of ex::to_rational() for symbols. This returns the
* unmodified symbol. */
return *this;
}
+/** Implementation of ex::to_polynomial() for symbols. This returns the
+ * unmodified symbol. */
+ex symbol::to_polynomial(lst &repl_lst) const
+{
+ return *this;
+}
+
/** Implementation of ex::to_rational() for a numeric. It splits complex
* numbers into re+I*im and replaces I and non-rational real numbers with a
return *this;
}
+/** Implementation of ex::to_polynomial() for a numeric. It splits complex
+ * numbers into re+I*im and replaces I and non-integer real numbers with a
+ * temporary symbol. */
+ex numeric::to_polynomial(lst &repl_lst) const
+{
+ if (is_real()) {
+ if (!is_integer())
+ return replace_with_symbol(*this, repl_lst);
+ } else { // complex
+ numeric re = real();
+ numeric im = imag();
+ ex re_ex = re.is_integer() ? re : replace_with_symbol(re, repl_lst);
+ ex im_ex = im.is_integer() ? im : replace_with_symbol(im, repl_lst);
+ return re_ex + im_ex * replace_with_symbol(I, repl_lst);
+ }
+ return *this;
+}
+
/** Implementation of ex::to_rational() for powers. It replaces non-integer
* powers by temporary symbols. */
return replace_with_symbol(*this, repl_lst);
}
+/** Implementation of ex::to_polynomial() for powers. It replaces non-posint
+ * powers by temporary symbols. */
+ex power::to_polynomial(lst &repl_lst) const
+{
+ if (exponent.info(info_flags::posint))
+ return power(basis.to_rational(repl_lst), exponent);
+ else
+ return replace_with_symbol(*this, repl_lst);
+}
+
/** Implementation of ex::to_rational() for expairseqs. */
ex expairseq::to_rational(lst &repl_lst) const
return thisexpairseq(s, default_overall_coeff());
}
+/** Implementation of ex::to_polynomial() for expairseqs. */
+ex expairseq::to_polynomial(lst &repl_lst) const
+{
+ epvector s;
+ s.reserve(seq.size());
+ epvector::const_iterator i = seq.begin(), end = seq.end();
+ while (i != end) {
+ s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_polynomial(repl_lst)));
+ ++i;
+ }
+ ex oc = overall_coeff.to_polynomial(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());
+}
+
+
+/** 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_polynomial(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)) {
+
+ return e.to_polynomial(repl);
+
+ } 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