ex eval(int level=0) const;
ex evalm() const;
ex series(const relational & r, int order, unsigned options = 0) const;
- ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
+ ex normal(exmap & repl, int level=0) const;
numeric integer_content() const;
ex smod(const numeric &xi) const;
numeric max_coefficient() const;
}
/** Helper function for subs(). Does not recurse into subexpressions. */
-ex basic::subs_one_level(const lst & ls, const lst & lr, unsigned options) const
+ex basic::subs_one_level(const exmap & m, unsigned options) const
{
- GINAC_ASSERT(ls.nops() == lr.nops());
-
- lst::const_iterator its, itr;
+ exmap::const_iterator it;
if (options & subs_options::subs_no_pattern) {
- for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
- if (is_equal(ex_to<basic>(*its)))
- return *itr;
- }
+ it = m.find(*this);
+ if (it != m.end())
+ return it->second;
} else {
- for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+ for (it = m.begin(); it != m.end(); ++it) {
lst repl_lst;
- if (match(ex_to<basic>(*its), repl_lst))
- return itr->subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards
+ if (match(ex_to<basic>(it->first), repl_lst))
+ return it->second.subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards
}
}
/** Substitute a set of objects by arbitrary expressions. The ex returned
* will already be evaluated. */
-ex basic::subs(const lst & ls, const lst & lr, unsigned options) const
+ex basic::subs(const exmap & m, unsigned options) const
{
size_t num = nops();
if (num) {
// Substitute in subexpressions
for (size_t i=0; i<num; i++) {
const ex & orig_op = op(i);
- const ex & subsed_op = orig_op.subs(ls, lr, options);
+ const ex & subsed_op = orig_op.subs(m, options);
if (!are_ex_trivially_equal(orig_op, subsed_op)) {
// Something changed, clone the object
// Substitute the other operands
for (; i<num; i++)
- copy->let_op(i) = op(i).subs(ls, lr, options);
+ copy->let_op(i) = op(i).subs(m, options);
// Perform substitutions on the new object as a whole
- return copy->subs_one_level(ls, lr, options);
+ return copy->subs_one_level(m, options);
}
}
}
// Nothing changed or no subexpressions
- return subs_one_level(ls, lr, options);
+ return subs_one_level(m, options);
}
/** Default interface of nth derivative ex::diff(s, n). It should be called
// public
-/** Substitute objects in an expression (syntactic substitution) and return
- * the result as a new expression. There are two valid types of
- * replacement arguments: 1) a relational like object==ex and 2) a list of
- * relationals lst(object1==ex1,object2==ex2,...), which is converted to
- * subs(lst(object1,object2,...),lst(ex1,ex2,...)). */
-ex basic::subs(const ex & e, unsigned options) const
-{
- if (e.info(info_flags::relation_equal)) {
- return subs(lst(e.lhs()), lst(e.rhs()), options);
- }
- if (!e.info(info_flags::list)) {
- throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
- }
-
- // Split list into two
- lst ls;
- lst lr;
- GINAC_ASSERT(is_a<lst>(e));
- for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
- ex r = *it;
- if (!r.info(info_flags::relation_equal)) {
- throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
- }
- ls.append(r.op(0));
- lr.append(r.op(1));
- }
- return subs(ls, lr, options);
-}
-
/** Compare objects syntactically to establish canonical ordering.
* All compare functions return: -1 for *this less than other, 0 equal,
* 1 greater. */
#include <cstddef> // for size_t
#include <vector>
+#include <map>
// CINT needs <algorithm> to work properly with <vector>
#include <algorithm>
template <class> class ptr;
typedef std::vector<ex> exvector;
+typedef std::map<ex, ex, ex_is_less> exmap;
/** Function object for map(). */
public:
// substitutions
- virtual ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+ virtual ex subs(const exmap & m, unsigned options = 0) const;
// function mapping
virtual ex map(map_function & f) const;
virtual ex series(const relational & r, int order, unsigned options = 0) const;
// rational functions
- virtual ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+ virtual ex normal(exmap & repl, int level = 0) const;
virtual ex to_rational(lst &repl_lst) const;
virtual ex to_polynomial(lst &repl_lst) const;
// non-virtual functions in this class
public:
- ex subs(const ex & e, unsigned options = 0) const;
- ex subs_one_level(const lst & ls, const lst & lr, unsigned options) const;
+ ex subs_one_level(const exmap & m, unsigned options) const;
ex diff(const symbol & s, unsigned nth = 1) const;
int compare(const basic & other) const;
bool is_equal(const basic & other) const;
ex op(size_t i) const;
ex & let_op(size_t i);
ex eval(int level = 0) const;
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+ ex subs(const exmap & m, unsigned options = 0) const;
protected:
bool is_equal_same_type(const basic & other) const;
protected:
STLT evalchildren(int level) const;
- STLT *subschildren(const lst & ls, const lst & lr, unsigned options = 0) const;
+ STLT *subschildren(const exmap & m, unsigned options = 0) const;
};
/** Default constructor */
}
template <template <class> class C>
-ex container<C>::subs(const lst & ls, const lst & lr, unsigned options) const
+ex container<C>::subs(const exmap & m, unsigned options) const
{
- STLT *vp = subschildren(ls, lr, options);
+ STLT *vp = subschildren(m, options);
if (vp)
- return ex_to<basic>(thiscontainer(vp)).subs_one_level(ls, lr, options);
+ return ex_to<basic>(thiscontainer(vp)).subs_one_level(m, options);
else
- return subs_one_level(ls, lr, options);
+ return subs_one_level(m, options);
}
/** Compare two containers of the same type. */
}
template <template <class> class C>
-typename container<C>::STLT *container<C>::subschildren(const lst & ls, const lst & lr, unsigned options) const
+typename container<C>::STLT *container<C>::subschildren(const exmap & m, unsigned options) const
{
// returns a NULL pointer if nothing had to be substituted
// returns a pointer to a newly created epvector otherwise
const_iterator cit = seq.begin(), end = seq.end();
while (cit != end) {
- const ex & subsed_ex = cit->subs(ls, lr, options);
+ const ex & subsed_ex = cit->subs(m, options);
if (!are_ex_trivially_equal(*cit, subsed_ex)) {
// copy first part of seq which hasn't changed
// copy rest
while (cit != end) {
- s->push_back(cit->subs(ls, lr, options));
+ s->push_back(cit->subs(m, options));
++cit;
}
return any_found;
}
+/** Substitute objects in an expression (syntactic substitution) and return
+ * the result as a new expression. */
+ex ex::subs(const lst & ls, const lst & lr, unsigned options) const
+{
+ GINAC_ASSERT(ls.nops() == lr.nops());
+
+ // Convert the lists to a map
+ exmap m;
+ for (lst::const_iterator its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+ m[*its] = *itr;
+
+ // Search for products and powers in the expressions to be substituted
+ // (for an optimization in expairseq::subs())
+ if (is_exactly_a<mul>(*its) || is_exactly_a<power>(*its))
+ options |= subs_options::pattern_is_product;
+ }
+ if (!(options & subs_options::pattern_is_product))
+ options |= subs_options::pattern_is_not_product;
+
+ return bp->subs(m, options);
+}
+
+/** Substitute objects in an expression (syntactic substitution) and return
+ * the result as a new expression. There are two valid types of
+ * replacement arguments: 1) a relational like object==ex and 2) a list of
+ * relationals lst(object1==ex1,object2==ex2,...). */
+ex ex::subs(const ex & e, unsigned options) const
+{
+ if (e.info(info_flags::relation_equal)) {
+ exmap m;
+ const ex & s = e.lhs();
+ m[s] = e.rhs();
+ if (is_exactly_a<mul>(s) || is_exactly_a<power>(s))
+ options |= subs_options::pattern_is_product;
+ return bp->subs(m, options);
+ } else if (!e.info(info_flags::list))
+ throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
+
+ // Convert the list to a map
+ exmap m;
+ GINAC_ASSERT(is_a<lst>(e));
+ for (lst::const_iterator it = ex_to<lst>(e).begin(); it != ex_to<lst>(e).end(); ++it) {
+ ex r = *it;
+ if (!r.info(info_flags::relation_equal))
+ throw(std::invalid_argument("basic::subs(ex): argument must be a list of equations"));
+ const ex & s = r.lhs();
+ m[s] = r.rhs();
+
+ // Search for products and powers in the expressions to be substituted
+ // (for an optimization in expairseq::subs())
+ if (is_exactly_a<mul>(s) || is_exactly_a<power>(s))
+ options |= subs_options::pattern_is_product;
+ }
+ if (!(options & subs_options::pattern_is_product))
+ options |= subs_options::pattern_is_not_product;
+
+ return bp->subs(m, options);
+}
+
/** Traverse expression tree with given visitor, preorder traversal. */
void ex::traverse_preorder(visitor & v) const
{
bool match(const ex & pattern, lst & repl_lst) const { return bp->match(pattern, repl_lst); }
// substitutions
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return bp->subs(ls, lr, options); }
- ex subs(const ex & e, unsigned options = 0) const { return bp->subs(e, options); }
+ ex subs(const exmap & m, unsigned options = 0) const;
+ ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+ ex subs(const ex & e, unsigned options = 0) const;
// function mapping
ex map(map_function & f) const { return bp->map(f); }
inline bool match(const ex & thisex, const ex & pattern, lst & repl_lst)
{ return thisex.match(pattern, repl_lst); }
-inline ex subs(const ex & thisex, const ex & e, unsigned options = 0)
-{ return thisex.subs(e, options); }
-
-inline ex subs(const ex & thisex, const lst & ls, const lst & lr, unsigned options = 0)
-{ return thisex.subs(ls, lr, options); }
-
inline ex simplify_indexed(const ex & thisex)
{ return thisex.simplify_indexed(); }
void operator() (ex &lh, ex &rh) const { lh.swap(rh); }
};
+inline ex ex::subs(const exmap & m, unsigned options) const
+{
+ return bp->subs(m, options);
+}
+
+inline ex subs(const ex & thisex, const exmap & m, unsigned options = 0)
+{ return thisex.subs(m, options); }
+
+inline ex subs(const ex & thisex, const lst & ls, const lst & lr, unsigned options = 0)
+{ return thisex.subs(ls, lr, options); }
+
+inline ex subs(const ex & thisex, const ex & e, unsigned options = 0)
+{ return thisex.subs(e, options); }
+
/* Convert function pointer to function object suitable for map(). */
class pointer_to_map_function : public map_function {
protected:
ex (*ptr)(const ex &);
public:
- explicit pointer_to_map_function(ex (*x)(const ex &)) : ptr(x) {}
+ explicit pointer_to_map_function(ex x(const ex &)) : ptr(x) {}
ex operator()(const ex & e) { return ptr(e); }
};
ex (*ptr)(const ex &, T1);
T1 arg1;
public:
- explicit pointer_to_map_function_1arg(ex (*x)(const ex &, T1), T1 a1) : ptr(x), arg1(a1) {}
+ explicit pointer_to_map_function_1arg(ex x(const ex &, T1), T1 a1) : ptr(x), arg1(a1) {}
ex operator()(const ex & e) { return ptr(e, arg1); }
};
T1 arg1;
T2 arg2;
public:
- explicit pointer_to_map_function_2args(ex (*x)(const ex &, T1, T2), T1 a1, T2 a2) : ptr(x), arg1(a1), arg2(a2) {}
+ explicit pointer_to_map_function_2args(ex x(const ex &, T1, T2), T1 a1, T2 a2) : ptr(x), arg1(a1), arg2(a2) {}
ex operator()(const ex & e) { return ptr(e, arg1, arg2); }
};
T2 arg2;
T3 arg3;
public:
- explicit pointer_to_map_function_3args(ex (*x)(const ex &, T1, T2, T3), T1 a1, T2 a2, T3 a3) : ptr(x), arg1(a1), arg2(a2), arg3(a3) {}
+ explicit pointer_to_map_function_3args(ex x(const ex &, T1, T2, T3), T1 a1, T2 a2, T3 a3) : ptr(x), arg1(a1), arg2(a2), arg3(a3) {}
ex operator()(const ex & e) { return ptr(e, arg1, arg2, arg3); }
};
-inline ex ex::map(ex (*f)(const ex & e)) const
+inline ex ex::map(ex f(const ex &)) const
{
pointer_to_map_function fcn(f);
return bp->map(fcn);
return inherited::match(pattern, repl_lst);
}
-ex expairseq::subs(const lst &ls, const lst &lr, unsigned options) const
+ex expairseq::subs(const exmap & m, unsigned options) const
{
- epvector *vp = subschildren(ls, lr, options);
+ epvector *vp = subschildren(m, options);
if (vp)
return ex_to<basic>(thisexpairseq(vp, overall_coeff));
else if ((options & subs_options::subs_algebraic) && is_exactly_a<mul>(*this))
- return static_cast<const mul *>(this)->algebraic_subs_mul(ls, lr, options);
+ return static_cast<const mul *>(this)->algebraic_subs_mul(m, options);
else
- return subs_one_level(ls, lr, options);
+ return subs_one_level(m, options);
}
// protected
* @see expairseq::subs()
* @return pointer to epvector containing pairs after application of subs,
* or NULL pointer if no members were changed. */
-epvector * expairseq::subschildren(const lst &ls, const lst &lr, unsigned options) const
+epvector * expairseq::subschildren(const exmap & m, unsigned options) const
{
- GINAC_ASSERT(ls.nops()==lr.nops());
-
- // The substitution is "complex" when any of the objects to be substituted
- // is a product or power. In this case we have to recombine the pairs
- // because the numeric coefficients may be part of the search pattern.
- bool complex_subs = false;
- for (lst::const_iterator it = ls.begin(); it != ls.end(); ++it) {
- if (is_exactly_a<mul>(*it) || is_exactly_a<power>(*it)) {
- complex_subs = true;
- break;
+ // When any of the objects to be substituted is a product or power
+ // we have to recombine the pairs because the numeric coefficients may
+ // be part of the search pattern.
+ if (!(options & (subs_options::pattern_is_product | subs_options::pattern_is_not_product))) {
+
+ // Search the list of substitutions and cache our findings
+ for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
+ if (is_exactly_a<mul>(it->first) || is_exactly_a<power>(it->first)) {
+ options |= subs_options::pattern_is_product;
+ break;
+ }
}
+ if (!(options & subs_options::pattern_is_product))
+ options |= subs_options::pattern_is_not_product;
}
- if (complex_subs) {
+ if (options & subs_options::pattern_is_product) {
// Substitute in the recombined pairs
epvector::const_iterator cit = seq.begin(), last = seq.end();
while (cit != last) {
const ex &orig_ex = recombine_pair_to_ex(*cit);
- const ex &subsed_ex = orig_ex.subs(ls, lr, options);
+ const ex &subsed_ex = orig_ex.subs(m, options);
if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
// Something changed, copy seq, subs and return it
// Copy rest
while (cit != last) {
- s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(ls, lr, options)));
+ s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(m, options)));
++cit;
}
return s;
epvector::const_iterator cit = seq.begin(), last = seq.end();
while (cit != last) {
- const ex &subsed_ex = cit->rest.subs(ls, lr, options);
+ const ex &subsed_ex = cit->rest.subs(m, options);
if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
// Something changed, copy seq, subs and return it
// Copy rest
while (cit != last) {
- s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(ls, lr, options),
+ s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(m, options),
cit->coeff));
++cit;
}
ex to_rational(lst &repl_lst) const;
ex to_polynomial(lst &repl_lst) const;
bool match(const ex & pattern, lst & repl_lst) const;
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+ ex subs(const exmap & m, unsigned options = 0) const;
protected:
bool is_equal_same_type(const basic & other) const;
unsigned return_type() const;
bool is_canonical() const;
epvector * expandchildren(unsigned options) const;
epvector * evalchildren(int level) const;
- epvector * subschildren(const lst & ls, const lst & lr, unsigned options = 0) const;
+ epvector * subschildren(const exmap & m, unsigned options = 0) const;
// member variables
class expand_options {
public:
enum {
- expand_trigonometric = 0x0001,
- expand_indexed = 0x0002,
- expand_function_args = 0x0004
+ expand_indexed = 0x0001, ///< expands (a+b).i to a.i+b.i
+ expand_function_args = 0x0002 ///< expands the arguments of functions
};
};
class subs_options {
public:
enum {
- subs_no_pattern = 0x0001,
- subs_algebraic = 0x0002
+ // this should be called "no_pattern"...
+ subs_no_pattern = 0x0001, ///< disable pattern matching
+ // this should be called "algebraic"...
+ subs_algebraic = 0x0002, ///< enable algebraic substitutions
+ pattern_is_product = 0x0004, ///< used internally by expairseq::subschildren()
+ pattern_is_not_product = 0x0008 ///< used internally by expairseq::subschildren()
};
};
class status_flags {
public:
enum {
- dynallocated = 0x0001, ///< Heap-allocated (i.e. created by new if we want to be clever and bypass the stack, @see ex::construct_from_basic() )
- evaluated = 0x0002, ///< .eval() has already done its job
- expanded = 0x0004, ///< .expand(0) has already done its job (other expand() options ignore this flag)
- hash_calculated = 0x0008 ///< .calchash() has already done its job
+ dynallocated = 0x0001, ///< heap-allocated (i.e. created by new if we want to be clever and bypass the stack, @see ex::construct_from_basic() )
+ evaluated = 0x0002, ///< .eval() has already done its job
+ expanded = 0x0004, ///< .expand(0) has already done its job (other expand() options ignore this flag)
+ hash_calculated = 0x0008 ///< .calchash() has already done its job
};
};
return *this;
}
-ex idx::subs(const lst & ls, const lst & lr, unsigned options) const
+ex idx::subs(const exmap & m, unsigned options) const
{
- GINAC_ASSERT(ls.nops() == lr.nops());
-
// First look for index substitutions
- lst::const_iterator its, itr;
- for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
- if (is_equal(ex_to<basic>(*its))) {
-
- // Substitution index->index
- if (is_a<idx>(*itr))
- return *itr;
-
- // Otherwise substitute value
- idx *i_copy = duplicate();
- i_copy->value = *itr;
- i_copy->clearflag(status_flags::hash_calculated);
- return i_copy->setflag(status_flags::dynallocated);
- }
+ exmap::const_iterator it = m.find(*this);
+ if (it != m.end()) {
+
+ // Substitution index->index
+ if (is_a<idx>(it->second))
+ return it->second;
+
+ // Otherwise substitute value
+ idx *i_copy = duplicate();
+ i_copy->value = it->second;
+ i_copy->clearflag(status_flags::hash_calculated);
+ return i_copy->setflag(status_flags::dynallocated);
}
// None, substitute objects in value (not in dimension)
- const ex &subsed_value = value.subs(ls, lr, options);
+ const ex &subsed_value = value.subs(m, options);
if (are_ex_trivially_equal(value, subsed_value))
return *this;
ex op(size_t i) const;
ex map(map_function & f) const;
ex evalf(int level = 0) const;
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+ ex subs(const exmap & m, unsigned options = 0) const;
protected:
ex derivative(const symbol & s) const;
status_flags::evaluated);
}
-ex matrix::subs(const lst & ls, const lst & lr, unsigned options) const
-{
+ex matrix::subs(const exmap & mp, unsigned options) const
+{
exvector m2(row * col);
for (unsigned r=0; r<row; ++r)
for (unsigned c=0; c<col; ++c)
- m2[r*col+c] = m[r*col+c].subs(ls, lr, options);
+ m2[r*col+c] = m[r*col+c].subs(mp, options);
- return matrix(row, col, m2).subs_one_level(ls, lr, options);
+ return matrix(row, col, m2).subs_one_level(mp, options);
}
// protected
ex & let_op(size_t i);
ex eval(int level=0) const;
ex evalm() const {return *this;}
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+ ex subs(const exmap & m, unsigned options = 0) const;
ex eval_indexed(const basic & i) const;
ex add_indexed(const ex & self, const ex & other) const;
ex scalar_mul_indexed(const ex & self, const numeric & other) const;
return true;
}
-ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) const
+ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const
{
std::vector<bool> subsed(seq.size(), false);
exvector subsresult(seq.size());
- lst::const_iterator its, itr;
- for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+ for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
- if (is_exactly_a<mul>(*its)) {
+ if (is_exactly_a<mul>(it->first)) {
int nummatches = std::numeric_limits<int>::max();
std::vector<bool> currsubsed(seq.size(), false);
bool succeed = true;
lst repls;
- for (size_t j=0; j<its->nops(); j++) {
+ for (size_t j=0; j<it->first.nops(); j++) {
bool found=false;
for (size_t k=0; k<nops(); k++) {
if (currsubsed[k] || subsed[k])
continue;
- if (tryfactsubs(op(k), its->op(j), nummatches, repls)) {
+ if (tryfactsubs(op(k), it->first.op(j), nummatches, repls)) {
currsubsed[k] = true;
found = true;
break;
subsresult[j] = op(j);
else {
foundfirstsubsedfactor = true;
- subsresult[j] = op(j) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches);
+ subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::subs_no_pattern) / it->first.subs(ex(repls), subs_options::subs_no_pattern), nummatches);
}
subsed[j] = true;
}
lst repls;
for (size_t j=0; j<this->nops(); j++) {
- if (!subsed[j] && tryfactsubs(op(j), *its, nummatches, repls)) {
+ if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)) {
subsed[j] = true;
- subsresult[j] = op(j) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches);
+ subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::subs_no_pattern) / it->first.subs(ex(repls), subs_options::subs_no_pattern), nummatches);
}
}
}
}
}
if (!subsfound)
- return subs_one_level(ls, lr, options | subs_options::subs_algebraic);
+ return subs_one_level(m, options | subs_options::subs_algebraic);
exvector ev; ev.reserve(nops());
for (size_t i=0; i<nops(); i++) {
ex evalf(int level=0) const;
ex evalm() const;
ex series(const relational & s, int order, unsigned options = 0) const;
- ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+ ex normal(exmap & repl, int level = 0) const;
numeric integer_content() const;
ex smod(const numeric &xi) const;
numeric max_coefficient() const;
// non-virtual functions in this class
public:
- ex algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) const;
+ ex algebraic_subs_mul(const exmap & m, unsigned options) const;
protected:
epvector * expandchildren(unsigned options) const;
};
/** Create a symbol for replacing the expression "e" (or return a previously
- * assigned symbol). The symbol is appended to sym_lst and returned, the
- * expression is appended to repl_lst.
+ * assigned symbol). The symbol and expression are appended to repl, for
+ * a later application of subs().
* @see ex::normal */
-static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst)
+static ex replace_with_symbol(const ex & e, exmap & repl)
{
- // Expression already in repl_lst? Then return the assigned symbol
- lst::const_iterator its, itr;
- for (its = sym_lst.begin(), itr = repl_lst.begin(); itr != repl_lst.end(); ++its, ++itr)
- if (itr->is_equal(e))
- return *its;
+ // Expression already in repl? Then return the assigned symbol
+ for (exmap::const_iterator it = repl.begin(); it != repl.end(); ++it)
+ if (it->second.is_equal(e))
+ return it->first;
// Otherwise create new symbol and add to list, taking care that the
- // replacement expression doesn't contain symbols from the sym_lst
+ // replacement expression doesn't itself contain symbols from repl,
// because subs() is not recursive
- symbol s;
- ex es(s);
- ex e_replaced = e.subs(sym_lst, repl_lst);
- sym_lst.append(es);
- repl_lst.append(e_replaced);
+ ex es = (new symbol)->setflag(status_flags::dynallocated);
+ ex e_replaced = e.subs(repl);
+ repl[es] = e_replaced;
return es;
}
* to repl_lst and the symbol is returned.
* @see basic::to_rational
* @see basic::to_polynomial */
-static ex replace_with_symbol(const ex &e, lst &repl_lst)
+static ex replace_with_symbol(const ex & e, lst & repl_lst)
{
// Expression already in repl_lst? Then return the assigned symbol
for (lst::const_iterator it = repl_lst.begin(); it != repl_lst.end(); ++it)
return it->op(0);
// Otherwise create new symbol and add to list, taking care that the
- // replacement expression doesn't contain symbols from the sym_lst
+ // replacement expression doesn't itself contain symbols from the repl_lst,
// because subs() is not recursive
- symbol s;
- ex es(s);
+ ex es = (new symbol)->setflag(status_flags::dynallocated);
ex e_replaced = e.subs(repl_lst);
repl_lst.append(es == e_replaced);
return es;
/** Default implementation of ex::normal(). It normalizes the children and
* replaces the object with a temporary symbol.
* @see ex::normal */
-ex basic::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex basic::normal(exmap & repl, int level) const
{
if (nops() == 0)
- return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
else {
if (level == 1)
- return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
else if (level == -max_recursion_level)
throw(std::runtime_error("max recursion level reached"));
else {
normal_map_function map_normal(level - 1);
- return (new lst(replace_with_symbol(map(map_normal), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(map(map_normal), repl), _ex1))->setflag(status_flags::dynallocated);
}
}
}
/** Implementation of ex::normal() for symbols. This returns the unmodified symbol.
* @see ex::normal */
-ex symbol::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex symbol::normal(exmap & repl, int level) const
{
return (new lst(*this, _ex1))->setflag(status_flags::dynallocated);
}
* into re+I*im and replaces I and non-rational real numbers with a temporary
* symbol.
* @see ex::normal */
-ex numeric::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex numeric::normal(exmap & repl, int level) const
{
numeric num = numer();
ex numex = num;
if (num.is_real()) {
if (!num.is_integer())
- numex = replace_with_symbol(numex, sym_lst, repl_lst);
+ numex = replace_with_symbol(numex, repl);
} else { // complex
numeric re = num.real(), im = num.imag();
- ex re_ex = re.is_rational() ? re : replace_with_symbol(re, sym_lst, repl_lst);
- ex im_ex = im.is_rational() ? im : replace_with_symbol(im, sym_lst, repl_lst);
- numex = re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst);
+ ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl);
+ ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl);
+ numex = re_ex + im_ex * replace_with_symbol(I, repl);
}
// 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(lst &sym_lst, lst &repl_lst, int level) const
+ex add::normal(exmap & repl, int level) const
{
if (level == 1)
- return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
else if (level == -max_recursion_level)
throw(std::runtime_error("max recursion level reached"));
dens.reserve(seq.size()+1);
epvector::const_iterator it = seq.begin(), itend = seq.end();
while (it != itend) {
- ex n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(sym_lst, repl_lst, level-1);
+ ex n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(repl, level-1);
nums.push_back(n.op(0));
dens.push_back(n.op(1));
it++;
}
- ex n = ex_to<numeric>(overall_coeff).normal(sym_lst, repl_lst, level-1);
+ ex n = ex_to<numeric>(overall_coeff).normal(repl, level-1);
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(lst &sym_lst, lst &repl_lst, int level) const
+ex mul::normal(exmap & repl, int level) const
{
if (level == 1)
- return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
else if (level == -max_recursion_level)
throw(std::runtime_error("max recursion level reached"));
ex n;
epvector::const_iterator it = seq.begin(), itend = seq.end();
while (it != itend) {
- n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(sym_lst, repl_lst, level-1);
+ n = ex_to<basic>(recombine_pair_to_ex(*it)).normal(repl, level-1);
num.push_back(n.op(0));
den.push_back(n.op(1));
it++;
}
- n = ex_to<numeric>(overall_coeff).normal(sym_lst, repl_lst, level-1);
+ n = ex_to<numeric>(overall_coeff).normal(repl, level-1);
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(lst &sym_lst, lst &repl_lst, int level) const
+ex power::normal(exmap & repl, int level) const
{
if (level == 1)
- return (new lst(replace_with_symbol(*this, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(*this, repl), _ex1))->setflag(status_flags::dynallocated);
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(sym_lst, repl_lst, level-1);
- ex n_exponent = ex_to<basic>(exponent).normal(sym_lst, repl_lst, level-1);
+ ex n_basis = ex_to<basic>(basis).normal(repl, level-1);
+ ex n_exponent = ex_to<basic>(exponent).normal(repl, level-1);
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 (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), repl), _ex1))->setflag(status_flags::dynallocated);
} else if (n_exponent.info(info_flags::negative)) {
if (n_basis.op(1).is_equal(_ex1)) {
// a^-x -> {1, sym(a^x)}
- return (new lst(_ex1, replace_with_symbol(power(n_basis.op(0), -n_exponent), sym_lst, repl_lst)))->setflag(status_flags::dynallocated);
+ return (new lst(_ex1, replace_with_symbol(power(n_basis.op(0), -n_exponent), repl)))->setflag(status_flags::dynallocated);
} else {
// (a/b)^-x -> {sym((b/a)^x), 1}
- return (new lst(replace_with_symbol(power(n_basis.op(1) / n_basis.op(0), -n_exponent), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(power(n_basis.op(1) / n_basis.op(0), -n_exponent), repl), _ex1))->setflag(status_flags::dynallocated);
}
}
}
// (a/b)^x -> {sym((a/b)^x, 1}
- return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), repl), _ex1))->setflag(status_flags::dynallocated);
}
/** Implementation of ex::normal() for pseries. It normalizes each coefficient
* and replaces the series by a temporary symbol.
* @see ex::normal */
-ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const
+ex pseries::normal(exmap & repl, int level) const
{
epvector newseq;
epvector::const_iterator i = seq.begin(), end = seq.end();
++i;
}
ex n = pseries(relational(var,point), newseq);
- return (new lst(replace_with_symbol(n, sym_lst, repl_lst), _ex1))->setflag(status_flags::dynallocated);
+ return (new lst(replace_with_symbol(n, repl), _ex1))->setflag(status_flags::dynallocated);
}
* @return normalized expression */
ex ex::normal(int level) const
{
- lst sym_lst, repl_lst;
+ exmap repl;
- ex e = bp->normal(sym_lst, repl_lst, level);
+ ex e = bp->normal(repl, level);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
- if (sym_lst.nops() > 0)
- e = e.subs(sym_lst, repl_lst);
+ if (!repl.empty())
+ e = e.subs(repl);
// Convert {numerator, denominator} form back to fraction
return e.op(0) / e.op(1);
* @return numerator */
ex ex::numer() const
{
- lst sym_lst, repl_lst;
+ exmap repl;
- ex e = bp->normal(sym_lst, repl_lst, 0);
+ ex e = bp->normal(repl, 0);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
- if (sym_lst.nops() > 0)
- return e.op(0).subs(sym_lst, repl_lst);
- else
+ if (repl.empty())
return e.op(0);
+ else
+ return e.op(0).subs(repl);
}
/** Get denominator of an expression. If the expression is not of the normal
* @return denominator */
ex ex::denom() const
{
- lst sym_lst, repl_lst;
+ exmap repl;
- ex e = bp->normal(sym_lst, repl_lst, 0);
+ ex e = bp->normal(repl, 0);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
- if (sym_lst.nops() > 0)
- return e.op(1).subs(sym_lst, repl_lst);
- else
+ if (repl.empty())
return e.op(1);
+ else
+ return e.op(1).subs(repl);
}
/** Get numerator and denominator of an expression. If the expresison is not
* @return a list [numerator, denominator] */
ex ex::numer_denom() const
{
- lst sym_lst, repl_lst;
+ exmap repl;
- ex e = bp->normal(sym_lst, repl_lst, 0);
+ ex e = bp->normal(repl, 0);
GINAC_ASSERT(is_a<lst>(e));
// Re-insert replaced symbols
- if (sym_lst.nops() > 0)
- return e.subs(sym_lst, repl_lst);
- else
+ if (repl.empty())
return e;
+ else
+ return e.subs(repl);
}
bool has(const ex &other) const;
ex eval(int level = 0) const;
ex evalf(int level = 0) const;
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return subs_one_level(ls, lr, options); } // overwrites basic::subs() for performance reasons
- ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+ ex subs(const exmap & m, unsigned options = 0) const { return subs_one_level(m, options); } // overwrites basic::subs() for performance reasons
+ ex normal(exmap & repl, int level = 0) const;
ex to_rational(lst &repl_lst) const;
ex to_polynomial(lst &repl_lst) const;
numeric integer_content() const;
// from mul.cpp
extern bool tryfactsubs(const ex &, const ex &, int &, lst &);
-ex power::subs(const lst & ls, const lst & lr, unsigned options) const
+ex power::subs(const exmap & m, unsigned options) const
{
- const ex &subsed_basis = basis.subs(ls, lr, options);
- const ex &subsed_exponent = exponent.subs(ls, lr, options);
+ const ex &subsed_basis = basis.subs(m, options);
+ const ex &subsed_exponent = exponent.subs(m, options);
if (!are_ex_trivially_equal(basis, subsed_basis)
|| !are_ex_trivially_equal(exponent, subsed_exponent))
- return power(subsed_basis, subsed_exponent).subs_one_level(ls, lr, options);
+ return power(subsed_basis, subsed_exponent).subs_one_level(m, options);
if (!(options & subs_options::subs_algebraic))
- return subs_one_level(ls, lr, options);
+ return subs_one_level(m, options);
- lst::const_iterator its, itr;
- for (its = ls.begin(), itr = lr.begin(); its != ls.end(); ++its, ++itr) {
+ for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
int nummatches = std::numeric_limits<int>::max();
lst repls;
- if (tryfactsubs(*this, *its, nummatches, repls))
- return (ex_to<basic>((*this) * power(itr->subs(ex(repls), subs_options::subs_no_pattern) / its->subs(ex(repls), subs_options::subs_no_pattern), nummatches))).subs_one_level(ls, lr, options);
+ if (tryfactsubs(*this, it->first, nummatches, repls))
+ return (ex_to<basic>((*this) * power(it->second.subs(ex(repls), subs_options::subs_no_pattern) / it->first.subs(ex(repls), subs_options::subs_no_pattern), nummatches))).subs_one_level(m, options);
}
- return subs_one_level(ls, lr, options);
+ return subs_one_level(m, options);
}
ex power::eval_ncmul(const exvector & v) const
ex evalf(int level=0) const;
ex evalm() const;
ex series(const relational & s, int order, unsigned options = 0) const;
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
- ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+ ex subs(const exmap & m, unsigned options = 0) const;
+ ex normal(exmap & repl, int level = 0) const;
ex to_rational(lst &repl_lst) const;
ex to_polynomial(lst &repl_lst) const;
exvector get_free_indices() const;
return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
}
-ex pseries::subs(const lst & ls, const lst & lr, unsigned options) const
+ex pseries::subs(const exmap & m, unsigned options) const
{
// If expansion variable is being substituted, convert the series to a
// polynomial and do the substitution there because the result might
// no longer be a power series
- if (ls.has(var))
- return convert_to_poly(true).subs(ls, lr, options);
+ if (m.find(var) != m.end())
+ return convert_to_poly(true).subs(m, options);
// Otherwise construct a new series with substituted coefficients and
// expansion point
newseq.reserve(seq.size());
epvector::const_iterator it = seq.begin(), itend = seq.end();
while (it != itend) {
- newseq.push_back(expair(it->rest.subs(ls, lr, options), it->coeff));
+ newseq.push_back(expair(it->rest.subs(m, options), it->coeff));
++it;
}
- return (new pseries(relational(var,point.subs(ls, lr, options)), newseq))->setflag(status_flags::dynallocated);
+ return (new pseries(relational(var,point.subs(m, options)), newseq))->setflag(status_flags::dynallocated);
}
/** Implementation of ex::expand() for a power series. It expands all the
ex eval(int level=0) const;
ex evalf(int level=0) const;
ex series(const relational & r, int order, unsigned options = 0) const;
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
- ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+ ex subs(const exmap & m, unsigned options = 0) const;
+ ex normal(exmap & repl, int level = 0) const;
ex expand(unsigned options = 0) const;
protected:
ex derivative(const symbol & s) const;
return (new relational(lh.eval(level-1),rh.eval(level-1),o))->setflag(status_flags::dynallocated | status_flags::evaluated);
}
-ex relational::subs(const lst & ls, const lst & lr, unsigned options) const
+ex relational::subs(const exmap & m, unsigned options) const
{
- const ex & subsed_lh = lh.subs(ls, lr, options);
- const ex & subsed_rh = rh.subs(ls, lr, options);
+ const ex & subsed_lh = lh.subs(m, options);
+ const ex & subsed_rh = rh.subs(m, options);
if (!are_ex_trivially_equal(lh, subsed_lh) || !are_ex_trivially_equal(rh, subsed_rh))
- return relational(subsed_lh, subsed_rh, o).subs_one_level(ls, lr, options);
+ return relational(subsed_lh, subsed_rh, o).subs_one_level(m, options);
else
- return subs_one_level(ls, lr, options);
+ return subs_one_level(m, options);
}
ex relational::eval_ncmul(const exvector & v) const
size_t nops() const;
ex op(size_t i) const;
ex map(map_function & f) const;
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+ ex subs(const exmap & m, unsigned options = 0) const;
ex eval(int level=0) const;
protected:
public:
// substitutions
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return inherited::subs(ls, lr, options); }
+ ex subs(const exmap & m, unsigned options = 0) const { return inherited::subs(m, options); }
// function mapping
ex map(map_function & f) const { return inherited::map(f); }
ex series(const relational & r, int order, unsigned options = 0) const { return inherited::series(r, order, options); }
// rational functions
- ex normal(lst & sym_lst, lst & repl_lst, int level = 0) const { return inherited::normal(sym_lst, repl_lst, level); }
+ ex normal(exmap & repl, int level = 0) const { return inherited::normal(repl, level); }
ex to_rational(lst & repl_lst) const { return inherited::to_rational(repl_lst); }
ex to_polynomial(lst & repl_lst) const { return inherited::to_polynomial(repl_lst); }
ex eval(int level = 0) const;
ex evalf(int level = 0) const { return *this; } // overwrites basic::evalf() for performance reasons
ex series(const relational & s, int order, unsigned options = 0) const;
- ex subs(const lst & ls, const lst & lr, unsigned options = 0) const { return subs_one_level(ls, lr, options); } // overwrites basic::subs() for performance reasons
- ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
+ ex subs(const exmap & m, unsigned options = 0) const { return subs_one_level(m, options); } // overwrites basic::subs() for performance reasons
+ ex normal(exmap & repl, int level = 0) const;
ex to_rational(lst &repl_lst) const;
ex to_polynomial(lst &repl_lst) const;
protected:
// dimension
if (!i1.get_dim().is_equal(i2.get_dim())) {
ex min_dim = i1.minimal_dim(i2);
- return i.subs(lst(i1 == i1.replace_dim(min_dim), i2 == i2.replace_dim(min_dim)));
+ exmap m;
+ m[i1] = i1.replace_dim(min_dim);
+ m[i2] = i2.replace_dim(min_dim);
+ return i.subs(m);
}
// Trace of delta tensor is the (effective) dimension of the space
// dimension
if (!i1.get_dim().is_equal(i2.get_dim())) {
ex min_dim = i1.minimal_dim(i2);
- return i.subs(lst(i1 == i1.replace_dim(min_dim), i2 == i2.replace_dim(min_dim)));
+ exmap m;
+ m[i1] = i1.replace_dim(min_dim);
+ m[i2] = i2.replace_dim(min_dim);
+ return i.subs(m);
}
// A metric tensor with one covariant and one contravariant index gets
git://www.ginac.de / ginac.git / commitdiff