Index: ginac/integral.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/integral.h,v
retrieving revision 1.1
diff -c -r1.1 integral.h
*** ginac/integral.h	12 Oct 2004 09:32:11 -0000	1.1
--- ginac/integral.h	27 Oct 2004 07:57:04 -0000
***************
*** 56,61 ****
--- 56,62 ----
  	ex eval_integ() const;
  protected:
  	ex derivative(const symbol & s) const;
+ 	ex series(const relational & r, int order, unsigned options = 0) const;
  
  	// new virtual functions which can be overridden by derived classes
  	// none
Index: ginac/pseries.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/pseries.cpp,v
retrieving revision 1.82
diff -c -r1.82 pseries.cpp
*** ginac/pseries.cpp	6 Oct 2004 11:59:22 -0000	1.82
--- ginac/pseries.cpp	27 Oct 2004 07:57:04 -0000
***************
*** 33,38 ****
--- 33,39 ----
  #include "relational.h"
  #include "operators.h"
  #include "symbol.h"
+ #include "integral.h"
  #include "archive.h"
  #include "utils.h"
  
***************
*** 266,271 ****
--- 267,274 ----
  	if (i >= seq.size())
  		throw (std::out_of_range("op() out of range"));
  
+ 	if(is_order_function(seq[i].rest))
+ 		return Order(power(var-point, seq[i].coeff));
  	return seq[i].rest * power(var - point, seq[i].coeff);
  }
  
***************
*** 424,429 ****
--- 427,457 ----
  	return result;
  }
  
+ ex pseries::eval_integ() const
+ {
+ 	epvector *newseq = 0;
+ 	for(epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+ 		if(newseq)
+ 		{	newseq->push_back(expair(i->rest.eval_integ(), i->coeff));
+ 			continue;
+ 		}
+ 		ex newterm = i->rest.eval_integ();
+ 		if(!are_ex_trivially_equal(newterm, i->rest)) {
+ 			newseq = new epvector;
+ 			newseq->reserve(seq.size());
+ 			for(epvector::const_iterator j=seq.begin(); j!=i; ++j)
+ 				newseq->push_back(*j);
+ 			newseq->push_back(expair(newterm, i->coeff));
+ 		}
+ 	}
+ 
+ 	ex newpoint = point.eval_integ();
+ 	if(newseq || !are_ex_trivially_equal(newpoint, point))
+ 		return (new pseries(var==newpoint, *newseq))
+ 			->setflag(status_flags::dynallocated);
+ 	return *this;
+ }
+ 
  ex pseries::subs(const exmap & m, unsigned options) const
  {
  	// If expansion variable is being substituted, convert the series to a
***************
*** 519,524 ****
--- 547,566 ----
  	return seq.empty() || !is_order_function((seq.end()-1)->rest);
  }
  
