* computation, square-free factorization and rational function normalization. */
/*
- * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2002 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
static struct _stat_print {
_stat_print() {}
~_stat_print() {
- cout << "gcd() called " << gcd_called << " times\n";
- cout << "sr_gcd() called " << sr_gcd_called << " times\n";
- cout << "heur_gcd() called " << heur_gcd_called << " times\n";
- cout << "heur_gcd() failed " << heur_gcd_failed << " times\n";
+ std::cout << "gcd() called " << gcd_called << " times\n";
+ std::cout << "sr_gcd() called " << sr_gcd_called << " times\n";
+ std::cout << "heur_gcd() called " << heur_gcd_called << " times\n";
+ std::cout << "heur_gcd() failed " << heur_gcd_failed << " times\n";
}
} stat_print;
#endif
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);
+ exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
while (rdeg >= bdeg) {
ex term, rcoeff = r.coeff(x, rdeg);
if (blcoeff_is_numeric)
}
-/** 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())
* @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)) {
q = a / b;
return true;
// 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);
+ exvector v; v.reserve(std::max(rdeg - bdeg + 1, 0));
while (rdeg >= bdeg) {
ex term, rcoeff = r.coeff(*x, rdeg);
if (blcoeff_is_numeric)
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))
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 polynomail 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 <EM>square-free</EM>
+ * 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 <EM>square-free factorization</EM> 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_ex_of_type(a,numeric) || // algorithm does not trap a==0
- is_ex_of_type(a,symbol)) // shortcut
+ if (is_a<numeric>(a) || // algorithm does not trap a==0
+ is_a<symbol>(a)) // shortcut
return a;
// If no lst of variables to factorize in was specified we have to
// Find the symbol to factor in at this stage
if (!is_ex_of_type(args.op(0), symbol))
throw (std::runtime_error("sqrfree(): invalid factorization variable"));
- const symbol x = ex_to<symbol>(args.op(0));
+ const symbol &x = ex_to<symbol>(args.op(0));
// convert the argument from something in Q[X] to something in Z[X]
- numeric lcm = lcm_of_coefficients_denominators(a);
- ex tmp = multiply_lcm(a,lcm);
+ const numeric lcm = lcm_of_coefficients_denominators(a);
+ const ex tmp = multiply_lcm(a,lcm);
// find the factors
exvector factors = sqrfree_yun(tmp,x);
lst newargs = args;
newargs.remove_first();
- // recurse down the factors in remaining vars
+ // recurse down the factors in remaining variables
if (newargs.nops()>0) {
- exvector::iterator i = factors.begin(), end = factors.end();
- while (i != end) {
+ exvector::iterator i = factors.begin();
+ while (i != factors.end()) {
*i = sqrfree(*i, newargs);
++i;
}
for (int p = 1; it!=itend; ++it, ++p)
result *= power(*it, p);
- // Yun's algorithm does not account for constant factors. (For
- // univariate polynomials it works only in the monic case.) We can
- // correct this by inserting what has been lost back into the result:
- result = result * quo(tmp, result, x);
+ // Yun's algorithm does not account for constant factors. (For univariate
+ // polynomials it works only in the monic case.) We can correct this by
+ // inserting what has been lost back into the result. For completeness
+ // we'll also have to recurse down that factor in the remaining variables.
+ if (newargs.nops()>0)
+ result *= sqrfree(quo(tmp, result, x), newargs);
+ else
+ result *= quo(tmp, result, x);
+
+ // Put in the reational overall factor again and return
return result * lcm.inverse();
}
}
+/** 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)
+{
+ if (is_a<add>(e)) {
+
+ unsigned num = e.nops();
+ exvector terms; terms.reserve(num);
+ lst repl;
+ ex gc;
+
+ // Find the common GCD
+ for (unsigned i=0; i<num; i++) {
+ ex x = e.op(i).to_rational(repl);
+ if (is_a<add>(x) || is_a<mul>(x)) {
+ ex f = 1;
+ x = find_common_factor(x, f);
+ 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.subs(repl);
+
+ // 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_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).subs(repl);
+ goto term_done;
+ }
+ }
+ }
+
+ divide(t, gc, x);
+ t = x.subs(repl);
+term_done: ;
+ }
+ return (new add(terms))->setflag(status_flags::dynallocated);
+
+ } else if (is_a<mul>(e)) {
+
+ exvector v;
+ for (unsigned i=0; i<e.nops(); i++)
+ v.push_back(find_common_factor(e.op(i), factor));
+ return (new mul(v))->setflag(status_flags::dynallocated);
+
+ } 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_a<add>(e) || is_a<mul>(e)) {
+
+ ex factor = 1;
+ ex r = find_common_factor(e, factor);
+ return factor * r;
+
+ } else
+ return e;
+}
+
+
/*
* Normal form of rational functions
*/
git://www.ginac.de / ginac.git / blobdiff