* computation, square-free factorization and rational function normalization. */
/*
- * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2004 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
/** 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.
*
* @return integer content */
numeric ex::integer_content() const
{
epvector::const_iterator it = seq.begin();
epvector::const_iterator itend = seq.end();
- numeric c = _num0;
+ numeric c = _num0, l = _num1;
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 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
/** 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
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);
return c;
// GCD of all coefficients
- int deg = e.degree(x);
- int ldeg = e.ldegree(x);
+ int deg = r.degree(x);
+ int ldeg = r.ldegree(x);
if (deg == ldeg)
- return e.lcoeff(x) / e.unit(x);
- c = _ex0;
+ return lcoeff * c;
+ 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].
+/** Compute primitive part of a multivariate polynomial in Q[x].
* The product of unit part, content part, and primitive part is the
* polynomial itself.
*
}
-/** 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.
*
}
+/** 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