## "Fossies" - the Fresh Open Source Software Archive

### Source code changes of the file "ginac/normal.cpp" betweenginac-1.7.11.tar.bz2 and ginac-1.8.0.tar.bz2

About: GiNaC (GiNaC is Not a CAS (Computer Algebra System)) is a C++ library for symbolic calculations.

normal.cpp  (ginac-1.7.11.tar.bz2):normal.cpp  (ginac-1.8.0.tar.bz2)
skipping to change at line 1776 skipping to change at line 1776
/** Compute square-free factorization of multivariate polynomial a(x) using /** Compute square-free factorization of multivariate polynomial a(x) using
* Yun's algorithm. Used internally by sqrfree(). * Yun's algorithm. Used internally by sqrfree().
* *
* @param a multivariate polynomial over Z[X], treated here as univariate * @param a multivariate polynomial over Z[X], treated here as univariate
* polynomial in x (needs not be expanded). * polynomial in x (needs not be expanded).
* @param x variable to factor in * @param x variable to factor in
* @return vector of expairs (factor, exponent), sorted by exponents */ * @return vector of expairs (factor, exponent), sorted by exponents */
static epvector sqrfree_yun(const ex &a, const symbol &x) static epvector sqrfree_yun(const ex &a, const symbol &x)
{ {
ex w = a;
ex w =
ex z = w.diff(x); ex z = w.diff(x);
ex g = gcd(w, z); ex g = gcd(w, z);
if (g.is_zero()) { if (g.is_zero()) {
// manifest zero or hidden zero // manifest zero or hidden zero
return {}; return {};
} }
if (g.is_equal(_ex1)) { if (g.is_equal(_ex1)) {
// w(x) and w'(x) share no factors: w(x) is square-free // w(x) and w'(x) share no factors: w(x) is square-free
return {expair(a, _ex1)}; return {expair(a, _ex1)};
} }
epvector factors; epvector factors;
ex i = 0; // exponent ex i = 0; // exponent
do { do {
w = quo(w, g, x); w = quo(w, g, x);
if (w.is_zero()) {
// hidden zero
break;
}
z = quo(z, g, x) - w.diff(x); z = quo(z, g, x) - w.diff(x);
i += 1; i += 1;
if (w.is_equal(x)) { if (w.is_equal(x)) {
// shortcut for x^n with n ∈ ℕ // shortcut for x^n with n ∈ ℕ
i += quo(z, w.diff(x), x); i += quo(z, w.diff(x), x);
factors.push_back(expair(w, i)); factors.push_back(expair(w, i));
break; break;
} }
g = gcd(w, z); g = gcd(w, z);
if (!g.is_equal(_ex1)) { if (!g.is_equal(_ex1)) {
factors.push_back(expair(g, i)); factors.push_back(expair(g, i));
} }
} while (!z.is_zero()); } while (!z.is_zero());
// correct for lost factor // correct for lost factor
// // (being based on GCDs, Yun's algorithm only finds factors up to a unit)
Yun's algorithm only finds
factors up to a
const ex lost_factor = quo(a, mul{factors}, x); const ex lost_factor = quo(a, mul{factors}, x);
if (lost_factor.is_equal(_ex1)) { if (lost_factor.is_equal(_ex1)) {
// trivial lost factor // trivial lost factor
return factors; return factors;
} }
if (factors[0].coeff.is_equal(1)) { if (!factors.empty() && factors[0].coeff.is_equal(1)) {
// multiply factor^1 with lost_factor // multiply factor^1 with lost_factor
factors[0].rest *= lost_factor; factors[0].rest *= lost_factor;
return factors; return factors;
} }
// s // no factor^1: prepend lost_factor^1 to the results
epvector results = {expair(lost_factor, 1)}; epvector results = {expair(lost_factor, 1)};
std::move(factors.begin(), factors.end(), std::back_inserter(results)); std::move(factors.begin(), factors.end(), std::back_inserter(results));
return results; return results;
} }
/** Compute a square-free factorization of a multivariate polynomial in Q[X]. /** Compute a square-free factorization of a multivariate polynomial in Q[X].
* *
* @param a multivariate polynomial over Q[X] (needs not be expanded) * @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 * @param l lst of variables to factor in, may be left empty for autodetection
* @return a square-free factorization of \p a. * @return a square-free factorization of \p a.
skipping to change at line 1994 skipping to change at line 1995
* Note: The internal normal() functions (= basic::normal() and overloaded * Note: The internal normal() functions (= basic::normal() and overloaded
* functions) all return lists of the form {numerator, denominator}. This * functions) all return lists of the form {numerator, denominator}. This
* is to get around mul::eval()'s automatic expansion of numeric coefficients. * is to get around mul::eval()'s automatic expansion of numeric coefficients.
* E.g. (a+b)/3 is automatically converted to a/3+b/3 but we want to keep * E.g. (a+b)/3 is automatically converted to a/3+b/3 but we want to keep
* the information that (a+b) is the numerator and 3 is the denominator. * the information that (a+b) is the numerator and 3 is the denominator.
*/ */
/** Create a symbol for replacing the expression "e" (or return a previously /** Create a symbol for replacing the expression "e" (or return a previously
* assigned symbol). The symbol and expression are appended to repl, for * assigned symbol). The symbol and expression are appended to repl, for
* a later application of subs(). * 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 do
ne 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 expres
sion
* 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 si
gnals
* 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 res
tore
* such iterative substitutions in the right way (this solves 2).
* @see ex::normal */ * @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, ls t & modifier)
{ {
// Since the repl contains replaced expressions we should search for them // Since the repl contains replaced expressions we should search for them
ex e_replaced = e.subs(repl, subs_options::no_pattern); ex e_replaced = e.subs(repl, subs_options::no_pattern);
// Expression already replaced? Then return the assigned symbol // Expression already replaced? Then return the assigned symbol
auto it = rev_lookup.find(e_replaced); auto it = rev_lookup.find(e_replaced);
if (it != rev_lookup.end()) if (it != rev_lookup.end())
return it->second; return it->second;
// The expression can be the base of substituted power, which requires a
more careful search
if (! is_a<numeric>(e_replaced))
for (auto & it : repl)
if (is_a<power>(it.second) && e_replaced.is_equal(it.seco
nd.op(0))) {
ex degree = pow(it.second.op(1), _ex_1);
if (is_a<numeric>(degree) && ex_to<numeric>(degre
e).is_integer())
return pow(it.first, degree);
}
// 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_rep
laced, 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 th
is is simply the power of the existing symbol.
// std::clog << e_replaced << " i
s a " << ratio << " power of " << it.first << std::endl;
return dynallocate<power>(it.firs
t, 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 << " i
s power " << Num << " and "
// << it.first <<
" is power " << denom(ratio) << " of the common base "
// << exp(e_replac
ed.op(0)/Num) << std::endl;
// ... and modify the replacemen
t 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.o
p(0)/Num)});
return dynallocate<power>(es, Num
);
}
}
}
}
} else if (is_a<power>(e_replaced) && !is_a<numeric>(e_replaced.op(0)) //
We do not replace simple monomials like x^3 or sqrt(2)
&& ! (is_a<symbol>(e_replaced.op(0))
&& is_a<numeric>(e_replaced.op(1)) && ex_to<numeric>(e_rep
laced.op(1)).is_integer())) {
for (auto & it : repl) {
if (e_replaced.op(0).is_equal(it.second) // The base is a
n allocated symbol or base of power
|| (is_a<power>(it.second) && e_replaced.op(0).is_equ
al(it.second.op(0)))) {
ex ratio; // We bind together two above cases
if (is_a<power>(it.second))
ratio = normal(e_replaced.op(1) / it.seco
nd.op(1));
else
ratio = e_replaced.op(1);
if (is_a<numeric>(ratio) && ex_to<numeric>(ratio)
.is_rational()) {
// Different powers can be treated as pow
ers of the same basic equation
if (ex_to<numeric>(ratio).is_integer()) {
// If ratio is an integer then th
is 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.firs
t, 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 "
// << pow(e_replac
ed.op(0), e_replaced.op(1)/Num) << std::endl;
// ... and modify the replacemen
t tables
rev_lookup.erase(it.second);
rev_lookup.insert({pow(e_replaced
.op(0), e_replaced.op(1)/Num), es});
repl.erase(it.first);
repl.insert({es, pow(e_replaced.o
p(0), e_replaced.op(1)/Num)});
return dynallocate<power>(es, Num
);
}
}
}
}
// There is no existing substitution, thus we are creating a new
one.
// This needs to be done separately to treat possible occurrences
of
// b = e_replaced.op(0) elsewhere in the expression as pow(b, 1).
ex degree = pow(e_replaced.op(1), _ex_1);
if (is_a<numeric>(degree) && ex_to<numeric>(degree).is_integer())
{
ex es = dynallocate<symbol>();
modifier.append(e_replaced.op(0) == power(es, degree));
repl.insert({es, e_replaced});
rev_lookup.insert({e_replaced, es});
return es;
}
}
// Otherwise create new symbol and add to list, taking care that the // Otherwise create new symbol and add to list, taking care that the
// replacement expression doesn't itself contain symbols from repl, // replacement expression doesn't itself contain symbols from repl,
// because subs() is not recursive // because subs() is not recursive
ex es = dynallocate<symbol>(); ex es = dynallocate<symbol>();
repl.insert(std::make_pair(es, e_replaced)); repl.insert(std::make_pair(es, e_replaced));
rev_lookup.insert(std::make_pair(e_replaced, es)); rev_lookup.insert(std::make_pair(e_replaced, es));
return es; return es;
} }
/** Create a symbol for replacing the expression "e" (or return a previously /** Create a symbol for replacing the expression "e" (or return a previously
skipping to change at line 2045 skipping to change at line 2163
} }
/** Function object to be applied by basic::normal(). */ /** Function object to be applied by basic::normal(). */
struct normal_map_function : public map_function { struct normal_map_function : public map_function {
ex operator()(const ex & e) override { return normal(e); } ex operator()(const ex & e) override { return normal(e); }
}; };
/** Default implementation of ex::normal(). It normalizes the children and /** Default implementation of ex::normal(). It normalizes the children and
* replaces the object with a temporary symbol. * replaces the object with a temporary symbol.
* @see ex::normal */ * @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) if (nops() == 0)
return dynallocate<lst>({replace_with_symbol(*this, repl, rev_loo kup), _ex1}); return dynallocate<lst>({replace_with_symbol(*this, repl, rev_loo kup, modifier), _ex1});
normal_map_function map_normal; normal_map_function map_normal;
dynallocate<lst>({replace_with_symbol(map(map_normal), repl, int nmod = modifier.nops(); // To watch new modifiers to the replacement
_ex1}); 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_patte
rn));
}
return result;
} }
/** Implementation of ex::normal() for symbols. This returns the unmodified symb ol. /** Implementation of ex::normal() for symbols. This returns the unmodified symb ol.
* @see ex::normal */ * @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}); return dynallocate<lst>({*this, _ex1});
} }
/** Implementation of ex::normal() for a numeric. It splits complex numbers /** Implementation of ex::normal() for a numeric. It splits complex numbers
* into re+I*im and replaces I and non-rational real numbers with a temporary * into re+I*im and replaces I and non-rational real numbers with a temporary
* symbol. * symbol.
* @see ex::normal */ * @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(); numeric num = numer();
ex numex = num; ex numex = num;
if (num.is_real()) { if (num.is_real()) {
if (!num.is_integer()) if (!num.is_integer())
numex = replace_with_symbol(numex, repl, rev_lookup); numex = replace_with_symbol(numex, repl, rev_lookup, modi fier);
} else { // complex } else { // complex
numeric re = num.real(), im = num.imag(); numeric re = num.real(), im = num.imag();
ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl, 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, 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, numex = re_ex + im_ex * replace_with_symbol(I, repl, rev_lookup,
modifier);
} }
// Denominator is always a real integer (see numeric::denom()) // Denominator is always a real integer (see numeric::denom())
return dynallocate<lst>({numex, denom()}); return dynallocate<lst>({numex, denom()});
} }
/** Fraction cancellation. /** Fraction cancellation.
* @param n numerator * @param n numerator
* @param d denominator * @param d denominator
* @return cancelled fraction {n, d} as a list */ * @return cancelled fraction {n, d} as a list */
skipping to change at line 2147 skipping to change at line 2273
} }
// Return result as list // Return result as list
//std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << std::endl; //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << std::endl;
return dynallocate<lst>({num * pre_factor.numer(), den * pre_factor.denom ()}); return dynallocate<lst>({num * pre_factor.numer(), den * pre_factor.denom ()});
} }
/** Implementation of ex::normal() for a sum. It expands terms and performs /** Implementation of ex::normal() for a sum. It expands terms and performs
* @see ex::normal */ * @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 // Normalize children and split each one into numerator and denominator
exvector nums, dens; exvector nums, dens;
nums.reserve(seq.size()+1); nums.reserve(seq.size()+1);
dens.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) { for (auto & it : seq) {
ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lo okup); ex n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_lo okup, modifier);
nums.push_back(n.op(0)); nums.push_back(n.op(0));
dens.push_back(n.op(1)); 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)); nums.push_back(n.op(0));
dens.push_back(n.op(1)); dens.push_back(n.op(1));
GINAC_ASSERT(nums.size() == dens.size()); GINAC_ASSERT(nums.size() == dens.size());
// Now, nums is a vector of all numerators and dens is a vector of // Now, nums is a vector of all numerators and dens is a vector of
// all denominators // all denominators
//std::clog << "add::normal uses " << nums.size() << " summands:\n"; //std::clog << "add::normal uses " << nums.size() << " summands:\n";
auto num_it = nums.begin(), num_itend = nums.end(); auto num_it = nums.begin(), num_itend = nums.end();
auto den_it = dens.begin(), den_itend = dens.end(); auto den_it = dens.begin(), den_itend = dens.end();
//std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl; //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::n
o_pattern);
++num_it;
*den_it = den_it->subs(modifier.op(imod), subs_options::n
o_pattern);
++den_it;
}
// Reset iterators for the next round
num_it = nums.begin();
den_it = dens.begin();
}
ex num = *num_it++, den = *den_it++; ex num = *num_it++, den = *den_it++;
while (num_it != num_itend) { while (num_it != num_itend) {
//std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl; //std::clog << " num = " << *num_it << ", den = " << *den_it << std::endl;
ex next_num = *num_it++, next_den = *den_it++; ex next_num = *num_it++, next_den = *den_it++;
// Trivially add sequences of fractions with identical denominato rs // Trivially add sequences of fractions with identical denominato rs
while ((den_it != den_itend) && next_den.is_equal(*den_it)) { while ((den_it != den_itend) && next_den.is_equal(*den_it)) {
next_num += *num_it; next_num += *num_it;
num_it++; den_it++; num_it++; den_it++;
} }
skipping to change at line 2198 skipping to change at line 2337
} }
//std::clog << " common denominator = " << den << std::endl; //std::clog << " common denominator = " << den << std::endl;
// Cancel common factors from num/den // Cancel common factors from num/den
return frac_cancel(num, den); return frac_cancel(num, den);
} }
/** Implementation of ex::normal() for a product. It cancels common factors /** Implementation of ex::normal() for a product. It cancels common factors
* from fractions. * from fractions.
* @see ex::normal() */ * @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 // Normalize children, separate into numerator and denominator
exvector num; num.reserve(seq.size()); exvector num; num.reserve(seq.size());
exvector den; den.reserve(seq.size()); exvector den; den.reserve(seq.size());
ex n; ex n;
int nmod = modifier.nops(); // To watch new modifiers to the replacement list
for (auto & it : seq) { for (auto & it : seq) {
n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_looku p); n = ex_to<basic>(recombine_pair_to_ex(it)).normal(repl, rev_looku p, modifier);
num.push_back(n.op(0)); num.push_back(n.op(0));
den.push_back(n.op(1)); 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)); num.push_back(n.op(0));
den.push_back(n.op(1)); 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::n
o_pattern);
++num_it;
*den_it = den_it->subs(modifier.op(imod), subs_options::n
o_pattern);
++den_it;
}
num_it = num.begin();
den_it = den.begin();
}
// Perform fraction cancellation // Perform fraction cancellation
return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den)); return frac_cancel(dynallocate<mul>(num), dynallocate<mul>(den));
} }
/** Implementation of ex::normal() 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 * distributes integer exponents to numerator and denominator, and replaces
* non-integer powers by temporary symbols. * non-integer powers by temporary symbols.
* @see ex::normal */ * @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) // Normalize basis and exponent (exponent gets reassembled)
ex n_basis = ex_to<basic>(basis).normal(repl, int nmod = modifier.nops(); // To watch new modifiers to the replacement
ex n_exponent = ex_to<basic>(exponent).normal(repl, list
ex n_basis = ex_to<basic>(basis).normal(repl, rev_lookup, modifier);
for (int imod = nmod; imod < modifier.nops(); ++imod)
n_basis = n_basis.subs(modifier.op(imod), subs_options::no_patter
n);
nmod = modifier.nops();
ex n_exponent = ex_to<basic>(exponent).normal(repl, rev_lookup, modifier)
;
for (int imod = nmod; imod < modifier.nops(); ++imod)
n_exponent = n_exponent.subs(modifier.op(imod), subs_options::no_
pattern);
n_exponent = n_exponent.op(0) / n_exponent.op(1); n_exponent = n_exponent.op(0) / n_exponent.op(1);
if (n_exponent.info(info_flags::integer)) { if (n_exponent.info(info_flags::integer)) {
if (n_exponent.info(info_flags::positive)) { if (n_exponent.info(info_flags::positive)) {
// (a/b)^n -> {a^n, b^n} // (a/b)^n -> {a^n, b^n}
return dynallocate<lst>({pow(n_basis.op(0), n_exponent), pow(n_basis.op(1), n_exponent)}); return dynallocate<lst>({pow(n_basis.op(0), n_exponent), pow(n_basis.op(1), n_exponent)});
} else if (n_exponent.info(info_flags::negative)) { } else if (n_exponent.info(info_flags::negative)) {
// (a/b)^-n -> {b^n, a^n} // (a/b)^-n -> {b^n, a^n}
return dynallocate<lst>({pow(n_basis.op(1), -n_exponent), pow(n_basis.op(0), -n_exponent)}); return dynallocate<lst>({pow(n_basis.op(1), -n_exponent), pow(n_basis.op(0), -n_exponent)});
} }
} else { } else {
if (n_exponent.info(info_flags::positive)) { if (n_exponent.info(info_flags::positive)) {
// (a/b)^x -> {sym((a/b)^x), 1} // (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)) { } else if (n_exponent.info(info_flags::negative)) {
if (n_basis.op(1).is_equal(_ex1)) { if (n_basis.op(1).is_equal(_ex1)) {
// a^-x -> {1, sym(a^x)} // a^-x -> {1, sym(a^x)}
return dynallocate<lst>({_ex1, replace_with_symbo l(pow(n_basis.op(0), -n_exponent), repl, rev_lookup)}); return dynallocate<lst>({_ex1, replace_with_symbo l(pow(n_basis.op(0), -n_exponent), repl, rev_lookup, modifier)});
} else { } else {
// (a/b)^-x -> {sym((b/a)^x), 1} // (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; 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} // (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 /** Implementation of ex::normal() for pseries. It normalizes each coefficient
* and replaces the series by a temporary symbol. * and replaces the series by a temporary symbol.
* @see ex::normal */ * @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; epvector newseq;
for (auto & it : seq) { for (auto & it : seq) {
ex restexp = it.rest.normal(); ex restexp = it.rest.normal();
if (!restexp.is_zero()) if (!restexp.is_zero())
newseq.push_back(expair(restexp, it.coeff)); newseq.push_back(expair(restexp, it.coeff));
} }
ex n = pseries(relational(var,point), std::move(newseq)); ex n = pseries(relational(var,point), std::move(newseq));
return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup; return dynallocate<lst>({replace_with_symbol(n, repl, rev_lookup, modifie r), _ex1});
} }
/** Normalization of rational functions. /** Normalization of rational functions.
* This function converts an expression to its normal form * This function converts an expression to its normal form
* "numerator/denominator", where numerator and denominator are (relatively * "numerator/denominator", where numerator and denominator are (relatively
* prime) polynomials. Any subexpressions which are not rational functions * prime) polynomials. Any subexpressions which are not rational functions
* (like non-rational numbers, non-integer powers or functions like sin(), * (like non-rational numbers, non-integer powers or functions like sin(),
* cos() etc.) are replaced by temporary symbols which are re-substituted by * cos() etc.) are replaced by temporary symbols which are re-substituted by
* the (normalized) subexpressions before normal() returns (this way, any * the (normalized) subexpressions before normal() returns (this way, any
* expression can be treated as a rational function). normal() is applied * expression can be treated as a rational function). normal() is applied
* recursively to arguments of functions etc. * recursively to arguments of functions etc.
* *
* @return normalized expression */ * @return normalized expression */
ex ex::normal() const ex ex::normal() const
{ {
exmap repl, rev_lookup; 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)); GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols // 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); e = e.subs(repl, subs_options::no_pattern);
}
// Convert {numerator, denominator} form back to fraction // Convert {numerator, denominator} form back to fraction
return e.op(0) / e.op(1); return e.op(0) / e.op(1);
} }
/** Get numerator of an expression. If the expression is not of the normal /** Get numerator of an expression. If the expression is not of the normal
* form "numerator/denominator", it is first converted to this form and * form "numerator/denominator", it is first converted to this form and
* then the numerator is returned. * then the numerator is returned.
* *
* @see ex::normal * @see ex::normal
* @return numerator */ * @return numerator */
ex ex::numer() const ex ex::numer() const
{ {
exmap repl, rev_lookup; 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)); GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols // Re-insert replaced symbols
if (repl.empty()) if (repl.empty())
return e.op(0); 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); return e.op(0).subs(repl, subs_options::no_pattern);
}
} }
/** Get denominator of an expression. If the expression is not of the normal /** Get denominator of an expression. If the expression is not of the normal
* form "numerator/denominator", it is first converted to this form and * form "numerator/denominator", it is first converted to this form and
* then the denominator is returned. * then the denominator is returned.
* *
* @see ex::normal * @see ex::normal
* @return denominator */ * @return denominator */
ex ex::denom() const ex ex::denom() const
{ {
exmap repl, rev_lookup; 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)); GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols // Re-insert replaced symbols
if (repl.empty()) if (repl.empty())
return e.op(1); 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); return e.op(1).subs(repl, subs_options::no_pattern);
}
} }
/** Get numerator and denominator of an expression. If the expression is not /** Get numerator and denominator of an expression. If the expression is not
* of the normal form "numerator/denominator", it is first converted to this * of the normal form "numerator/denominator", it is first converted to this
* form and then a list [numerator, denominator] is returned. * form and then a list [numerator, denominator] is returned.
* *
* @see ex::normal * @see ex::normal
* @return a list [numerator, denominator] */ * @return a list [numerator, denominator] */
ex ex::numer_denom() const ex ex::numer_denom() const
{ {
exmap repl, rev_lookup; 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)); GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols // Re-insert replaced symbols
if (repl.empty()) if (repl.empty())
return e; 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 e.subs(repl, subs_options::no_pattern);
}
} }
/** Rationalization of non-rational functions. /** Rationalization of non-rational functions.
* This function converts a general expression to a rational function * This function converts a general expression to a rational function
* by replacing all non-rational subexpressions (like non-rational numbers, * by replacing all non-rational subexpressions (like non-rational numbers,
* non-integer powers or functions like sin(), cos() etc.) to temporary * non-integer powers or functions like sin(), cos() etc.) to temporary
* symbols. This makes it possible to use functions like gcd() and divide() * symbols. This makes it possible to use functions like gcd() and divide()
* on non-rational functions by applying to_rational() on the arguments, * on non-rational functions by applying to_rational() on the arguments,
* calling the desired function and re-substituting the temporary symbols * calling the desired function and re-substituting the temporary symbols
* in the result. To make the last step possible, all temporary symbols and * in the result. To make the last step possible, all temporary symbols and
End of changes. 50 change blocks.
44 lines changed or deleted 282 lines changed or added