From 4275fab61d542fc9b67ef5489361830434693b71 Mon Sep 17 00:00:00 2001 From: Christian Bauer Date: Wed, 4 Dec 2002 20:44:10 +0000 Subject: [PATCH] implemented 'smartsubs' algebraic substitution, donated by Chris Dams (tutorial not yet in sync, will update later) --- doc/tutorial/ginac.texi | 68 ++++++++++++++++++++-- ginac/basic.cpp | 12 ++-- ginac/basic.h | 4 +- ginac/container.pl | 20 +++---- ginac/ex.h | 12 ++-- ginac/expairseq.cpp | 24 ++++---- ginac/expairseq.h | 4 +- ginac/flags.h | 10 ++++ ginac/idx.cpp | 4 +- ginac/idx.h | 2 +- ginac/matrix.cpp | 8 +-- ginac/matrix.h | 2 +- ginac/mul.cpp | 126 ++++++++++++++++++++++++++++++++++++++++ ginac/mul.h | 2 + ginac/power.cpp | 24 ++++++-- ginac/power.h | 2 +- ginac/pseries.cpp | 8 +-- ginac/pseries.h | 2 +- 18 files changed, 275 insertions(+), 59 deletions(-) diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index 113cf491..0b9acc36 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -734,7 +734,7 @@ containers of expressions and so on. @cindex container @cindex atom -To get an idea about what kinds of symbolic composits may be built we +To get an idea about what kinds of symbolic composites may be built we have a look at the most important classes in the class hierarchy and some of the relations among the classes: @@ -895,7 +895,7 @@ can use the expression's @code{.subs()} method (@pxref{Substituting Expressions} For storing numerical things, GiNaC uses Bruno Haible's library CLN. The classes therein serve as foundation classes for GiNaC. CLN stands for Class Library for Numbers or alternatively for Common Lisp Numbers. -In order to find out more about CLN's internals the reader is refered to +In order to find out more about CLN's internals, the reader is referred to the documentation of that library. @inforef{Introduction, , cln}, for more information. Suffice to say that it is by itself build on top of another library, the GNU Multiple Precision library GMP, which is an @@ -1009,7 +1009,7 @@ naively expect it to be rounded in the decimal system. But note also, that in both cases you got a couple of extra digits. This is because numbers are internally stored by CLN as chunks of binary digits in order to match your machine's word size and to not waste precision. Thus, on -architectures with differnt word size, the above output might even +architectures with different word size, the above output might even differ with regard to actually computed digits. It should be clear that objects of class @code{numeric} should be used @@ -3061,6 +3061,7 @@ Notes: are also valid patterns. @end itemize +@subsection Matching expressions @cindex @code{match()} The most basic application of patterns is to check whether an expression matches a given pattern. This is done by the function @@ -3170,6 +3171,7 @@ FAIL @{$0==x^2@} @end example +@subsection Matching parts of expressions @cindex @code{has()} A more general way to look for patterns in expressions is provided by the member function @@ -3238,6 +3240,7 @@ sin(y)*a+sin(x)*b+sin(x)*a+sin(y)*b @{sin(y),sin(x)@} @end example +@subsection Substituting expressions @cindex @code{subs()} Probably the most useful application of patterns is to use them for substituting expressions with the @code{subs()} method. Wildcards can be @@ -3280,6 +3283,63 @@ The last example would be written in C++ in this way: @} @end example +@subsection Algebraic substitutions +This subsection is not written by one of the original authors of GiNaC but +by Chris Dams. All errors in the here described feature can be attributed to +him (but mailed to the GiNaC mailing list, which he reads). + +The @code{subs()} method has an extra, optional, argument. This +argument can be used to pass one of the @code{subs_options} to it. +The only option that is currently +available is the @code{subs_smart} option and it affects multiplications +and powers. If you want to substitute some factors of a multiplication, you +only need to list these factors in your pattern. Furthermore if a(n) (integer) +power of some expression occurs in your pattern and in the expression that +you want the substitution to occur in, it can be substituted as many times +as possible, without getting negative powers. + +An example clarifies it all (hopefully). + +@example +cout << (a*a*a*a+b*b*b*b+pow(x+y,4)).subs(wild()*wild()==pow(wild(),3), + subs_options::subs_smart) << endl; +// --> (y+x)^6+b^6+a^6 + +cout << ((a+b+c)*(a+b+c)).subs(a+b==x,subs_options::subs_smart) << endl; +// --> (c+b+a)^2 +// Powers and multiplications are smart, but addition is just the same. + +cout << ((a+b+c)*(a+b+c)).subs(a+b+wild()==x+wild(), subs_options::subs_smart) + << endl; +// --> (x+c)^2 +// As I said: addition is just the same. + +cout << (pow(a,5)*pow(b,7)+2*b).subs(b*b*a==x,subs_options::subs_smart) << endl; +// --> x^3*b*a^2+2*b + +cout << (pow(a,-5)*pow(b,-7)+2*b).subs(1/(b*b*a)==x,subs_options::subs_smart) + << endl; +// --> 2*b+x^3*b^(-1)*a^(-2) + +cout << (4*x*x*x-2*x*x+5*x-1).subs(x==a,subs_options::subs_smart) + << endl; +// --> -1-2*a^2+4*a^3+5*a + +cout << (4*x*x*x-2*x*x+5*x-1).subs(pow(x,wild())==pow(a,wild()), + subs_options::subs_smart) << endl; +// --> -1+5*x+4*x^3-2*x^2 +// You should not really need this kind of patterns very often now. +// But perhaps this it's-not-a-bug-it's-a-feature (c/sh)ould still change. + +cout << ex(sin(1+sin(x))).subs(sin(wild())==cos(wild()), + subs_options::subs_smart) << endl; +// --> cos(1+cos(x)) + +cout << expand((a*sin(x+y)*sin(x+y)+a*cos(x+y)*cos(x+y)+b) + .subs((pow(cos(wild()),2)==1-pow(sin(wild()),2)), + subs_options::subs_smart)) << endl; +// --> b+a +@end example @node Applying a Function on Subexpressions, Polynomial Arithmetic, Pattern Matching and Advanced Substitutions, Methods and Functions @c node-name, next, previous, up @@ -5793,7 +5853,7 @@ simple_SOURCES = simple.cpp @end example This @file{Makefile.am}, says that we are building a single executable, -from a single sourcefile @file{simple.cpp}. Since every program +from a single source file @file{simple.cpp}. Since every program we are building uses GiNaC we simply added the GiNaC options to @env{$LIBS} and @env{$CPPFLAGS}, but in other circumstances, we might want to specify them on a per-program basis: for instance by diff --git a/ginac/basic.cpp b/ginac/basic.cpp index 26db7b34..bdf37987 100644 --- a/ginac/basic.cpp +++ b/ginac/basic.cpp @@ -503,11 +503,11 @@ bool basic::match(const ex & pattern, lst & repl_lst) const /** Substitute a set of objects by arbitrary expressions. The ex returned * will already be evaluated. */ -ex basic::subs(const lst & ls, const lst & lr, bool no_pattern) const +ex basic::subs(const lst & ls, const lst & lr, unsigned options) const { GINAC_ASSERT(ls.nops() == lr.nops()); - if (no_pattern) { + if (options & subs_options::subs_no_pattern) { for (unsigned i=0; i(ls.op(i)))) return lr.op(i); @@ -516,7 +516,7 @@ ex basic::subs(const lst & ls, const lst & lr, bool no_pattern) const for (unsigned i=0; i(ls.op(i)), repl_lst)) - return lr.op(i).subs(repl_lst, true); // avoid infinite recursion when re-substituting the wildcards + return lr.op(i).subs(repl_lst, options | subs_options::subs_no_pattern); // avoid infinite recursion when re-substituting the wildcards } } @@ -684,10 +684,10 @@ ex basic::expand(unsigned options) const * 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, bool no_pattern) const +ex basic::subs(const ex & e, unsigned options) const { if (e.info(info_flags::relation_equal)) { - return subs(lst(e), no_pattern); + return subs(lst(e), options); } if (!e.info(info_flags::list)) { throw(std::invalid_argument("basic::subs(ex): argument must be a list")); @@ -702,7 +702,7 @@ ex basic::subs(const ex & e, bool no_pattern) const ls.append(r.op(0)); lr.append(r.op(1)); } - return subs(ls, lr, no_pattern); + return subs(ls, lr, options); } /** Compare objects syntactically to establish canonical ordering. diff --git a/ginac/basic.h b/ginac/basic.h index cd97a835..fe99c88c 100644 --- a/ginac/basic.h +++ b/ginac/basic.h @@ -112,7 +112,7 @@ public: // only const functions please (may break reference counting) virtual ex evalm(void) const; virtual ex series(const relational & r, int order, unsigned options = 0) const; virtual bool match(const ex & pattern, lst & repl_lst) const; - virtual ex subs(const lst & ls, const lst & lr, bool no_pattern = false) const; + virtual ex subs(const lst & ls, const lst & lr, unsigned options = 0) const; virtual ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const; virtual ex to_rational(lst &repl_lst) const; virtual numeric integer_content(void) const; @@ -135,7 +135,7 @@ protected: // functions that should be called from class ex only // non-virtual functions in this class public: - ex subs(const ex & e, bool no_pattern = false) const; + ex subs(const ex & e, unsigned options = 0) const; ex diff(const symbol & s, unsigned nth=1) const; int compare(const basic & other) const; bool is_equal(const basic & other) const; diff --git a/ginac/container.pl b/ginac/container.pl index 2860ccf1..44703735 100755 --- a/ginac/container.pl +++ b/ginac/container.pl @@ -241,7 +241,7 @@ public: ex & let_op(int i); ex map(map_function & f) const; ex eval(int level=0) const; - ex subs(const lst & ls, const lst & lr, bool no_pattern = false) const; + ex subs(const lst & ls, const lst & lr, unsigned options = 0) const; protected: bool is_equal_same_type(const basic & other) const; unsigned return_type(void) const; @@ -263,7 +263,7 @@ protected: protected: bool is_canonical() const; ${STLT} evalchildren(int level) const; - ${STLT} * subschildren(const lst & ls, const lst & lr, bool no_pattern = false) const; + ${STLT} * subschildren(const lst & ls, const lst & lr, unsigned options = 0) const; protected: ${STLT} seq; @@ -474,13 +474,13 @@ ex ${CONTAINER}::eval(int level) const return this${CONTAINER}(evalchildren(level)); } -ex ${CONTAINER}::subs(const lst & ls, const lst & lr, bool no_pattern) const -{ - ${STLT} *vp = subschildren(ls, lr, no_pattern); +ex ${CONTAINER}::subs(const lst & ls, const lst & lr, unsigned options) const +{ + ${STLT} *vp = subschildren(ls, lr, options); if (vp) - return ex_to(this${CONTAINER}(vp)).basic::subs(ls, lr, no_pattern); + return ex_to(this${CONTAINER}(vp)).basic::subs(ls, lr, options); else - return basic::subs(ls, lr, no_pattern); + return basic::subs(ls, lr, options); } // protected @@ -641,7 +641,7 @@ ${STLT} ${CONTAINER}::evalchildren(int level) const return s; } -${STLT} * ${CONTAINER}::subschildren(const lst & ls, const lst & lr, bool no_pattern) const +${STLT} * ${CONTAINER}::subschildren(const lst & ls, const lst & lr, unsigned options) const { // returns a NULL pointer if nothing had to be substituted // returns a pointer to a newly created epvector otherwise @@ -649,7 +649,7 @@ ${STLT} * ${CONTAINER}::subschildren(const lst & ls, const lst & lr, bool no_pat ${STLT}::const_iterator cit = seq.begin(), end = seq.end(); while (cit != end) { - const ex & subsed_ex = cit->subs(ls, lr, no_pattern); + const ex & subsed_ex = cit->subs(ls, lr, options); if (!are_ex_trivially_equal(*cit, subsed_ex)) { // something changed, copy seq, subs and return it @@ -669,7 +669,7 @@ ${STLT} * ${CONTAINER}::subschildren(const lst & ls, const lst & lr, bool no_pat // copy rest while (cit2 != end) { - s->push_back(cit2->subs(ls, lr, no_pattern)); + s->push_back(cit2->subs(ls, lr, options)); ++cit2; } return s; diff --git a/ginac/ex.h b/ginac/ex.h index dfab1461..2432666c 100644 --- a/ginac/ex.h +++ b/ginac/ex.h @@ -146,8 +146,8 @@ public: ex series(const ex & r, int order, unsigned options = 0) const; bool match(const ex & pattern) const; bool match(const ex & pattern, lst & repl_lst) const { return bp->match(pattern, repl_lst); } - ex subs(const lst & ls, const lst & lr, bool no_pattern = false) const { return bp->subs(ls, lr, no_pattern); } - ex subs(const ex & e, bool no_pattern = false) const { return bp->subs(e, no_pattern); } + 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); } exvector get_free_indices(void) const { return bp->get_free_indices(); } ex simplify_indexed(void) const; ex simplify_indexed(const scalar_products & sp) const; @@ -425,11 +425,11 @@ inline ex series(const ex & thisex, const ex & r, int order, unsigned options = 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) -{ return thisex.subs(e); } +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) -{ return thisex.subs(ls, lr); } +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(); } diff --git a/ginac/expairseq.cpp b/ginac/expairseq.cpp index 51af5f88..c0a14850 100644 --- a/ginac/expairseq.cpp +++ b/ginac/expairseq.cpp @@ -401,13 +401,16 @@ found: ; return inherited::match(pattern, repl_lst); } -ex expairseq::subs(const lst &ls, const lst &lr, bool no_pattern) const -{ - epvector *vp = subschildren(ls, lr, no_pattern); +ex expairseq::subs(const lst &ls, const lst &lr, unsigned options) const +{ + if ((options & subs_options::subs_algebraic) && is_exactly_a(*this)) + return static_cast(this)->algebraic_subs_mul(ls, lr, options); + + epvector *vp = subschildren(ls, lr, options); if (vp) return ex_to(thisexpairseq(vp, overall_coeff)); else - return basic::subs(ls, lr, no_pattern); + return basic::subs(ls, lr, options); } // protected @@ -1559,7 +1562,7 @@ epvector * expairseq::evalchildren(int level) const * @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, bool no_pattern) const +epvector * expairseq::subschildren(const lst &ls, const lst &lr, unsigned options) const { GINAC_ASSERT(ls.nops()==lr.nops()); @@ -1567,11 +1570,12 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, bool no_pattern // 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 (unsigned i=0; i(ls.op(i)) || is_exactly_a(ls.op(i))) { complex_subs = true; break; } + } if (complex_subs) { @@ -1580,7 +1584,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, bool no_pattern while (cit != last) { const ex &orig_ex = recombine_pair_to_ex(*cit); - const ex &subsed_ex = orig_ex.subs(ls, lr, no_pattern); + const ex &subsed_ex = orig_ex.subs(ls, lr, options); if (!are_ex_trivially_equal(orig_ex, subsed_ex)) { // Something changed, copy seq, subs and return it @@ -1596,7 +1600,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, bool no_pattern // Copy rest while (cit != last) { - s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(ls, lr, no_pattern))); + s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(ls, lr, options))); ++cit; } return s; @@ -1611,7 +1615,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, bool no_pattern epvector::const_iterator cit = seq.begin(), last = seq.end(); while (cit != last) { - const ex &subsed_ex = cit->rest.subs(ls, lr, no_pattern); + const ex &subsed_ex = cit->rest.subs(ls, lr, options); if (!are_ex_trivially_equal(cit->rest, subsed_ex)) { // Something changed, copy seq, subs and return it @@ -1627,7 +1631,7 @@ epvector * expairseq::subschildren(const lst &ls, const lst &lr, bool no_pattern // Copy rest while (cit != last) { - s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(ls, lr, no_pattern), + s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(ls, lr, options), cit->coeff)); ++cit; } diff --git a/ginac/expairseq.h b/ginac/expairseq.h index 23534978..dc519f10 100644 --- a/ginac/expairseq.h +++ b/ginac/expairseq.h @@ -96,7 +96,7 @@ public: ex eval(int level=0) const; ex to_rational(lst &repl_lst) const; bool match(const ex & pattern, lst & repl_lst) const; - ex subs(const lst & ls, const lst & lr, bool no_pattern = false) const; + ex subs(const lst & ls, const lst & lr, unsigned options = 0) const; protected: int compare_same_type(const basic & other) const; bool is_equal_same_type(const basic & other) const; @@ -163,7 +163,7 @@ protected: bool is_canonical() const; epvector * expandchildren(unsigned options) const; epvector * evalchildren(int level) const; - epvector * subschildren(const lst & ls, const lst & lr, bool no_pattern = false) const; + epvector * subschildren(const lst & ls, const lst & lr, unsigned options = 0) const; // member variables diff --git a/ginac/flags.h b/ginac/flags.h index 5dd9ae0d..ac40e6e4 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -25,6 +25,7 @@ namespace GiNaC { +/** Flags to control the behavior of expand(). */ class expand_options { public: enum { @@ -34,6 +35,15 @@ public: }; }; +/** Flags to control the behavior of subs(). */ +class subs_options { +public: + enum { + subs_no_pattern = 0x0001, + subs_algebraic = 0x0002 + }; +}; + /** Flags to control series expansion. */ class series_options { public: diff --git a/ginac/idx.cpp b/ginac/idx.cpp index b60da944..4ba62bc2 100644 --- a/ginac/idx.cpp +++ b/ginac/idx.cpp @@ -351,7 +351,7 @@ ex idx::evalf(int level) const return *this; } -ex idx::subs(const lst & ls, const lst & lr, bool no_pattern) const +ex idx::subs(const lst & ls, const lst & lr, unsigned options) const { GINAC_ASSERT(ls.nops() == lr.nops()); @@ -372,7 +372,7 @@ ex idx::subs(const lst & ls, const lst & lr, bool no_pattern) const } // None, substitute objects in value (not in dimension) - const ex &subsed_value = value.subs(ls, lr, no_pattern); + const ex &subsed_value = value.subs(ls, lr, options); if (are_ex_trivially_equal(value, subsed_value)) return *this; diff --git a/ginac/idx.h b/ginac/idx.h index b05d555b..d498d32c 100644 --- a/ginac/idx.h +++ b/ginac/idx.h @@ -53,7 +53,7 @@ public: unsigned nops() const; ex & let_op(int i); ex evalf(int level = 0) const; - ex subs(const lst & ls, const lst & lr, bool no_pattern = false) const; + ex subs(const lst & ls, const lst & lr, unsigned options = 0) const; protected: ex derivative(const symbol & s) const; diff --git a/ginac/matrix.cpp b/ginac/matrix.cpp index 1b6d1b72..c240832d 100644 --- a/ginac/matrix.cpp +++ b/ginac/matrix.cpp @@ -235,14 +235,14 @@ ex matrix::eval(int level) const status_flags::evaluated); } -ex matrix::subs(const lst & ls, const lst & lr, bool no_pattern) const -{ +ex matrix::subs(const lst & ls, const lst & lr, unsigned options) const +{ exvector m2(row * col); for (unsigned r=0; r #include #include +#include #include "mul.h" #include "add.h" #include "power.h" #include "operators.h" #include "matrix.h" +#include "lst.h" #include "archive.h" #include "utils.h" @@ -516,6 +518,130 @@ ex mul::eval_ncmul(const exvector & v) const return inherited::eval_ncmul(v); } +bool tryfactsubs(const ex & origfactor, const ex & patternfactor, unsigned & nummatches, lst & repls) +{ + ex origbase; + int origexponent; + int origexpsign; + + if (is_exactly_a(origfactor) && origfactor.op(1).info(info_flags::integer)) { + origbase = origfactor.op(0); + int expon = ex_to(origfactor.op(1)).to_int(); + origexponent = expon > 0 ? expon : -expon; + origexpsign = expon > 0 ? 1 : -1; + } else { + origbase = origfactor; + origexponent = 1; + origexpsign = 1; + } + + ex patternbase; + int patternexponent; + int patternexpsign; + + if (is_exactly_a(patternfactor) && patternfactor.op(1).info(info_flags::integer)) { + patternbase = patternfactor.op(0); + int expon = ex_to(patternfactor.op(1)).to_int(); + patternexponent = expon > 0 ? expon : -expon; + patternexpsign = expon > 0 ? 1 : -1; + } else { + patternbase = patternfactor; + patternexponent = 1; + patternexpsign = 1; + } + + lst saverepls = repls; + if (origexponent < patternexponent || origexpsign != patternexpsign || !origbase.match(patternbase,saverepls)) + return false; + repls = saverepls; + + int newnummatches = origexponent / patternexponent; + if (newnummatches < nummatches) + nummatches = newnummatches; + return true; +} + +ex mul::algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) const +{ + std::vector subsed(seq.size(), false); + exvector subsresult(seq.size()); + + for (int i=0; i(ls.op(i))) { + + unsigned nummatches = std::numeric_limits::max(); + std::vector currsubsed(seq.size(), false); + bool succeed = true; + lst repls; + + for (int j=0; j::max(); + lst repls; + + for (int j=0; jnops(); j++) { + if (!subsed[j] && tryfactsubs(op(j), ls.op(i), nummatches, repls)) { + subsed[j] = true; + subsresult[j] = op(j) * power(lr.op(i).subs(ex(repls), subs_options::subs_no_pattern) / ls.op(i).subs(ex(repls), subs_options::subs_no_pattern), nummatches); + } + } + } + } + + bool subsfound = false; + for (int i=0; isetflag(status_flags::dynallocated); +} + // protected /** Implementation of ex::diff() for a product. It applies the product rule. diff --git a/ginac/mul.h b/ginac/mul.h index 4f2ae500..33cf2dd1 100644 --- a/ginac/mul.h +++ b/ginac/mul.h @@ -84,6 +84,8 @@ protected: // none // non-virtual functions in this class +public: + ex algebraic_subs_mul(const lst & ls, const lst & lr, unsigned options) const; protected: epvector * expandchildren(unsigned options) const; }; diff --git a/ginac/power.cpp b/ginac/power.cpp index ed26a4f7..43696a67 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -23,6 +23,7 @@ #include #include #include +#include #include "power.h" #include "expairseq.h" @@ -36,6 +37,7 @@ #include "matrix.h" #include "indexed.h" #include "symbol.h" +#include "lst.h" #include "print.h" #include "archive.h" #include "utils.h" @@ -523,16 +525,28 @@ ex power::evalm(void) const return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated); } -ex power::subs(const lst & ls, const lst & lr, bool no_pattern) const +extern bool tryfactsubs(const ex &, const ex &, unsigned &, lst &); + +ex power::subs(const lst & ls, const lst & lr, unsigned options) const { - const ex &subsed_basis = basis.subs(ls, lr, no_pattern); - const ex &subsed_exponent = exponent.subs(ls, lr, no_pattern); + if (options & subs_options::subs_algebraic) { + for (int i=0; i::max(); + lst repls; + if (tryfactsubs(*this, ls.op(i), nummatches, repls)) + return (ex_to((*this) * power(lr.op(i).subs(ex(repls), subs_options::subs_no_pattern) / ls.op(i).subs(ex(repls), subs_options::subs_no_pattern), nummatches))).basic::subs(ls, lr, options); + } + return basic::subs(ls, lr, options); + } + + const ex &subsed_basis = basis.subs(ls, lr, options); + const ex &subsed_exponent = exponent.subs(ls, lr, options); if (are_ex_trivially_equal(basis, subsed_basis) && are_ex_trivially_equal(exponent, subsed_exponent)) - return basic::subs(ls, lr, no_pattern); + return basic::subs(ls, lr, options); else - return power(subsed_basis, subsed_exponent).basic::subs(ls, lr, no_pattern); + return power(subsed_basis, subsed_exponent).basic::subs(ls, lr, options); } ex power::eval_ncmul(const exvector & v) const diff --git a/ginac/power.h b/ginac/power.h index e9ca134c..0386c60e 100644 --- a/ginac/power.h +++ b/ginac/power.h @@ -61,7 +61,7 @@ public: ex evalf(int level=0) const; ex evalm(void) const; ex series(const relational & s, int order, unsigned options = 0) const; - ex subs(const lst & ls, const lst & lr, bool no_pattern = false) 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 to_rational(lst &repl_lst) const; exvector get_free_indices(void) const; diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 3eedd8f8..9ac23966 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -410,13 +410,13 @@ ex pseries::evalf(int level) 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, bool no_pattern) const +ex pseries::subs(const lst & ls, const lst & lr, 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, no_pattern); + return convert_to_poly(true).subs(ls, lr, options); // Otherwise construct a new series with substituted coefficients and // expansion point @@ -424,10 +424,10 @@ ex pseries::subs(const lst & ls, const lst & lr, bool no_pattern) const 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, no_pattern), it->coeff)); + newseq.push_back(expair(it->rest.subs(ls, lr, options), it->coeff)); ++it; } - return (new pseries(relational(var,point.subs(ls, lr, no_pattern)), newseq))->setflag(status_flags::dynallocated); + return (new pseries(relational(var,point.subs(ls, lr, options)), newseq))->setflag(status_flags::dynallocated); } /** Implementation of ex::expand() for a power series. It expands all the diff --git a/ginac/pseries.h b/ginac/pseries.h index f2fc7591..d5c3ab99 100644 --- a/ginac/pseries.h +++ b/ginac/pseries.h @@ -54,7 +54,7 @@ public: 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, bool no_pattern = false) 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 expand(unsigned options = 0) const; protected: -- 2.47.0