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: