improvement (hopefully) to subs and a patch included
Chris Dams
chrisd at sci.kun.nl
Mon Nov 18 13:11:28 CET 2002
Hello everybody,
Here is a patch that is supposed to make substitutions, optionally, a bit
more semantic (as opposed to syntactic). I also wrote a bit of
documentation for it, so hava a look at the patch for a description. (Also
corrected some spelling errors in the manual).
This disadvantage of this patch is that it breaks compatibility of the
interface. The subs method gets an extra optional argument like this
ex subs(const lst & ls, const lst & lr, unsigned options = 0, bool
no_pattern = false) const
This is inconvenient if people have written subclasses. They will need to
include this extra option everywhere. Furthermore if you use the
no_pattern optional argument, you will also hava a problem. However, I
think that the option-argument has the potential of being quite usefull,
and the no_pattern-argument is probably only for people getting really
deep into the internals of GiNaC, so I really think this is the correct
order. Note that further additions to the subs-method would go into new
subs_options, so that further breaking of the interface will not be
necesary. However, for this reason, I guess this patch should not be
included until the 1.1 release.
I do not know whether
diff -c -r --exclude=CVS
is the most convenient format for a patch for you. If this is going to
take you hours to do the patching, I will be happy to send some other
format.
Greetings,
Chris Dams
-------------- next part --------------
diff -c -r --exclude=CVS ginacorg/doc/tutorial/ginac.texi ginac/doc/tutorial/ginac.texi
*** ginacorg/doc/tutorial/ginac.texi Thu Oct 3 15:51:52 2002
--- ginac/doc/tutorial/ginac.texi Mon Nov 18 11:12:43 2002
***************
*** 734,740 ****
@cindex container
@cindex atom
! To get an idea about what kinds of symbolic composits may be built we
have a look at the most important classes in the class hierarchy and
some of the relations among the classes:
--- 734,740 ----
@cindex container
@cindex atom
! 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,901 ****
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
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
--- 895,901 ----
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 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,1015 ****
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
differ with regard to actually computed digits.
It should be clear that objects of class @code{numeric} should be used
--- 1009,1015 ----
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 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
***************
*** 1507,1513 ****
The @samp{algo} argument of @code{determinant()} allows to select
between different algorithms for calculating the determinant. The
! asymptotic speed (as parametrized by the matrix size) can greatly differ
between those algorithms, depending on the nature of the matrix'
entries. The possible values are defined in the @file{flags.h} header
file. By default, GiNaC uses a heuristic to automatically select an
--- 1507,1513 ----
The @samp{algo} argument of @code{determinant()} allows to select
between different algorithms for calculating the determinant. The
! asymptotic speed (as parameterized by the matrix size) can greatly differ
between those algorithms, depending on the nature of the matrix'
entries. The possible values are defined in the @file{flags.h} header
file. By default, GiNaC uses a heuristic to automatically select an
***************
*** 3040,3045 ****
--- 3040,3046 ----
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
***************
*** 3149,3154 ****
--- 3150,3156 ----
@{$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
***************
*** 3217,3222 ****
--- 3219,3225 ----
@{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
***************
*** 3259,3264 ****
--- 3262,3324 ----
@}
@end example
+ @subsection Even Smarter 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
***************
*** 5752,5758 ****
@end example
This @file{Makefile.am}, says that we are building a single executable,
! from a single sourcefile @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
--- 5812,5818 ----
@end example
This @file{Makefile.am}, says that we are building a single executable,
! 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 -c -r --exclude=CVS ginacorg/ginac/basic.cpp ginac/ginac/basic.cpp
*** ginacorg/ginac/basic.cpp Wed Jul 3 17:31:51 2002
--- ginac/ginac/basic.cpp Sat Nov 16 13:26:55 2002
***************
*** 503,509 ****
/** 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
{
GINAC_ASSERT(ls.nops() == lr.nops());
--- 503,509 ----
/** 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, bool no_pattern) const
{
GINAC_ASSERT(ls.nops() == lr.nops());
***************
*** 516,522 ****
for (unsigned i=0; i<ls.nops(); i++) {
lst repl_lst;
if (match(ex_to<basic>(ls.op(i)), repl_lst))
! return lr.op(i).subs(repl_lst, true); // avoid infinite recursion when re-substituting the wildcards
}
}
--- 516,522 ----
for (unsigned i=0; i<ls.nops(); i++) {
lst repl_lst;
if (match(ex_to<basic>(ls.op(i)), repl_lst))
! return lr.op(i).subs(repl_lst, options, true); // avoid infinite recursion when re-substituting the wildcards
}
}
***************
*** 684,693 ****
* 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
{
if (e.info(info_flags::relation_equal)) {
! return subs(lst(e), no_pattern);
}
if (!e.info(info_flags::list)) {
throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
--- 684,693 ----
* 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, bool no_pattern) const
{
if (e.info(info_flags::relation_equal)) {
! return subs(lst(e), options, no_pattern);
}
if (!e.info(info_flags::list)) {
throw(std::invalid_argument("basic::subs(ex): argument must be a list"));
***************
*** 702,708 ****
ls.append(r.op(0));
lr.append(r.op(1));
}
! return subs(ls, lr, no_pattern);
}
/** Compare objects syntactically to establish canonical ordering.
--- 702,708 ----
ls.append(r.op(0));
lr.append(r.op(1));
}
! return subs(ls, lr, options, no_pattern);
}
/** Compare objects syntactically to establish canonical ordering.
diff -c -r --exclude=CVS ginacorg/ginac/basic.h ginac/ginac/basic.h
*** ginacorg/ginac/basic.h Sun Jan 27 13:32:31 2002
--- ginac/ginac/basic.h Tue Nov 12 16:03:39 2002
***************
*** 112,118 ****
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 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;
--- 112,118 ----
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, unsigned options = 0, bool no_pattern = false) 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,141 ****
// non-virtual functions in this class
public:
! ex subs(const ex & e, bool no_pattern = false) const;
ex diff(const symbol & s, unsigned nth=1) const;
int compare(const basic & other) const;
bool is_equal(const basic & other) const;
--- 135,141 ----
// non-virtual functions in this class
public:
! ex subs(const ex & e, unsigned options = 0, bool no_pattern = false) const;
ex diff(const symbol & s, unsigned nth=1) const;
int compare(const basic & other) const;
bool is_equal(const basic & other) const;
diff -c -r --exclude=CVS ginacorg/ginac/container.pl ginac/ginac/container.pl
*** ginacorg/ginac/container.pl Mon Aug 5 18:03:02 2002
--- ginac/ginac/container.pl Thu Nov 14 22:31:16 2002
***************
*** 241,247 ****
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;
protected:
bool is_equal_same_type(const basic & other) const;
unsigned return_type(void) const;
--- 241,247 ----
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, unsigned options = 0, bool no_pattern = false) const;
protected:
bool is_equal_same_type(const basic & other) const;
unsigned return_type(void) const;
***************
*** 263,269 ****
protected:
bool is_canonical() const;
${STLT} evalchildren(int level) const;
! ${STLT} * subschildren(const lst & ls, const lst & lr, bool no_pattern = false) const;
protected:
${STLT} seq;
--- 263,269 ----
protected:
bool is_canonical() const;
${STLT} evalchildren(int level) const;
! ${STLT} * subschildren(const lst & ls, const lst & lr, unsigned options = 0, bool no_pattern = false) const;
protected:
${STLT} seq;
***************
*** 474,486 ****
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);
if (vp)
! return ex_to<basic>(this${CONTAINER}(vp)).basic::subs(ls, lr, no_pattern);
else
! return basic::subs(ls, lr, no_pattern);
}
// protected
--- 474,486 ----
return this${CONTAINER}(evalchildren(level));
}
! ex ${CONTAINER}::subs(const lst & ls, const lst & lr, unsigned options, bool no_pattern) const
! {
! ${STLT} *vp = subschildren(ls, lr, options, no_pattern);
if (vp)
! return ex_to<basic>(this${CONTAINER}(vp)).basic::subs(ls, lr, options, no_pattern);
else
! return basic::subs(ls, lr, options, no_pattern);
}
// protected
***************
*** 641,647 ****
return s;
}
! ${STLT} * ${CONTAINER}::subschildren(const lst & ls, const lst & lr, bool no_pattern) const
{
// returns a NULL pointer if nothing had to be substituted
// returns a pointer to a newly created epvector otherwise
--- 641,647 ----
return s;
}
! ${STLT} * ${CONTAINER}::subschildren(const lst & ls, const lst & lr, unsigned options, bool no_pattern) const
{
// returns a NULL pointer if nothing had to be substituted
// returns a pointer to a newly created epvector otherwise
***************
*** 649,655 ****
${STLT}::const_iterator cit = seq.begin(), end = seq.end();
while (cit != end) {
! const ex & subsed_ex = cit->subs(ls, lr, no_pattern);
if (!are_ex_trivially_equal(*cit, subsed_ex)) {
// something changed, copy seq, subs and return it
--- 649,655 ----
${STLT}::const_iterator cit = seq.begin(), end = seq.end();
while (cit != end) {
! const ex & subsed_ex = cit->subs(ls, lr, options, no_pattern);
if (!are_ex_trivially_equal(*cit, subsed_ex)) {
// something changed, copy seq, subs and return it
***************
*** 669,675 ****
// copy rest
while (cit2 != end) {
! s->push_back(cit2->subs(ls, lr, no_pattern));
++cit2;
}
return s;
--- 669,675 ----
// copy rest
while (cit2 != end) {
! s->push_back(cit2->subs(ls, lr, options, no_pattern));
++cit2;
}
return s;
diff -c -r --exclude=CVS ginacorg/ginac/ex.h ginac/ginac/ex.h
*** ginacorg/ginac/ex.h Thu Jul 18 16:08:52 2002
--- ginac/ginac/ex.h Tue Nov 12 16:05:25 2002
***************
*** 146,153 ****
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); }
exvector get_free_indices(void) const { return bp->get_free_indices(); }
ex simplify_indexed(void) const;
ex simplify_indexed(const scalar_products & sp) const;
--- 146,153 ----
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, unsigned options = 0, bool no_pattern = false) const { return bp->subs(ls, lr, options, no_pattern); }
! ex subs(const ex & e, unsigned options = 0, bool no_pattern = false) const { return bp->subs(e, options, no_pattern); }
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,435 ****
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 lst & ls, const lst & lr)
! { return thisex.subs(ls, lr); }
inline ex simplify_indexed(const ex & thisex)
{ return thisex.simplify_indexed(); }
--- 425,435 ----
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(); }
diff -c -r --exclude=CVS ginacorg/ginac/expairseq.cpp ginac/ginac/expairseq.cpp
*** ginacorg/ginac/expairseq.cpp Wed Jul 3 17:31:51 2002
--- ginac/ginac/expairseq.cpp Fri Nov 15 12:44:10 2002
***************
*** 401,413 ****
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);
if (vp)
return ex_to<basic>(thisexpairseq(vp, overall_coeff));
else
! return basic::subs(ls, lr, no_pattern);
}
// protected
--- 401,417 ----
return inherited::match(pattern, repl_lst);
}
! ex expairseq::subs(const lst &ls, const lst &lr, unsigned options, bool no_pattern) const
! {
! if ((options & subs_options::subs_smart) && is_a<mul>(*this))
! {
! return ((mul*)this) -> smartsubsmul(ls, lr, no_pattern);
! }
! epvector *vp = subschildren(ls, lr, options, no_pattern);
if (vp)
return ex_to<basic>(thisexpairseq(vp, overall_coeff));
else
! return basic::subs(ls, lr, options, no_pattern);
}
// protected
***************
*** 1559,1565 ****
* @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
{
GINAC_ASSERT(ls.nops()==lr.nops());
--- 1563,1569 ----
* @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, bool no_pattern) const
{
GINAC_ASSERT(ls.nops()==lr.nops());
***************
*** 1573,1578 ****
--- 1577,1583 ----
break;
}
+
if (complex_subs) {
// Substitute in the recombined pairs
***************
*** 1580,1586 ****
while (cit != last) {
const ex &orig_ex = recombine_pair_to_ex(*cit);
! const ex &subsed_ex = orig_ex.subs(ls, lr, no_pattern);
if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
// Something changed, copy seq, subs and return it
--- 1585,1591 ----
while (cit != last) {
const ex &orig_ex = recombine_pair_to_ex(*cit);
! const ex &subsed_ex = orig_ex.subs(ls, lr, options, no_pattern);
if (!are_ex_trivially_equal(orig_ex, subsed_ex)) {
// Something changed, copy seq, subs and return it
***************
*** 1596,1602 ****
// Copy rest
while (cit != last) {
! s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(ls, lr, no_pattern)));
++cit;
}
return s;
--- 1601,1607 ----
// Copy rest
while (cit != last) {
! s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit).subs(ls, lr, options, no_pattern)));
++cit;
}
return s;
***************
*** 1611,1617 ****
epvector::const_iterator cit = seq.begin(), last = seq.end();
while (cit != last) {
! const ex &subsed_ex = cit->rest.subs(ls, lr, no_pattern);
if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
// Something changed, copy seq, subs and return it
--- 1616,1622 ----
epvector::const_iterator cit = seq.begin(), last = seq.end();
while (cit != last) {
! const ex &subsed_ex = cit->rest.subs(ls, lr, options, no_pattern);
if (!are_ex_trivially_equal(cit->rest, subsed_ex)) {
// Something changed, copy seq, subs and return it
***************
*** 1627,1633 ****
// Copy rest
while (cit != last) {
! s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(ls, lr, no_pattern),
cit->coeff));
++cit;
}
--- 1632,1638 ----
// Copy rest
while (cit != last) {
! s->push_back(combine_ex_with_coeff_to_pair(cit->rest.subs(ls, lr, options, no_pattern),
cit->coeff));
++cit;
}
diff -c -r --exclude=CVS ginacorg/ginac/expairseq.h ginac/ginac/expairseq.h
*** ginacorg/ginac/expairseq.h Sun Jan 27 13:32:31 2002
--- ginac/ginac/expairseq.h Tue Nov 12 23:03:22 2002
***************
*** 96,102 ****
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;
protected:
int compare_same_type(const basic & other) const;
bool is_equal_same_type(const basic & other) const;
--- 96,102 ----
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, unsigned options = 0, bool no_pattern = false) const;
protected:
int compare_same_type(const basic & other) const;
bool is_equal_same_type(const basic & other) const;
***************
*** 163,169 ****
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;
// member variables
--- 163,169 ----
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, bool no_pattern = false) const;
// member variables
diff -c -r --exclude=CVS ginacorg/ginac/flags.h ginac/ginac/flags.h
*** ginacorg/ginac/flags.h Thu Oct 3 15:51:38 2002
--- ginac/ginac/flags.h Tue Nov 12 16:25:44 2002
***************
*** 34,39 ****
--- 34,46 ----
};
};
+ class subs_options {
+ public:
+ enum {
+ subs_smart = 0x0001
+ };
+ };
+
/** Flags to control series expansion. */
class series_options {
public:
diff -c -r --exclude=CVS ginacorg/ginac/idx.cpp ginac/ginac/idx.cpp
*** ginacorg/ginac/idx.cpp Wed Oct 9 16:12:16 2002
--- ginac/ginac/idx.cpp Tue Nov 12 16:16:08 2002
***************
*** 351,357 ****
return *this;
}
! ex idx::subs(const lst & ls, const lst & lr, bool no_pattern) const
{
GINAC_ASSERT(ls.nops() == lr.nops());
--- 351,357 ----
return *this;
}
! ex idx::subs(const lst & ls, const lst & lr, unsigned options, bool no_pattern) const
{
GINAC_ASSERT(ls.nops() == lr.nops());
***************
*** 372,378 ****
}
// None, substitute objects in value (not in dimension)
! const ex &subsed_value = value.subs(ls, lr, no_pattern);
if (are_ex_trivially_equal(value, subsed_value))
return *this;
--- 372,378 ----
}
// None, substitute objects in value (not in dimension)
! const ex &subsed_value = value.subs(ls, lr, options, no_pattern);
if (are_ex_trivially_equal(value, subsed_value))
return *this;
diff -c -r --exclude=CVS ginacorg/ginac/idx.h ginac/ginac/idx.h
*** ginacorg/ginac/idx.h Wed Oct 9 16:12:16 2002
--- ginac/ginac/idx.h Tue Nov 12 16:10:01 2002
***************
*** 53,59 ****
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;
protected:
ex derivative(const symbol & s) const;
--- 53,59 ----
unsigned nops() const;
ex & let_op(int i);
ex evalf(int level = 0) const;
! ex subs(const lst & ls, const lst & lr, unsigned options = 0, bool no_pattern = false) const;
protected:
ex derivative(const symbol & s) const;
diff -c -r --exclude=CVS ginacorg/ginac/matrix.cpp ginac/ginac/matrix.cpp
*** ginacorg/ginac/matrix.cpp Wed Jul 3 17:31:51 2002
--- ginac/ginac/matrix.cpp Fri Nov 15 12:45:25 2002
***************
*** 233,246 ****
status_flags::evaluated );
}
! ex matrix::subs(const lst & ls, const lst & lr, bool no_pattern) 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, no_pattern);
! return matrix(row, col, m2).basic::subs(ls, lr, no_pattern);
}
// protected
--- 233,246 ----
status_flags::evaluated );
}
! ex matrix::subs(const lst & ls, const lst & lr, unsigned options, bool no_pattern) 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, no_pattern);
! return matrix(row, col, m2).basic::subs(ls, lr, options, no_pattern);
}
// protected
diff -c -r --exclude=CVS ginacorg/ginac/matrix.h ginac/ginac/matrix.h
*** ginacorg/ginac/matrix.h Fri Mar 8 17:39:57 2002
--- ginac/ginac/matrix.h Tue Nov 12 22:01:07 2002
***************
*** 48,54 ****
ex & let_op(int i);
ex eval(int level=0) const;
ex evalm(void) const {return *this;}
! ex subs(const lst & ls, const lst & lr, bool no_pattern = false) 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;
--- 48,54 ----
ex & let_op(int i);
ex eval(int level=0) const;
ex evalm(void) const {return *this;}
! ex subs(const lst & ls, const lst & lr, unsigned options = 0, bool no_pattern = false) 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;
diff -c -r --exclude=CVS ginacorg/ginac/mul.cpp ginac/ginac/mul.cpp
*** ginacorg/ginac/mul.cpp Mon Sep 2 16:43:08 2002
--- ginac/ginac/mul.cpp Sat Nov 16 13:40:33 2002
***************
*** 31,36 ****
--- 31,37 ----
#include "matrix.h"
#include "archive.h"
#include "utils.h"
+ #include "lst.h"
namespace GiNaC {
***************
*** 484,489 ****
--- 485,595 ----
return inherited::simplify_ncmul(v);
}
+ bool tryfactsubs(const ex&origfactor, const ex&patternfactor, unsigned &nummatches, lst&repls)
+ { ex origbase;
+ int origexponent;
+ int origexpsign;
+ if( is_a<power>(origfactor) && origfactor.op(1).info(info_flags::integer) ) {
+ origbase = origfactor.op(0);
+ int expon = ex_to<numeric>(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_a<power>(patternfactor) && patternfactor.op(1).info(info_flags::integer) ) {
+ patternbase = patternfactor.op(0);
+ int expon = ex_to<numeric>(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::smartsubsmul(const lst &ls, const lst &lr, bool no_pattern) const
+ { std::vector<bool> subsed( seq.size(), false );
+ exvector subsresult(seq.size());
+ for( int i=0; i<ls.nops(); i++ ) {
+ if(is_a<mul>(ls.op(i))) {
+ unsigned nummatches=0xffff;
+ std::vector<bool>currsubsed ( seq.size(), false );
+ bool succeed=true;
+ lst repls;
+ for( int j=0; j<ls.op(i).nops(); j++ ) {
+ bool found=false;
+ for( int k=0; k<this->nops(); k++ ) {
+ if( currsubsed[k] || subsed[k] )
+ continue;
+ if( tryfactsubs(this->op(k), ls.op(i).op(j), nummatches, repls)) {
+ currsubsed[k] = true;
+ found = true;
+ break;
+ }
+ }
+ if(!found) {
+ succeed = false;
+ break;
+ }
+ }
+ if(!succeed)
+ continue;
+ bool foundfirstsubsedfactor = false;
+ for(int j=0; j<subsed.size(); j++) {
+ if(currsubsed[j]) {
+ if(foundfirstsubsedfactor)
+ subsresult[j] = this->op(j);
+ else {
+ foundfirstsubsedfactor = true;
+ subsresult[j] = this->op(j)*power( lr.op(i).subs(ex(repls), 0, true) / ls.op(i).subs(ex(repls), 0, true), nummatches);
+ }
+ subsed[j] = true;
+ }
+ }
+ } else {
+ unsigned nummatches = 0xffff;
+ lst repls;
+ for(int j=0; j<this->nops(); j++) {
+ if( !subsed[j] && tryfactsubs( this->op(j), ls.op(i), nummatches, repls )) {
+ subsed[j] = true;
+ subsresult[j] = this->op(j)*power(lr.op(i).subs(ex(repls), 0, true)/ls.op(i).subs(ex(repls), 0, true), nummatches);
+ }
+ }
+ }
+ }
+ bool subsfound = false;
+ for(int i=0; i<subsed.size(); i++) {
+ if(subsed[i]) {
+ subsfound = true;
+ break;
+ }
+ }
+ if(!subsfound)
+ return basic::subs(ls, lr, subs_options::subs_smart, no_pattern);
+ exvector ev( this->nops() );
+ for(int i=0; i<this->nops(); i++)
+ if(subsed[i])
+ ev[i] = subsresult[i];
+ else
+ ev[i] = this->op(i);
+ return ex( (new mul(ev)) -> setflag(status_flags::dynallocated) );
+ }
+
// protected
/** Implementation of ex::diff() for a product. It applies the product rule.
diff -c -r --exclude=CVS ginacorg/ginac/mul.h ginac/ginac/mul.h
*** ginacorg/ginac/mul.h Tue Apr 30 19:42:41 2002
--- ginac/ginac/mul.h Thu Nov 14 12:37:24 2002
***************
*** 63,68 ****
--- 63,69 ----
numeric max_coefficient(void) const;
exvector get_free_indices(void) const;
ex simplify_ncmul(const exvector & v) const;
+ ex smartsubsmul(const lst &ls, const lst&lr, bool no_pattern = false) const;
protected:
ex derivative(const symbol & s) const;
unsigned return_type(void) const;
diff -c -r --exclude=CVS ginacorg/ginac/power.cpp ginac/ginac/power.cpp
*** ginacorg/ginac/power.cpp Thu Oct 24 17:06:44 2002
--- ginac/ginac/power.cpp Sat Nov 16 13:40:35 2002
***************
*** 39,44 ****
--- 39,45 ----
#include "print.h"
#include "archive.h"
#include "utils.h"
+ #include "lst.h"
namespace GiNaC {
***************
*** 511,526 ****
return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated);
}
! ex power::subs(const lst & ls, const lst & lr, bool no_pattern) const
{
! const ex &subsed_basis = basis.subs(ls, lr, no_pattern);
! const ex &subsed_exponent = exponent.subs(ls, lr, no_pattern);
if (are_ex_trivially_equal(basis, subsed_basis)
&& are_ex_trivially_equal(exponent, subsed_exponent))
! return basic::subs(ls, lr, no_pattern);
else
! return power(subsed_basis, subsed_exponent).basic::subs(ls, lr, no_pattern);
}
ex power::simplify_ncmul(const exvector & v) const
--- 512,539 ----
return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated);
}
! extern bool tryfactsubs(const ex&, const ex&, unsigned&, lst&);
!
! ex power::subs(const lst & ls, const lst & lr, unsigned options, bool no_pattern) const
{
! if( options && subs_options::subs_smart ) {
! for(int i=0; i<ls.nops(); i++ ) {
! unsigned nummatches=0xffff;
! lst repls;
! if(tryfactsubs(*this, ls.op(i), nummatches, repls))
! return (ex_to<basic>((*this)*power(lr.op(i).subs(ex(repls), 0, true)/ls.op(i).subs(ex(repls), 0, true),nummatches))).basic::subs(ls, lr, options, no_pattern);
! }
! return basic::subs(ls, lr, options, no_pattern);
! }
!
! const ex &subsed_basis = basis.subs(ls, lr, options, no_pattern);
! const ex &subsed_exponent = exponent.subs(ls, lr, options, no_pattern);
if (are_ex_trivially_equal(basis, subsed_basis)
&& are_ex_trivially_equal(exponent, subsed_exponent))
! return basic::subs(ls, lr, options, no_pattern);
else
! return power(subsed_basis, subsed_exponent).basic::subs(ls, lr, options, no_pattern);
}
ex power::simplify_ncmul(const exvector & v) const
diff -c -r --exclude=CVS ginacorg/ginac/power.h ginac/ginac/power.h
*** ginacorg/ginac/power.h Sun Jan 27 13:32:32 2002
--- ginac/ginac/power.h Tue Nov 12 16:10:49 2002
***************
*** 61,67 ****
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 normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
ex to_rational(lst &repl_lst) const;
exvector get_free_indices(void) const;
--- 61,67 ----
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, unsigned options = 0, bool no_pattern = false) 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 -c -r --exclude=CVS ginacorg/ginac/pseries.cpp ginac/ginac/pseries.cpp
*** ginacorg/ginac/pseries.cpp Wed Aug 7 17:47:49 2002
--- ginac/ginac/pseries.cpp Tue Nov 12 16:20:37 2002
***************
*** 410,422 ****
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
{
// 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);
// Otherwise construct a new series with substituted coefficients and
// expansion point
--- 410,422 ----
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, bool no_pattern) 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, no_pattern);
// Otherwise construct a new series with substituted coefficients and
// expansion point
***************
*** 424,433 ****
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));
++it;
}
! return (new pseries(relational(var,point.subs(ls, lr, no_pattern)), newseq))->setflag(status_flags::dynallocated);
}
/** Implementation of ex::expand() for a power series. It expands all the
--- 424,433 ----
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, no_pattern), it->coeff));
++it;
}
! return (new pseries(relational(var,point.subs(ls, lr, options, no_pattern)), newseq))->setflag(status_flags::dynallocated);
}
/** Implementation of ex::expand() for a power series. It expands all the
diff -c -r --exclude=CVS ginacorg/ginac/pseries.h ginac/ginac/pseries.h
*** ginacorg/ginac/pseries.h Sun Jan 27 13:32:32 2002
--- ginac/ginac/pseries.h Tue Nov 12 16:11:08 2002
***************
*** 54,60 ****
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 normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
ex expand(unsigned options = 0) const;
protected:
--- 54,60 ----
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, bool no_pattern = false) const;
ex normal(lst &sym_lst, lst &repl_lst, int level = 0) const;
ex expand(unsigned options = 0) const;
protected:
More information about the GiNaC-devel
mailing list