* computation, square-free factorization and rational function normalization. */
/*
- * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2005 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
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
- * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include <algorithm>
#endif
-/** Return pointer to first symbol found in expression. Due to GiNaC´s
+/** Return pointer to first symbol found in expression. Due to GiNaC's
* internal ordering of terms, it may not be obvious which symbol this
* function returns for a given expression.
*
if (e.info(info_flags::rational))
return lcm(ex_to<numeric>(e).denom(), l);
else if (is_exactly_a<add>(e)) {
- numeric c = _num1;
+ numeric c = *_num1_p;
for (size_t i=0; i<e.nops(); i++)
c = lcmcoeff(e.op(i), c);
return lcm(c, l);
} else if (is_exactly_a<mul>(e)) {
- numeric c = _num1;
+ numeric c = *_num1_p;
for (size_t i=0; i<e.nops(); i++)
- c *= lcmcoeff(e.op(i), _num1);
+ c *= lcmcoeff(e.op(i), *_num1_p);
return lcm(c, l);
} else if (is_exactly_a<power>(e)) {
if (is_a<symbol>(e.op(0)))
* @return LCM of denominators of coefficients */
static numeric lcm_of_coefficients_denominators(const ex &e)
{
- return lcmcoeff(e, _num1);
+ return lcmcoeff(e, *_num1_p);
}
/** Bring polynomial from Q[X] to Z[X] by multiplying in the previously
if (is_exactly_a<mul>(e)) {
size_t num = e.nops();
exvector v; v.reserve(num + 1);
- numeric lcm_accum = _num1;
+ numeric lcm_accum = *_num1_p;
for (size_t i=0; i<num; i++) {
- numeric op_lcm = lcmcoeff(e.op(i), _num1);
+ numeric op_lcm = lcmcoeff(e.op(i), *_num1_p);
v.push_back(multiply_lcm(e.op(i), op_lcm));
lcm_accum *= op_lcm;
}
/** Compute the integer content (= GCD of all numeric coefficients) of an
- * expanded polynomial.
+ * expanded polynomial. For a polynomial with rational coefficients, this
+ * returns g/l where g is the GCD of the coefficients' numerators and l
+ * is the LCM of the coefficients' denominators.
*
- * @param e expanded polynomial
* @return integer content */
numeric ex::integer_content() const
{
numeric basic::integer_content() const
{
- return _num1;
+ return *_num1_p;
}
numeric numeric::integer_content() const
{
epvector::const_iterator it = seq.begin();
epvector::const_iterator itend = seq.end();
- numeric c = _num0;
+ numeric c = *_num0_p, l = *_num1_p;
while (it != itend) {
GINAC_ASSERT(!is_exactly_a<numeric>(it->rest));
GINAC_ASSERT(is_exactly_a<numeric>(it->coeff));
- c = gcd(ex_to<numeric>(it->coeff), c);
+ c = gcd(ex_to<numeric>(it->coeff).numer(), c);
+ l = lcm(ex_to<numeric>(it->coeff).denom(), l);
it++;
}
GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
- c = gcd(ex_to<numeric>(overall_coeff),c);
- return c;
+ c = gcd(ex_to<numeric>(overall_coeff).numer(), c);
+ l = lcm(ex_to<numeric>(overall_coeff).denom(), l);
+ return c/l;
}
numeric mul::integer_content() const
// Compute values at evaluation points 0..adeg
vector<numeric> alpha; alpha.reserve(adeg + 1);
exvector u; u.reserve(adeg + 1);
- numeric point = _num0;
+ numeric point = *_num0_p;
ex c;
for (i=0; i<=adeg; i++) {
ex bs = b.subs(x == point, subs_options::no_pattern);
while (bs.is_zero()) {
- point += _num1;
+ point += *_num1_p;
bs = b.subs(x == point, subs_options::no_pattern);
}
if (!divide_in_z(a.subs(x == point, subs_options::no_pattern), bs, c, var+1))
return false;
alpha.push_back(point);
u.push_back(c);
- point += _num1;
+ point += *_num1_p;
}
// Compute inverses
vector<numeric> rcp; rcp.reserve(adeg + 1);
- rcp.push_back(_num0);
+ rcp.push_back(*_num0_p);
for (k=1; k<=adeg; k++) {
numeric product = alpha[k] - alpha[0];
for (i=1; i<k; i++)
*/
/** Compute unit part (= sign of leading coefficient) of a multivariate
- * polynomial in Z[x]. The product of unit part, content part, and primitive
+ * polynomial in Q[x]. The product of unit part, content part, and primitive
* part is the polynomial itself.
*
- * @param x variable in which to compute the unit part
+ * @param x main variable
* @return unit part
- * @see ex::content, ex::primpart */
+ * @see ex::content, ex::primpart, ex::unitcontprim */
ex ex::unit(const ex &x) const
{
ex c = expand().lcoeff(x);
if (is_exactly_a<numeric>(c))
- return c < _ex0 ? _ex_1 : _ex1;
+ return c.info(info_flags::negative) ?_ex_1 : _ex1;
else {
ex y;
if (get_first_symbol(c, y))
/** Compute content part (= unit normal GCD of all coefficients) of a
- * multivariate polynomial in Z[x]. The product of unit part, content part,
+ * multivariate polynomial in Q[x]. The product of unit part, content part,
* and primitive part is the polynomial itself.
*
- * @param x variable in which to compute the content part
+ * @param x main variable
* @return content part
- * @see ex::unit, ex::primpart */
+ * @see ex::unit, ex::primpart, ex::unitcontprim */
ex ex::content(const ex &x) const
{
- if (is_zero())
- return _ex0;
if (is_exactly_a<numeric>(*this))
return info(info_flags::negative) ? -*this : *this;
+
ex e = expand();
if (e.is_zero())
return _ex0;
- // First, try the integer content
+ // First, divide out the integer content (which we can calculate very efficiently).
+ // If the leading coefficient of the quotient is an integer, we are done.
ex c = e.integer_content();
ex r = e / c;
- ex lcoeff = r.lcoeff(x);
+ int deg = r.degree(x);
+ ex lcoeff = r.coeff(x, deg);
if (lcoeff.info(info_flags::integer))
return c;
// GCD of all coefficients
- int deg = e.degree(x);
- int ldeg = e.ldegree(x);
+ int ldeg = r.ldegree(x);
if (deg == ldeg)
- return e.lcoeff(x) / e.unit(x);
- c = _ex0;
+ return lcoeff * c / lcoeff.unit(x);
+ ex cont = _ex0;
for (int i=ldeg; i<=deg; i++)
- c = gcd(e.coeff(x, i), c, NULL, NULL, false);
- return c;
+ cont = gcd(r.coeff(x, i), cont, NULL, NULL, false);
+ return cont * c;
}
-/** Compute primitive part of a multivariate polynomial in Z[x].
- * The product of unit part, content part, and primitive part is the
- * polynomial itself.
+/** Compute primitive part of a multivariate polynomial in Q[x]. The result
+ * will be a unit-normal polynomial with a content part of 1. The product
+ * of unit part, content part, and primitive part is the polynomial itself.
*
- * @param x variable in which to compute the primitive part
+ * @param x main variable
* @return primitive part
- * @see ex::unit, ex::content */
+ * @see ex::unit, ex::content, ex::unitcontprim */
ex ex::primpart(const ex &x) const
{
- if (is_zero())
- return _ex0;
- if (is_exactly_a<numeric>(*this))
- return _ex1;
-
- ex c = content(x);
- if (c.is_zero())
- return _ex0;
- ex u = unit(x);
- if (is_exactly_a<numeric>(c))
- return *this / (c * u);
- else
- return quo(*this, c * u, x, false);
+ // We need to compute the unit and content anyway, so call unitcontprim()
+ ex u, c, p;
+ unitcontprim(x, u, c, p);
+ return p;
}
-/** Compute primitive part of a multivariate polynomial in Z[x] when the
+/** Compute primitive part of a multivariate polynomial in Q[x] when the
* content part is already known. This function is faster in computing the
* primitive part than the previous function.
*
- * @param x variable in which to compute the primitive part
+ * @param x main variable
* @param c previously computed content part
* @return primitive part */
ex ex::primpart(const ex &x, const ex &c) const
{
- if (is_zero())
- return _ex0;
- if (c.is_zero())
+ if (is_zero() || c.is_zero())
return _ex0;
if (is_exactly_a<numeric>(*this))
return _ex1;
+ // Divide by unit and content to get primitive part
ex u = unit(x);
if (is_exactly_a<numeric>(c))
return *this / (c * u);
}
+/** Compute unit part, content part, and primitive part of a multivariate
+ * polynomial in Q[x]. The product of the three parts is the polynomial
+ * itself.
+ *
+ * @param x main variable
+ * @param u unit part (returned)
+ * @param c content part (returned)
+ * @param p primitive part (returned)
+ * @see ex::unit, ex::content, ex::primpart */
+void ex::unitcontprim(const ex &x, ex &u, ex &c, ex &p) const
+{
+ // Quick check for zero (avoid expanding)
+ if (is_zero()) {
+ u = _ex1;
+ c = p = _ex0;
+ return;
+ }
+
+ // Special case: input is a number
+ if (is_exactly_a<numeric>(*this)) {
+ if (info(info_flags::negative)) {
+ u = _ex_1;
+ c = abs(ex_to<numeric>(*this));
+ } else {
+ u = _ex1;
+ c = *this;
+ }
+ p = _ex1;
+ return;
+ }
+
+ // Expand input polynomial
+ ex e = expand();
+ if (e.is_zero()) {
+ u = _ex1;
+ c = p = _ex0;
+ return;
+ }
+
+ // Compute unit and content
+ u = unit(x);
+ c = content(x);
+
+ // Divide by unit and content to get primitive part
+ if (c.is_zero()) {
+ p = _ex0;
+ return;
+ }
+ if (is_exactly_a<numeric>(c))
+ p = *this / (c * u);
+ else
+ p = quo(e, c * u, x, false);
+}
+
+
/*
* GCD of multivariate polynomials
*/
/** Return maximum (absolute value) coefficient of a polynomial.
* This function is used internally by heur_gcd().
*
- * @param e expanded multivariate polynomial
* @return maximum coefficient
* @see heur_gcd */
numeric ex::max_coefficient() const
* @see heur_gcd */
numeric basic::max_coefficient() const
{
- return _num1;
+ return *_num1_p;
}
numeric numeric::max_coefficient() const
numeric mq = q.max_coefficient();
numeric xi;
if (mp > mq)
- xi = mq * _num2 + _num2;
+ xi = mq * (*_num2_p) + (*_num2_p);
else
- xi = mp * _num2 + _num2;
+ xi = mp * (*_num2_p) + (*_num2_p);
// 6 tries maximum
for (int t=0; t<6; t++) {
ex dummy;
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_exactly_a<numeric>(lc) && ex_to<numeric>(lc).is_negative())
- return -g;
- else
- return g;
+ return g;
}
}
/** Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X)
- * and b(X) in Z[X].
+ * and b(X) in Z[X]. Optionally also compute the cofactors of a and b,
+ * defined by a = ca * gcd(a, b) and b = cb * gcd(a, b).
*
* @param a first multivariate polynomial
* @param b second multivariate polynomial
+ * @param ca pointer to expression that will receive the cofactor of a, or NULL
+ * @param cb pointer to expression that will receive the cofactor of b, or NULL
* @param check_args check whether a and b are polynomials with rational
* coefficients (defaults to "true")
* @return the GCD as a new expression */
// Input polynomials of the form poly^n are sometimes also trivial
if (is_exactly_a<power>(a)) {
ex p = a.op(0);
+ const ex& exp_a = a.op(1);
if (is_exactly_a<power>(b)) {
- if (p.is_equal(b.op(0))) {
+ ex pb = b.op(0);
+ const ex& exp_b = b.op(1);
+ if (p.is_equal(pb)) {
// a = p^n, b = p^m, gcd = p^min(n, m)
- ex exp_a = a.op(1), exp_b = b.op(1);
if (exp_a < exp_b) {
if (ca)
*ca = _ex1;
*cb = _ex1;
return power(p, exp_b);
}
- }
+ } else {
+ ex p_co, pb_co;
+ ex p_gcd = gcd(p, pb, &p_co, &pb_co, check_args);
+ if (p_gcd.is_equal(_ex1)) {
+ // a(x) = p(x)^n, b(x) = p_b(x)^m, gcd (p, p_b) = 1 ==>
+ // gcd(a,b) = 1
+ if (ca)
+ *ca = a;
+ if (cb)
+ *cb = b;
+ return _ex1;
+ // XXX: do I need to check for p_gcd = -1?
+ } else {
+ // there are common factors:
+ // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==>
+ // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m
+ if (exp_a < exp_b) {
+ return power(p_gcd, exp_a)*
+ gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false);
+ } else {
+ return power(p_gcd, exp_b)*
+ gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false);
+ }
+ } // p_gcd.is_equal(_ex1)
+ } // p.is_equal(pb)
+
} else {
if (p.is_equal(b)) {
// a = p^n, b = p, gcd = p
if (cb)
*cb = _ex1;
return p;
+ }
+
+ ex p_co, bpart_co;
+ ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
+
+ if (p_gcd.is_equal(_ex1)) {
+ // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
+ if (ca)
+ *ca = a;
+ if (cb)
+ *cb = b;
+ return _ex1;
+ } else {
+ // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
+ return p_gcd*gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false);
}
- }
+ } // is_exactly_a<power>(b)
+
} else if (is_exactly_a<power>(b)) {
ex p = b.op(0);
if (p.is_equal(a)) {
*cb = power(p, b.op(1) - 1);
return p;
}
+
+ ex p_co, apart_co;
+ const ex& exp_b(b.op(1));
+ ex p_gcd = gcd(a, p, &apart_co, &p_co, false);
+ if (p_gcd.is_equal(_ex1)) {
+ // b=p(x)^n, gcd(a, p) = 1 ==> gcd(a, b) == 1
+ if (ca)
+ *ca = a;
+ if (cb)
+ *cb = b;
+ return _ex1;
+ } else {
+ // there are common factors:
+ // a(x) = g(x) A(x), b(x) = g(x)^n B(x)^n ==> gcd = g(x) gcd(g(x)^(n-1) A(x)^n, B(x))
+
+ return p_gcd*gcd(apart_co, power(p_gcd, exp_b-1)*power(p_co, exp_b), ca, cb, false);
+ } // p_gcd.is_equal(_ex1)
}
#endif
// Try to eliminate variables
if (var->deg_a == 0) {
- ex c = bex.content(x);
- ex g = gcd(aex, c, ca, cb, false);
+ 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.unit(x) * bex.primpart(x, c);
+ *cb *= bex_u * bex_p;
return g;
} else if (var->deg_b == 0) {
- ex c = aex.content(x);
- ex g = gcd(c, bex, ca, cb, false);
+ 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);
if (ca)
- *ca *= aex.unit(x) * aex.primpart(x, c);
+ *ca *= aex_u * aex_p;
return g;
}
*/
/** Compute square-free factorization of multivariate polynomial a(x) using
- * Yun´s algorithm. Used internally by sqrfree().
+ * Yun's algorithm. Used internally by sqrfree().
*
* @param a multivariate polynomial over Z[X], treated here as univariate
* polynomial in 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
+ * @param l lst of variables to factor in, may be left empty for autodetection
* @return a square-free factorization of \p a.
*
* \note
}
/** 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.
+ * assigned symbol). The symbol and expression are appended to repl, and the
+ * symbol is returned.
* @see basic::to_rational
* @see basic::to_polynomial */
-static ex replace_with_symbol(const ex & e, lst & repl_lst)
+static ex replace_with_symbol(const ex & e, exmap & repl)
{
- // Expression already in repl_lst? Then return the assigned symbol
- for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
- if (it->op(1).is_equal(e))
- return it->op(0);
+ // Expression already replaced? Then return the assigned symbol
+ for (exmap::const_iterator it = repl.begin(); it != repl.end(); ++it)
+ if (it->second.is_equal(e))
+ return it->first;
// Otherwise create new symbol and add to list, taking care that the
- // replacement expression doesn't itself contain symbols from the repl_lst,
+ // replacement expression doesn't itself contain symbols from repl,
// because subs() is not recursive
ex es = (new symbol)->setflag(status_flags::dynallocated);
- ex e_replaced = e.subs(repl_lst, subs_options::no_pattern);
- repl_lst.append(es == e_replaced);
+ ex e_replaced = e.subs(repl, subs_options::no_pattern);
+ repl.insert(std::make_pair(es, e_replaced));
return es;
}
{
ex num = n;
ex den = d;
- numeric pre_factor = _num1;
+ numeric pre_factor = *_num1_p;
//std::clog << "frac_cancel num = " << num << ", den = " << den << std::endl;
* on non-rational functions by applying to_rational() on the arguments,
* calling the desired function and re-substituting the temporary symbols
* in the result. To make the last step possible, all temporary symbols and
- * their associated expressions are collected in the list specified by the
- * repl_lst parameter in the form {symbol == expression}, ready to be passed
- * as an argument to ex::subs().
+ * their associated expressions are collected in the map specified by the
+ * repl parameter, ready to be passed as an argument to ex::subs().
*
- * @param repl_lst collects a list of all temporary symbols and their replacements
+ * @param repl collects all temporary symbols and their replacements
* @return rationalized expression */
-ex ex::to_rational(lst &repl_lst) const
+ex ex::to_rational(exmap & repl) const
{
- return bp->to_rational(repl_lst);
+ return bp->to_rational(repl);
+}
+
+// GiNaC 1.1 compatibility function
+ex ex::to_rational(lst & repl_lst) const
+{
+ // Convert lst to exmap
+ exmap m;
+ for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
+ m.insert(std::make_pair(it->op(0), it->op(1)));
+
+ ex ret = bp->to_rational(m);
+
+ // Convert exmap back to lst
+ repl_lst.remove_all();
+ for (exmap::const_iterator it = m.begin(); it != m.end(); ++it)
+ repl_lst.append(it->first == it->second);
+
+ return ret;
}
-ex ex::to_polynomial(lst &repl_lst) const
+ex ex::to_polynomial(exmap & repl) const
{
- return bp->to_polynomial(repl_lst);
+ return bp->to_polynomial(repl);
}
+// GiNaC 1.1 compatibility function
+ex ex::to_polynomial(lst & repl_lst) const
+{
+ // Convert lst to exmap
+ exmap m;
+ for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
+ m.insert(std::make_pair(it->op(0), it->op(1)));
+
+ ex ret = bp->to_polynomial(m);
+
+ // Convert exmap back to lst
+ repl_lst.remove_all();
+ for (exmap::const_iterator it = m.begin(); it != m.end(); ++it)
+ repl_lst.append(it->first == it->second);
+
+ return ret;
+}
/** Default implementation of ex::to_rational(). This replaces the object with
* a temporary symbol. */
-ex basic::to_rational(lst &repl_lst) const
+ex basic::to_rational(exmap & repl) const
{
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
}
-ex basic::to_polynomial(lst &repl_lst) const
+ex basic::to_polynomial(exmap & repl) const
{
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
}
/** Implementation of ex::to_rational() for symbols. This returns the
* unmodified symbol. */
-ex symbol::to_rational(lst &repl_lst) const
+ex symbol::to_rational(exmap & repl) const
{
return *this;
}
/** Implementation of ex::to_polynomial() for symbols. This returns the
* unmodified symbol. */
-ex symbol::to_polynomial(lst &repl_lst) const
+ex symbol::to_polynomial(exmap & repl) 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
* temporary symbol. */
-ex numeric::to_rational(lst &repl_lst) const
+ex numeric::to_rational(exmap & repl) const
{
if (is_real()) {
if (!is_rational())
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
} else { // complex
numeric re = real();
numeric im = imag();
- ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst);
- ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst);
- return re_ex + im_ex * replace_with_symbol(I, repl_lst);
+ ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl);
+ ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl);
+ return re_ex + im_ex * replace_with_symbol(I, repl);
}
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
+ex numeric::to_polynomial(exmap & repl) const
{
if (is_real()) {
if (!is_integer())
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
} 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);
+ ex re_ex = re.is_integer() ? re : replace_with_symbol(re, repl);
+ ex im_ex = im.is_integer() ? im : replace_with_symbol(im, repl);
+ return re_ex + im_ex * replace_with_symbol(I, repl);
}
return *this;
}
/** Implementation of ex::to_rational() for powers. It replaces non-integer
* powers by temporary symbols. */
-ex power::to_rational(lst &repl_lst) const
+ex power::to_rational(exmap & repl) const
{
if (exponent.info(info_flags::integer))
- return power(basis.to_rational(repl_lst), exponent);
+ return power(basis.to_rational(repl), exponent);
else
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
}
/** Implementation of ex::to_polynomial() for powers. It replaces non-posint
* powers by temporary symbols. */
-ex power::to_polynomial(lst &repl_lst) const
+ex power::to_polynomial(exmap & repl) const
{
if (exponent.info(info_flags::posint))
- return power(basis.to_rational(repl_lst), exponent);
+ return power(basis.to_rational(repl), exponent);
else
- return replace_with_symbol(*this, repl_lst);
+ return replace_with_symbol(*this, repl);
}
/** Implementation of ex::to_rational() for expairseqs. */
-ex expairseq::to_rational(lst &repl_lst) const
+ex expairseq::to_rational(exmap & repl) 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_rational(repl_lst)));
+ s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_rational(repl)));
++i;
}
- ex oc = overall_coeff.to_rational(repl_lst);
+ ex oc = overall_coeff.to_rational(repl);
if (oc.info(info_flags::numeric))
return thisexpairseq(s, overall_coeff);
else
}
/** Implementation of ex::to_polynomial() for expairseqs. */
-ex expairseq::to_polynomial(lst &repl_lst) const
+ex expairseq::to_polynomial(exmap & repl) 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)));
+ s.push_back(split_ex_to_pair(recombine_pair_to_ex(*i).to_polynomial(repl)));
++i;
}
- ex oc = overall_coeff.to_polynomial(repl_lst);
+ ex oc = overall_coeff.to_polynomial(repl);
if (oc.info(info_flags::numeric))
return thisexpairseq(s, overall_coeff);
else
/** 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)
+static ex find_common_factor(const ex & e, ex & factor, exmap & repl)
{
if (is_exactly_a<add>(e)) {
{
if (is_exactly_a<add>(e) || is_exactly_a<mul>(e)) {
- lst repl;
+ exmap repl;
ex factor = 1;
ex r = find_common_factor(e, factor, repl);
return factor.subs(repl, subs_options::no_pattern) * r.subs(repl, subs_options::no_pattern);
}
+/** Resultant of two expressions e1,e2 with respect to symbol s.
+ * Method: Compute determinant of Sylvester matrix of e1,e2,s. */
+ex resultant(const ex & e1, const ex & e2, const ex & s)
+{
+ const ex ee1 = e1.expand();
+ const ex ee2 = e2.expand();
+ if (!ee1.info(info_flags::polynomial) ||
+ !ee2.info(info_flags::polynomial))
+ throw(std::runtime_error("resultant(): arguments must be polynomials"));
+
+ const int h1 = ee1.degree(s);
+ const int l1 = ee1.ldegree(s);
+ const int h2 = ee2.degree(s);
+ const int l2 = ee2.ldegree(s);
+
+ const int msize = h1 + h2;
+ matrix m(msize, msize);
+
+ for (int l = h1; l >= l1; --l) {
+ const ex e = ee1.coeff(s, l);
+ for (int k = 0; k < h2; ++k)
+ m(k, k+h1-l) = e;
+ }
+ for (int l = h2; l >= l2; --l) {
+ const ex e = ee2.coeff(s, l);
+ for (int k = 0; k < h1; ++k)
+ m(k+h2, k+h2-l) = e;
+ }
+
+ return m.determinant();
+}
+
+
} // namespace GiNaC
git://www.ginac.de / ginac.git / blobdiff