* Yun's algorithm. Used internally by sqrfree().
*
* @param a multivariate polynomial over Z[X], treated here as univariate
- * polynomial in x.
+ * polynomial in x (needs not be expanded).
* @param x variable to factor in
* @return vector of factors sorted in ascending degree */
static exvector sqrfree_yun(const ex &a, const symbol &x)
ex w = a;
ex z = w.diff(x);
ex g = gcd(w, z);
+ if (g.is_zero()) {
+ return res;
+ }
if (g.is_equal(_ex1)) {
res.push_back(a);
return res;
ex y;
do {
w = quo(w, g, x);
+ if (w.is_zero()) {
+ return res;
+ }
y = quo(z, g, x);
z = y - w.diff(x);
g = gcd(w, z);
/** Compute a square-free factorization of a multivariate polynomial in Q[X].
*
- * @param a multivariate polynomial over Q[X]
+ * @param a multivariate polynomial over Q[X] (needs not be expanded)
* @param l lst of variables to factor in, may be left empty for autodetection
* @return a square-free factorization of \p a.
*
*/
ex sqrfree(const ex &a, const lst &l)
{
- if (is_exactly_a<numeric>(a) || // algorithm does not trap a==0
- is_a<symbol>(a)) // shortcut
+ if (is_exactly_a<numeric>(a) ||
+ is_a<symbol>(a)) // shortcuts
return a;
// If no lst of variables to factorize in was specified we have to
/** Function object to be applied by basic::normal(). */
struct normal_map_function : public map_function {
- int level;
- normal_map_function(int l) : level(l) {}
- ex operator()(const ex & e) override { return normal(e, level); }
+ ex operator()(const ex & e) override { return normal(e); }
};
/** Default implementation of ex::normal(). It normalizes the children and
* replaces the object with a temporary symbol.
* @see ex::normal */
-ex basic::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex basic::normal(exmap & repl, exmap & rev_lookup) const
{
if (nops() == 0)
return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup), _ex1});
- else {
- if (level == 1)
- return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup), _ex1});
- else if (level == -max_recursion_level)
- throw(std::runtime_error("max recursion level reached"));
- else {
- normal_map_function map_normal(level - 1);
- return dynallocate<lst>({replace_with_symbol(map(map_normal), repl, rev_lookup), _ex1});
- }
- }
+
+ normal_map_function map_normal;
+ return dynallocate<lst>({replace_with_symbol(map(map_normal), repl, rev_lookup), _ex1});
}
/** Implementation of ex::normal() for symbols. This returns the unmodified symbol.
* @see ex::normal */
-ex symbol::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex symbol::normal(exmap & repl, exmap & rev_lookup) const
{
return dynallocate<lst>({*this, _ex1});
}
* into re+I*im and replaces I and non-rational real numbers with a temporary
* symbol.
* @see ex::normal */
-ex numeric::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex numeric::normal(exmap & repl, exmap & rev_lookup) const
{
numeric num = numer();
ex numex = num;
/** Implementation of ex::normal() for a sum. It expands terms and performs
* fractional addition.
* @see ex::normal */
-ex add::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex add::normal(exmap & repl, exmap & rev_lookup) const
{
- if (level == 1)
- return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup), _ex1});
- else if (level == -max_recursion_level)
- throw(std::runtime_error("max recursion level reached"));
-
// Normalize children and split each one into numerator and denominator
exvector nums, dens;
nums.reserve(seq.size()+1);
dens.reserve(seq.size()+1);
for (auto & it : seq) {
- ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, level-1);
+ ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup);
nums.push_back(n.op(0));
dens.push_back(n.op(1));
}
- ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, level-1);
+ ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup);
nums.push_back(n.op(0));
dens.push_back(n.op(1));
GINAC_ASSERT(nums.size() == dens.size());
/** Implementation of ex::normal() for a product. It cancels common factors
* from fractions.
* @see ex::normal() */
-ex mul::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex mul::normal(exmap & repl, exmap & rev_lookup) const
{
- if (level == 1)
- return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup), _ex1});
- else if (level == -max_recursion_level)
- throw(std::runtime_error("max recursion level reached"));
-
// Normalize children, separate into numerator and denominator
exvector num; num.reserve(seq.size());
exvector den; den.reserve(seq.size());
ex n;
for (auto & it : seq) {
- n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, level-1);
+ n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup);
num.push_back(n.op(0));
den.push_back(n.op(1));
}
- n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, level-1);
+ n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup);
num.push_back(n.op(0));
den.push_back(n.op(1));
* distributes integer exponents to numerator and denominator, and replaces
* non-integer powers by temporary symbols.
* @see ex::normal */
-ex power::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex power::normal(exmap & repl, exmap & rev_lookup) const
{
- if (level == 1)
- return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup), _ex1});
- else if (level == -max_recursion_level)
- throw(std::runtime_error("max recursion level reached"));
-
// Normalize basis and exponent (exponent gets reassembled)
- ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup, level-1);
- ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, level-1);
+ ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup);
+ ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup);
n_exponent = n_exponent.op(0) / n_exponent.op(1);
if (n_exponent.info(info_flags::integer)) {
/** Implementation of ex::normal() for pseries. It normalizes each coefficient
* and replaces the series by a temporary symbol.
* @see ex::normal */
-ex pseries::normal(exmap & repl, exmap & rev_lookup, int level) const
+ex pseries::normal(exmap & repl, exmap & rev_lookup) const
{
epvector newseq;
for (auto & it : seq) {
* expression can be treated as a rational function). normal() is applied
* recursively to arguments of functions etc.
*
- * @param level maximum depth of recursion
* @return normalized expression */
-ex ex::normal(int level) const
+ex ex::normal() const
{
exmap repl, rev_lookup;
- ex e = bp->normal(repl, rev_lookup, level);
+ ex e = bp->normal(repl, rev_lookup);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
{
exmap repl, rev_lookup;
- ex e = bp->normal(repl, rev_lookup, 0);
+ ex e = bp->normal(repl, rev_lookup);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
{
exmap repl, rev_lookup;
- ex e = bp->normal(repl, rev_lookup, 0);
+ ex e = bp->normal(repl, rev_lookup);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
{
exmap repl, rev_lookup;
- ex e = bp->normal(repl, rev_lookup, 0);
+ ex e = bp->normal(repl, rev_lookup);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
git://www.ginac.de / ginac.git / blobdiff