+ ex pseries::coeffop(size_t i) const
+ {
+ 	if(i > nops())
+ 		throw (std::out_of_range("coeffop() out of range"));
+ 	return seq[i].rest;
+ }
+ 
+ ex pseries::exponop(size_t i) const
+ {
+ 	if(i > nops())
+ 		throw (std::out_of_range("exponop() out of range"));
+ 	return seq[i].coeff;
+ }
+ 
  
  /*
   *  Implementations of series expansion
***************
*** 729,734 ****
--- 771,781 ----
  		nul.push_back(expair(Order(_ex1), _ex0));
  		return pseries(relational(var,point), nul);
  	}
+ 
+ 	if(seq.empty() || other.seq.empty()) {
+ 		return (new pseries(var==point, epvector()))
+ 			->setflag(status_flags::dynallocated);
+ 	}
  	
  	// Series multiplication
  	epvector new_seq;
***************
*** 880,886 ****
  		throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
  
  	// adjust number of coefficients
! 	deg = deg - (p*ldeg).to_int();
  	
  	// O(x^n)^(-m) is undefined
  	if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
--- 927,940 ----
  		throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
  
  	// adjust number of coefficients
! 	int numcoeff = deg - (p*ldeg).to_int();
! 	if(numcoeff <= 0) {
! 		epvector epv;
! 		epv.reserve(1);
! 		epv.push_back(expair(Order(_ex1), deg));
! 		return (new pseries(relational(var,point), epv))
! 			->setflag(status_flags::dynallocated);
! 	}
  	
  	// O(x^n)^(-m) is undefined
  	if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
***************
*** 888,896 ****
  	
  	// Compute coefficients of the powered series
  	exvector co;
! 	co.reserve(deg);
  	co.push_back(power(coeff(var, ldeg), p));
! 	for (int i=1; i<deg; ++i) {
  		ex sum = _ex0;
  		for (int j=1; j<=i; ++j) {
  			ex c = coeff(var, j + ldeg);
--- 942,950 ----
  	
  	// Compute coefficients of the powered series
  	exvector co;
! 	co.reserve(numcoeff);
  	co.push_back(power(coeff(var, ldeg), p));
! 	for (int i=1; i<numcoeff; ++i) {
  		ex sum = _ex0;
  		for (int j=1; j<=i; ++j) {
  			ex c = coeff(var, j + ldeg);
***************
*** 906,912 ****
  	// Construct new series (of non-zero coefficients)
  	epvector new_seq;
  	bool higher_order = false;
! 	for (int i=0; i<deg; ++i) {
  		if (!co[i].is_zero())
  			new_seq.push_back(expair(co[i], p * ldeg + i));
  		if (is_order_function(co[i])) {
--- 960,966 ----
  	// Construct new series (of non-zero coefficients)
  	epvector new_seq;
  	bool higher_order = false;
! 	for (int i=0; i<numcoeff; ++i) {
  		if (!co[i].is_zero())
  			new_seq.push_back(expair(co[i], p * ldeg + i));
  		if (is_order_function(co[i])) {
***************
*** 915,921 ****
  		}
  	}
  	if (!higher_order)
! 		new_seq.push_back(expair(Order(_ex1), p * ldeg + deg));
  
  	return pseries(relational(var,point), new_seq);
  }
--- 969,975 ----
  		}
  	}
  	if (!higher_order)
! 		new_seq.push_back(expair(Order(_ex1), p * ldeg + numcoeff));
  
  	return pseries(relational(var,point), new_seq);
  }
***************
*** 1033,1038 ****
--- 1087,1149 ----
  		}
  	} else
  		return convert_to_poly().series(r, order, options);
+ }
+ 
+ ex integral::series(const relational & r, int order, unsigned options) const
+ {
+ 	if(x.subs(r)!=x)
+ 		throw std::logic_error("Cannot series expand wrt dummy variable");
+ 	
+ 	// Expanding integrant with r substituted taken in boundaries.
+ 	ex fseries = f.series(r, order, options);
+ 	epvector fexpansion;
+ 	fexpansion.reserve(fseries.nops());
+ 	for(size_t i=0; i<fseries.nops(); ++i) {
+ 		ex currcoeff=ex_to<pseries>(fseries).coeffop(i);
+ 		fexpansion.push_back(expair(
+ 			currcoeff==Order(_ex1)
+ 				? currcoeff
+ 				: integral(x, a.subs(r), b.subs(r), currcoeff),
+ 			ex_to<pseries>(fseries).exponop(i)
+ 		));
+ 	}
+ 
+ 	// Expanding lower boundary
+ 	ex result=(new pseries(r, fexpansion))->setflag(status_flags::dynallocated);
+ 	ex aseries = (a-a.subs(r)).series(r, order, options);
+ 	fseries = f.series(x==(a.subs(r)), order, options);
+ 	for(size_t i=0; i<fseries.nops(); ++i) {
+ 		ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ 		if(is_order_function(currcoeff))
+ 			break;
+ 		ex currexpon = ex_to<pseries>(fseries).exponop(i);
+ 		int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
+ 		currcoeff=currcoeff.series(r, orderforf);
+ 		ex term = ex_to<pseries>(aseries)
+ 			.power_const(ex_to<numeric>(currexpon+1),order);
+ 		term=ex_to<pseries>(term).mul_const(ex_to<numeric>(-1/(currexpon+1)));
+ 		term=ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
+ 		result=ex_to<pseries>(result).add_series(ex_to<pseries>(term));
+ 	}
+ 
+ 	// Expanding upper boundary
+ 	ex bseries = (b-b.subs(r)).series(r, order, options);
+ 	fseries = f.series(x==(b.subs(r)), order, options);
+ 	for(size_t i=0; i<fseries.nops(); ++i) {
+ 		ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ 		if(is_order_function(currcoeff))
+ 			break;
+ 		ex currexpon = ex_to<pseries>(fseries).exponop(i);
+ 		int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
+ 		currcoeff=currcoeff.series(r, orderforf);
+ 		ex term = ex_to<pseries>(bseries)
+ 			.power_const(ex_to<numeric>(currexpon+1),order);
+ 		term=ex_to<pseries>(term).mul_const(ex_to<numeric>(1/(currexpon+1)));
+ 		term=ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
+ 		result=ex_to<pseries>(result).add_series(ex_to<pseries>(term));
+ 	}
+ 
+ 	return result;
  }
  
  
Index: ginac/pseries.h
===================================================================
RCS file: /home/cvs/GiNaC/ginac/pseries.h,v
retrieving revision 1.36
diff -c -r1.36 pseries.h
*** ginac/pseries.h	4 Jan 2004 13:11:45 -0000	1.36
--- ginac/pseries.h	27 Oct 2004 07:57:04 -0000
***************
*** 56,61 ****
--- 56,62 ----
  	ex normal(exmap & repl, exmap & rev_lookup, int level = 0) const;
  	ex expand(unsigned options = 0) const;
  	ex conjugate() const;
+ 	ex eval_integ() const;
  protected:
  	ex derivative(const symbol & s) const;
  
***************
*** 82,87 ****
--- 83,92 ----
  	/** Returns true if there is no order term, i.e. the series terminates and
  	 *  false otherwise. */
  	bool is_terminating() const;
+ 
+ 	/** Get coefficients and exponents. */
+ 	ex coeffop(size_t i) const;
+ 	ex exponop(size_t i) const;
  
  	ex add_series(const pseries &other) const;
  	ex mul_const(const numeric &other) const;