* computation, square-free factorization and rational function normalization. */
/*
- * GiNaC Copyright (C) 1999-2018 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2020 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
#if 0
std::clog << "Symbols:\n";
- it = v.begin(); itend = v.end();
+ auto it = v.begin(), itend = v.end();
while (it != itend) {
- std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << endl;
- std::clog << " lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << endl;
+ std::clog << " " << it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << ", max_lcnops=" << it->max_lcnops << std::endl;
+ std::clog << " lcoeff_a=" << a.lcoeff(it->sym) << ", lcoeff_b=" << b.lcoeff(it->sym) << std::endl;
++it;
}
#endif
*cb = b;
return _ex1;
}
- // move symbols which are contained only in one of the polynomials
- // to the end:
+ // move symbol contained only in one of the polynomials to the end:
rotate(sym_stats.begin(), vari, sym_stats.end());
sym_desc_vec::const_iterator var = sym_stats.begin();
if (cb)
*cb = b;
return _ex1;
- // XXX: do I need to check for p_gcd = -1?
}
// there are common factors:
ex z = w.diff(x);
ex g = gcd(w, z);
if (g.is_zero()) {
- return epvector{};
+ // manifest zero or hidden zero
+ return {};
}
if (g.is_equal(_ex1)) {
- return epvector{expair(a, _ex1)};
+ // w(x) and w'(x) share no factors: w(x) is square-free
+ return {expair(a, _ex1)};
}
- epvector results;
- ex exponent = _ex0;
+
+ epvector factors;
+ ex i = 0; // exponent
do {
w = quo(w, g, x);
if (w.is_zero()) {
- return results;
+ // hidden zero
+ break;
}
z = quo(z, g, x) - w.diff(x);
- exponent = exponent + 1;
+ i += 1;
if (w.is_equal(x)) {
// shortcut for x^n with n ∈ ℕ
- exponent += quo(z, w.diff(x), x);
- results.push_back(expair(w, exponent));
+ i += quo(z, w.diff(x), x);
+ factors.push_back(expair(w, i));
break;
}
g = gcd(w, z);
if (!g.is_equal(_ex1)) {
- results.push_back(expair(g, exponent));
+ factors.push_back(expair(g, i));
}
} while (!z.is_zero());
+
+ // correct for lost factor
+ // (being based on GCDs, Yun's algorithm only finds factors up to a unit)
+ const ex lost_factor = quo(a, mul{factors}, x);
+ if (lost_factor.is_equal(_ex1)) {
+ // trivial lost factor
+ return factors;
+ }
+ if (!factors.empty() && factors[0].coeff.is_equal(1)) {
+ // multiply factor^1 with lost_factor
+ factors[0].rest *= lost_factor;
+ return factors;
+ }
+ // no factor^1: prepend lost_factor^1 to the results
+ epvector results = {expair(lost_factor, 1)};
+ std::move(factors.begin(), factors.end(), std::back_inserter(results));
return results;
}
// convert the argument from something in Q[X] to something in Z[X]
const numeric lcm = lcm_of_coefficients_denominators(a);
- const ex tmp = multiply_lcm(a,lcm);
+ const ex tmp = multiply_lcm(a, lcm);
// find the factors
epvector factors = sqrfree_yun(tmp, x);
+ if (factors.empty()) {
+ // the polynomial was a hidden zero
+ return _ex0;
+ }
// remove symbol x and proceed recursively with the remaining symbols
args.remove_first();
}
// Done with recursion, now construct the final result
- ex result = _ex1;
- for (auto & it : factors)
- result *= pow(it.rest, it.coeff);
-
- // 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 (args.nops()>0)
- result *= sqrfree(quo(tmp, result, x), args);
- else
- result *= quo(tmp, result, x);
+ ex result = mul(factors);
// Put in the rational overall factor again and return
- return result * lcm.inverse();
+ return result * lcm.inverse();
}
// Find numerator and denominator
ex nd = numer_denom(a);
ex numer = nd.op(0), denom = nd.op(1);
-//clog << "numer = " << numer << ", denom = " << denom << endl;
+//std::clog << "numer = " << numer << ", denom = " << denom << std::endl;
// Convert N(x)/D(x) -> Q(x) + R(x)/D(x), so degree(R) < degree(D)
ex red_poly = quo(numer, denom, x), red_numer = rem(numer, denom, x).expand();
-//clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << endl;
+//std::clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << std::endl;
// Factorize denominator and compute cofactors
epvector yun = sqrfree_yun(denom, x);
}
}
size_t num_factors = factor.size();
-//clog << "factors : " << exprseq(factor) << endl;
-//clog << "cofactors: " << exprseq(cofac) << endl;
+//std::clog << "factors : " << exprseq(factor) << std::endl;
+//std::clog << "cofactors: " << exprseq(cofac) << std::endl;
// Construct coefficient matrix for decomposition
int max_denom_deg = denom.degree(x);
sys(i, j) = cofac[j].coeff(x, i);
rhs(i, 0) = red_numer.coeff(x, i);
}
-//clog << "coeffs: " << sys << endl;
-//clog << "rhs : " << rhs << endl;
+//std::clog << "coeffs: " << sys << std::endl;
+//std::clog << "rhs : " << rhs << std::endl;
// Solve resulting linear system
matrix vars(num_factors, 1);
x *= f;
}
- if (i == 0)
+ if (gc.is_zero())
gc = x;
else
gc = gcd(gc, x);
if (gc.is_equal(_ex1))
return e;
+ if (gc.is_zero())
+ return _ex0;
+
// The GCD is the factor we pull out
factor *= gc;
git://www.ginac.de / ginac.git / blobdiff