* computation, square-free factorization and rational function normalization. */
/*
- * GiNaC Copyright (C) 1999-2016 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
/** Maximum number of terms of leading coefficient of symbol in both polynomials */
size_t max_lcnops;
- /** Commparison operator for sorting */
+ /** Comparison operator for sorting */
bool operator<(const sym_desc &x) const
{
if (max_deg == x.max_deg)
#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
}
// Some trivial cases
- ex aex = a.expand(), bex = b.expand();
+ ex aex = a.expand();
if (aex.is_zero()) {
if (ca)
*ca = _ex0;
*cb = _ex1;
return b;
}
+ ex bex = b.expand();
if (bex.is_zero()) {
if (ca)
*ca = _ex1;
// The symbol with least degree which is contained in both polynomials
// is our main variable
- sym_desc_vec::iterator vari = sym_stats.begin();
+ auto vari = sym_stats.begin();
while ((vari != sym_stats.end()) &&
(((vari->ldeg_b == 0) && (vari->deg_b == 0)) ||
((vari->ldeg_a == 0) && (vari->deg_a == 0))))
*cb = b;
return _ex1;
}
- // move symbols which 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:
if (p.is_equal(b)) {
// a = p^n, b = p, gcd = p
if (ca)
- *ca = pow(p, a.op(1) - 1);
+ *ca = pow(p, exp_a - 1);
if (cb)
*cb = _ex1;
return p;
- }
+ }
+ if (is_a<symbol>(p)) {
+ // Cancel trivial common factor
+ int ldeg_a = ex_to<numeric>(exp_a).to_int();
+ int ldeg_b = b.ldegree(p);
+ int min_ldeg = std::min(ldeg_a, ldeg_b);
+ if (min_ldeg > 0) {
+ ex common = pow(p, min_ldeg);
+ return gcd(pow(p, ldeg_a - min_ldeg), (b / common).expand(), ca, cb, false) * common;
+ }
+ }
ex p_co, bpart_co;
ex p_gcd = gcd(p, b, &p_co, &bpart_co, false);
- // a(x) = p(x)^n, gcd(p, b) = 1 ==> gcd(a, b) = 1
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)
* @param a multivariate polynomial over Z[X], treated here as univariate
* 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)
+ * @return vector of expairs (factor, exponent), sorted by exponents */
+static epvector sqrfree_yun(const ex &a, const symbol &x)
{
- exvector res;
ex w = a;
ex z = w.diff(x);
ex g = gcd(w, z);
if (g.is_zero()) {
- return res;
+ // manifest zero or hidden zero
+ return {};
}
if (g.is_equal(_ex1)) {
- res.push_back(a);
- return res;
+ // w(x) and w'(x) share no factors: w(x) is square-free
+ return {expair(a, _ex1)};
}
- ex y;
+
+ epvector factors;
+ ex i = 0; // exponent
do {
w = quo(w, g, x);
if (w.is_zero()) {
- return res;
+ // hidden zero
+ break;
+ }
+ z = quo(z, g, x) - w.diff(x);
+ i += 1;
+ if (w.is_equal(x)) {
+ // shortcut for x^n with n ∈ ℕ
+ i += quo(z, w.diff(x), x);
+ factors.push_back(expair(w, i));
+ break;
}
- y = quo(z, g, x);
- z = y - w.diff(x);
g = gcd(w, z);
- res.push_back(g);
+ if (!g.is_equal(_ex1)) {
+ factors.push_back(expair(g, i));
+ }
} while (!z.is_zero());
- return res;
+
+ // 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
- exvector factors = sqrfree_yun(tmp, x);
+ epvector factors = sqrfree_yun(tmp, x);
+ if (factors.empty()) {
+ // the polynomial was a hidden zero
+ return _ex0;
+ }
- // construct the next list of symbols with the first element popped
- lst newargs = args;
- newargs.remove_first();
+ // remove symbol x and proceed recursively with the remaining symbols
+ args.remove_first();
// recurse down the factors in remaining variables
- if (newargs.nops()>0) {
+ if (args.nops()>0) {
for (auto & it : factors)
- it = sqrfree(it, newargs);
+ it.rest = sqrfree(it.rest, args);
}
// Done with recursion, now construct the final result
- ex result = _ex1;
- int p = 1;
- for (auto & it : factors)
- result *= pow(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. 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);
+ 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
- exvector yun = sqrfree_yun(denom, x);
-//clog << "yun factors: " << exprseq(yun) << endl;
- size_t num_yun = yun.size();
- exvector factor; factor.reserve(num_yun);
- exvector cofac; cofac.reserve(num_yun);
- for (size_t i=0; i<num_yun; i++) {
- if (!yun[i].is_equal(_ex1)) {
- for (size_t j=0; j<=i; j++) {
- factor.push_back(pow(yun[i], j+1));
- ex prod = _ex1;
- for (size_t k=0; k<num_yun; k++) {
- if (k == i)
- prod *= pow(yun[k], i-j);
- else
- prod *= pow(yun[k], k+1);
- }
- cofac.push_back(prod.expand());
+ epvector yun = sqrfree_yun(denom, x);
+ size_t yun_max_exponent = yun.empty() ? 0 : ex_to<numeric>(yun.back().coeff).to_int();
+ exvector factor, cofac;
+ for (size_t i=0; i<yun.size(); i++) {
+ numeric i_exponent = ex_to<numeric>(yun[i].coeff);
+ for (size_t j=0; j<i_exponent; j++) {
+ factor.push_back(pow(yun[i].rest, j+1));
+ ex prod = _ex1;
+ for (size_t k=0; k<yun.size(); k++) {
+ if (yun[k].coeff == i_exponent)
+ prod *= pow(yun[k].rest, i_exponent-1-j);
+ else
+ prod *= pow(yun[k].rest, yun[k].coeff);
}
+ cofac.push_back(prod.expand());
}
}
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);
/** Create a symbol for replacing the expression "e" (or return a previously
* assigned symbol). The symbol and expression are appended to repl, for
* a later application of subs().
+ * An entry in the replacement table repl can be changed in some cases.
+ * If it was altered, we need to provide the modifier for the previously build expressions.
+ * The modifier is an (ordered) list, because those substitutions need to be done in the
+ * incremental order.
+ * As an example let us consider a rationalisation of the expression
+ * e = exp(2*x)*cos(exp(2*x)+1)*exp(x)
+ * The first factor GiNaC denotes by something like symbol1 and will record:
+ * e =symbol1*cos(symbol1 + 1)*exp(x)
+ * repl = {symbol1 : exp(2*x)}
+ * Similarly, the second factor would be denoted as symbol2 and we will have
+ * e =symbol1*symbol2*exp(x)
+ * repl = {symbol1 : exp(2*x), symbol2 : cos(symbol1 + 1)}
+ * Denoting the third term as symbol3 GiNaC is willing to re-think exp(2*x) as
+ * symbol3^2 rather than just symbol1. Here are two issues:
+ * 1) The replacement "symbol1 -> symbol3^2" in the previous part of the expression
+ * needs to be done outside of the present routine;
+ * 2) The pair "symbol1 : exp(2*x)" shall be deleted from the replacement table repl.
+ * However, this will create illegal substitution "symbol2 : cos(symbol1 + 1)" with
+ * undefined symbol1.
+ * These both problems are mitigated through the additions of the record
+ * "symbol1==symbol3^2" to the list modifier. Changed length of the modifier signals
+ * to the calling code that the previous portion of the expression needs to be
+ * altered (it solves 1). Thus GiNaC can record now
+ * e =symbol3^2*symbol2*symbol3
+ * repl = {symbol2 : cos(symbol1 + 1), symbol3 : exp(x)}
+ * modifier = {symbol1==symbol3^2}
+ * Then, doing the backward substitutions the list modifier will be used to restore
+ * such iterative substitutions in the right way (this solves 2).
* @see ex::normal */
-static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup)
+static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup, lst & modifier)
{
// Since the repl contains replaced expressions we should search for them
ex e_replaced = e.subs(repl, subs_options::no_pattern);
if (it != rev_lookup.end())
return it->second;
+ // We treat powers and the exponent functions differently because
+ // they can be rationalised more efficiently
+ if (is_a<function>(e_replaced) && is_ex_the_function(e_replaced, exp)) {
+ for (auto & it : repl) {
+ if (is_a<function>(it.second) && is_ex_the_function(e_replaced, exp)) {
+ ex ratio = normal(e_replaced.op(0) / it.second.op(0));
+ if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).is_rational()) {
+ // Different exponents can be treated as powers of the same basic equation
+ if (ex_to<numeric>(ratio).is_integer()) {
+ // If ratio is an integer then this is simply the power of the existing symbol.
+ // std::clog << e_replaced << " is a " << ratio << " power of " << it.first << std::endl;
+ return dynallocate<power>(it.first, ratio);
+ } else {
+ // otherwise we need to give the replacement pattern to change
+ // the previous expression...
+ ex es = dynallocate<symbol>();
+ ex Num = numer(ratio);
+ modifier.append(it.first == power(es, denom(ratio)));
+ // std::clog << e_replaced << " is power " << Num << " and "
+ // << it.first << " is power " << denom(ratio) << " of the common base "
+ // << exp(e_replaced.op(0)/Num) << std::endl;
+ // ... and modify the replacement tables
+ rev_lookup.erase(it.second);
+ rev_lookup.insert({exp(e_replaced.op(0)/Num), es});
+ repl.erase(it.first);
+ repl.insert({es, exp(e_replaced.op(0)/Num)});
+ return dynallocate<power>(es, Num);
+ }
+ }
+ }
+ }
+ }
+
// Otherwise create new symbol and add to list, taking care that the
// replacement expression doesn't itself contain symbols from repl,
// because subs() is not recursive
/** 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) const
+ex basic::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
{
if (nops() == 0)
- return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup), _ex1});
+ return dynallocate<lst>({replace_with_symbol(*this, repl, rev_lookup, modifier), _ex1});
normal_map_function map_normal;
- return dynallocate<lst>({replace_with_symbol(map(map_normal), repl, rev_lookup), _ex1});
+ int nmod = modifier.nops(); // To watch new modifiers to the replacement list
+ lst result = dynallocate<lst>({replace_with_symbol(map(map_normal), repl, rev_lookup, modifier), _ex1});
+ for (int imod = nmod; imod < modifier.nops(); ++imod) {
+ exmap this_repl;
+ this_repl.insert(std::make_pair(modifier.op(imod).op(0), modifier.op(imod).op(1)));
+ result = ex_to<lst>(result.subs(this_repl, subs_options::no_pattern));
+ }
+
+ return result;
}
/** Implementation of ex::normal() for symbols. This returns the unmodified symbol.
* @see ex::normal */
-ex symbol::normal(exmap & repl, exmap & rev_lookup) const
+ex symbol::normal(exmap & repl, exmap & rev_lookup, lst & modifier) 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) const
+ex numeric::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
{
numeric num = numer();
ex numex = num;
if (num.is_real()) {
if (!num.is_integer())
- numex = replace_with_symbol(numex, repl, rev_lookup);
+ numex = replace_with_symbol(numex, repl, rev_lookup, modifier);
} else { // complex
numeric re = num.real(), im = num.imag();
- ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl, rev_lookup);
- ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl, rev_lookup);
- numex = re_ex + im_ex * replace_with_symbol(I, repl, rev_lookup);
+ ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl, rev_lookup, modifier);
+ ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl, rev_lookup, modifier);
+ numex = re_ex + im_ex * replace_with_symbol(I, repl, rev_lookup, modifier);
}
// Denominator is always a real integer (see numeric::denom())
/** 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) const
+ex add::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
{
// Normalize children and split each one into numerator and denominator
exvector nums, dens;
nums.reserve(seq.size()+1);
dens.reserve(seq.size()+1);
+ int nmod = modifier.nops(); // To watch new modifiers to the replacement list
for (auto & it : seq) {
- ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup);
+ ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
nums.push_back(n.op(0));
dens.push_back(n.op(1));
}
- ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup);
+ ex n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier);
nums.push_back(n.op(0));
dens.push_back(n.op(1));
GINAC_ASSERT(nums.size() == dens.size());
auto num_it = nums.begin(), num_itend = nums.end();
auto den_it = dens.begin(), den_itend = dens.end();
//std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
+ for (int imod = nmod; imod < modifier.nops(); ++imod) {
+ while (num_it != num_itend) {
+ *num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
+ ++num_it;
+ *den_it = den_it->subs(modifier.op(imod), subs_options::no_pattern);
+ ++den_it;
+ }
+ // Reset iterators for the next round
+ num_it = nums.begin();
+ den_it = dens.begin();
+ }
+
ex num = *num_it++, den = *den_it++;
while (num_it != num_itend) {
//std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
/** Implementation of ex::normal() for a product. It cancels common factors
* from fractions.
* @see ex::normal() */
-ex mul::normal(exmap & repl, exmap & rev_lookup) const
+ex mul::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
{
// Normalize children, separate into numerator and denominator
exvector num; num.reserve(seq.size());
exvector den; den.reserve(seq.size());
ex n;
+ int nmod = modifier.nops(); // To watch new modifiers to the replacement list
for (auto & it : seq) {
- n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup);
+ n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lookup, modifier);
num.push_back(n.op(0));
den.push_back(n.op(1));
}
- n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup);
+ n = ex_to<numeric>(overall_coeff).normal(repl, rev_lookup, modifier);
num.push_back(n.op(0));
den.push_back(n.op(1));
+ auto num_it = num.begin(), num_itend = num.end();
+ auto den_it = den.begin(), den_itend = den.end();
+ for (int imod = nmod; imod < modifier.nops(); ++imod) {
+ while (num_it != num_itend) {
+ *num_it = num_it->subs(modifier.op(imod), subs_options::no_pattern);
+ ++num_it;
+ *den_it = den_it->subs(modifier.op(imod), subs_options::no_pattern);
+ ++den_it;
+ }
+ num_it = num.begin();
+ den_it = den.begin();
+ }
// Perform fraction cancellation
return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den));
}
-/** Implementation of ex::normal([B) for powers. It normalizes the basis,
+/** Implementation of ex::normal() for powers. It normalizes the basis,
* 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) const
+ex power::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
{
// Normalize basis and exponent (exponent gets reassembled)
- ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup);
- ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup);
+ ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup, modifier);
+ ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, modifier);
n_exponent = n_exponent.op(0) / n_exponent.op(1);
if (n_exponent.info(info_flags::integer)) {
if (n_exponent.info(info_flags::positive)) {
// (a/b)^x -> {sym((a/b)^x), 1}
- return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup), _ex1});
+ return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup, modifier), _ex1});
} else if (n_exponent.info(info_flags::negative)) {
if (n_basis.op(1).is_equal(_ex1)) {
// a^-x -> {1, sym(a^x)}
- return dynallocate<lst>({_ex1, replace_with_symbol(pow(n_basis.op(0), -n_exponent), repl, rev_lookup)});
+ return dynallocate<lst>({_ex1, replace_with_symbol(pow(n_basis.op(0), -n_exponent), repl, rev_lookup, modifier)});
} else {
// (a/b)^-x -> {sym((b/a)^x), 1}
- return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup), _ex1});
+ return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup, modifier), _ex1});
}
}
}
// (a/b)^x -> {sym((a/b)^x, 1}
- return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup), _ex1});
+ return dynallocate<lst>({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup, modifier), _ex1});
}
/** 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) const
+ex pseries::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const
{
epvector newseq;
for (auto & it : seq) {
newseq.push_back(expair(restexp, it.coeff));
}
ex n = pseries(relational(var,point), std::move(newseq));
- return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup), _ex1});
+ return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup, modifier), _ex1});
}
ex ex::normal() const
{
exmap repl, rev_lookup;
+ lst modifier;
- ex e = bp->normal(repl, rev_lookup);
+ ex e = bp->normal(repl, rev_lookup, modifier);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
- if (!repl.empty())
+ if (!repl.empty()) {
+ for(int i=0; i < modifier.nops(); ++i)
+ e = e.subs(modifier.op(i), subs_options::no_pattern);
e = e.subs(repl, subs_options::no_pattern);
+ }
// Convert {numerator, denominator} form back to fraction
return e.op(0) / e.op(1);
ex ex::numer() const
{
exmap repl, rev_lookup;
+ lst modifier;
- ex e = bp->normal(repl, rev_lookup);
+ ex e = bp->normal(repl, rev_lookup, modifier);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
if (repl.empty())
return e.op(0);
- else
+ else {
+ for(int i=0; i < modifier.nops(); ++i)
+ e = e.subs(modifier.op(i), subs_options::no_pattern);
+
return e.op(0).subs(repl, subs_options::no_pattern);
+ }
}
/** Get denominator of an expression. If the expression is not of the normal
ex ex::denom() const
{
exmap repl, rev_lookup;
+ lst modifier;
- ex e = bp->normal(repl, rev_lookup);
+ ex e = bp->normal(repl, rev_lookup, modifier);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
if (repl.empty())
return e.op(1);
- else
+ else {
+ for(int i=0; i < modifier.nops(); ++i)
+ e = e.subs(modifier.op(i), subs_options::no_pattern);
+
return e.op(1).subs(repl, subs_options::no_pattern);
+ }
}
/** Get numerator and denominator of an expression. If the expression is not
ex ex::numer_denom() const
{
exmap repl, rev_lookup;
+ lst modifier;
- ex e = bp->normal(repl, rev_lookup);
+ ex e = bp->normal(repl, rev_lookup, modifier);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
if (repl.empty())
return e;
- else
+ else {
+ for(int i=0; i < modifier.nops(); ++i)
+ e = e.subs(modifier.op(i), subs_options::no_pattern);
+
return e.subs(repl, subs_options::no_pattern);
+ }
}
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 (auto & it : repl_lst)
- 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 (auto & it : m)
- repl_lst.append(it.first == it.second);
-
- return ret;
-}
-
ex ex::to_polynomial(exmap & repl) const
{
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 (auto & it : repl_lst)
- 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 (auto & it : m)
- 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(exmap & repl) const
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