X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpseries.cpp;h=c290fe0a1f7b09b0c73fa2cebdf4e9cb1f609c38;hp=bd8e4803a0feb57e9778ceda4cfe4d8e3cd3c56a;hb=1d73fed17debbc88819dc7f84a4a2e608d8eb978;hpb=fc6d6774bb3d125be622bd84edb7e2ed4f77655c

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index bd8e4803..c290fe0a 100644
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -4,7 +4,7 @@
  *  methods for series expansion. */
 
 /*
- *  GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2010 Johannes Gutenberg University Mainz, Germany
  *
  *  This program is free software; you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -21,9 +21,6 @@
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#include <numeric>
-#include <stdexcept>
-
 #include "pseries.h"
 #include "add.h"
 #include "inifcns.h" // for Order function
@@ -37,6 +34,10 @@
 #include "archive.h"
 #include "utils.h"
 
+#include <limits>
+#include <numeric>
+#include <stdexcept>
+
 namespace GiNaC {
 
 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(pseries, basic,
@@ -51,7 +52,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(pseries, basic,
  *  Default constructor
  */
 
-pseries::pseries() : inherited(&pseries::tinfo_static) { }
+pseries::pseries() { }
 
 
 /*
@@ -67,7 +68,7 @@ pseries::pseries() : inherited(&pseries::tinfo_static) { }
  *  @param rel_  expansion variable and point (must hold a relational)
  *  @param ops_  vector of {coefficient, power} pairs (coefficient must not be zero)
  *  @return newly constructed pseries */
-pseries::pseries(const ex &rel_, const epvector &ops_) : basic(&pseries::tinfo_static), seq(ops_)
+pseries::pseries(const ex &rel_, const epvector &ops_) : seq(ops_)
 {
 	GINAC_ASSERT(is_a<relational>(rel_));
 	GINAC_ASSERT(is_a<symbol>(rel_.lhs()));
@@ -80,8 +81,9 @@ pseries::pseries(const ex &rel_, const epvector &ops_) : basic(&pseries::tinfo_s
  *  Archiving
  */
 
-pseries::pseries(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
+void pseries::read_archive(const archive_node &n, lst &sym_lst) 
 {
+	inherited::read_archive(n, sym_lst);
 	archive_node::archive_node_cit first = n.find_first("coeff");
 	archive_node::archive_node_cit last = n.find_last("power");
 	++last;
@@ -112,7 +114,6 @@ void pseries::archive(archive_node &n) const
 	n.add_ex("point", point);
 }
 
-DEFAULT_UNARCHIVE(pseries)
 
 //////////
 // functions overriding virtual functions from base classes
@@ -292,7 +293,7 @@ int pseries::degree(const ex &s) const
 		epvector::const_iterator it = seq.begin(), itend = seq.end();
 		if (it == itend)
 			return 0;
-		int max_pow = INT_MIN;
+		int max_pow = std::numeric_limits<int>::min();
 		while (it != itend) {
 			int pow = it->rest.degree(s);
 			if (pow > max_pow)
@@ -320,7 +321,7 @@ int pseries::ldegree(const ex &s) const
 		epvector::const_iterator it = seq.begin(), itend = seq.end();
 		if (it == itend)
 			return 0;
-		int min_pow = INT_MAX;
+		int min_pow = std::numeric_limits<int>::max();
 		while (it != itend) {
 			int pow = it->rest.ldegree(s);
 			if (pow < min_pow)
@@ -489,6 +490,34 @@ ex pseries::eval_integ() const
 	return *this;
 }
 
+ex pseries::evalm() const
+{
+	// evalm each coefficient
+	epvector newseq;
+	bool something_changed = false;
+	for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+		if (something_changed) {
+			ex newcoeff = i->rest.evalm();
+			if (!newcoeff.is_zero())
+				newseq.push_back(expair(newcoeff, i->coeff));
+		}
+		else {
+			ex newcoeff = i->rest.evalm();
+			if (!are_ex_trivially_equal(newcoeff, i->rest)) {
+				something_changed = true;
+				newseq.reserve(seq.size());
+				std::copy(seq.begin(), i, std::back_inserter<epvector>(newseq));
+				if (!newcoeff.is_zero())
+					newseq.push_back(expair(newcoeff, i->coeff));
+			}
+		}
+	}
+	if (something_changed)
+		return (new pseries(var==point, newseq))->setflag(status_flags::dynallocated);
+	else
+		return *this;
+}
+
 ex pseries::subs(const exmap & m, unsigned options) const
 {
 	// If expansion variable is being substituted, convert the series to a
@@ -690,7 +719,7 @@ ex pseries::add_series(const pseries &other) const
 	epvector::const_iterator b = other.seq.begin();
 	epvector::const_iterator a_end = seq.end();
 	epvector::const_iterator b_end = other.seq.end();
-	int pow_a = INT_MAX, pow_b = INT_MAX;
+	int pow_a = std::numeric_limits<int>::max(), pow_b = std::numeric_limits<int>::max();
 	for (;;) {
 		// If a is empty, fill up with elements from b and stop
 		if (a == a_end) {
@@ -823,8 +852,8 @@ ex pseries::mul_series(const pseries &other) const
 	int cdeg_min = a_min + b_min;
 	int cdeg_max = a_max + b_max;
 	
-	int higher_order_a = INT_MAX;
-	int higher_order_b = INT_MAX;
+	int higher_order_a = std::numeric_limits<int>::max();
+	int higher_order_b = std::numeric_limits<int>::max();
 	if (is_order_function(coeff(var, a_max)))
 		higher_order_a = a_max + b_min;
 	if (is_order_function(other.coeff(var, b_max)))
@@ -845,7 +874,7 @@ ex pseries::mul_series(const pseries &other) const
 		if (!co.is_zero())
 			new_seq.push_back(expair(co, numeric(cdeg)));
 	}
-	if (higher_order_c < INT_MAX)
+	if (higher_order_c < std::numeric_limits<int>::max())
 		new_seq.push_back(expair(Order(_ex1), numeric(higher_order_c)));
 	return pseries(relational(var, point), new_seq);
 }
@@ -863,10 +892,13 @@ ex mul::series(const relational & r, int order, unsigned options) const
 		
 	// holds ldegrees of the series of individual factors
 	std::vector<int> ldegrees;
+	std::vector<bool> ldegree_redo;
 
 	// find minimal degrees
 	const epvector::const_iterator itbeg = seq.begin();
 	const epvector::const_iterator itend = seq.end();
+	// first round: obtain a bound up to which minimal degrees have to be
+	// considered
 	for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
 
 		ex expon = it->coeff;
@@ -880,19 +912,61 @@ ex mul::series(const relational & r, int order, unsigned options) const
 		}
 
 		int real_ldegree = 0;
+		bool flag_redo = false;
 		try {
 			real_ldegree = buf.expand().ldegree(sym-r.rhs());
 		} catch (std::runtime_error) {}
 
 		if (real_ldegree == 0) {
+			if ( factor < 0 ) {
+				// This case must terminate, otherwise we would have division by
+				// zero.
+				int orderloop = 0;
+				do {
+					orderloop++;
+					real_ldegree = buf.series(r, orderloop, options).ldegree(sym);
+				} while (real_ldegree == orderloop);
+			} else {
+				// Here it is possible that buf does not have a ldegree, therefore
+				// check only if ldegree is negative, otherwise reconsider the case
+				// in the second round.
+				real_ldegree = buf.series(r, 0, options).ldegree(sym);
+				if (real_ldegree == 0)
+					flag_redo = true;
+			}
+		}
+
+		ldegrees.push_back(factor * real_ldegree);
+		ldegree_redo.push_back(flag_redo);
+	}
+
+	int degbound = order-std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
+	// Second round: determine the remaining positive ldegrees by the series
+	// method.
+	// here we can ignore ldegrees larger than degbound
+	size_t j = 0;
+	for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
+		if ( ldegree_redo[j] ) {
+			ex expon = it->coeff;
+			int factor = 1;
+			ex buf;
+			if (expon.info(info_flags::integer)) {
+				buf = it->rest;
+				factor = ex_to<numeric>(expon).to_int();
+			} else {
+				buf = recombine_pair_to_ex(*it);
+			}
+			int real_ldegree = 0;
 			int orderloop = 0;
 			do {
 				orderloop++;
 				real_ldegree = buf.series(r, orderloop, options).ldegree(sym);
-			} while (real_ldegree == orderloop);
+			} while ((real_ldegree == orderloop)
+					&& ( factor*real_ldegree < degbound));
+			ldegrees[j] = factor * real_ldegree;
+			degbound -= factor * real_ldegree;
 		}
-
-		ldegrees.push_back(factor * real_ldegree);
+		j++;
 	}
 
 	int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
@@ -1042,6 +1116,29 @@ ex power::series(const relational & r, int order, unsigned options) const
 		must_expand_basis = true;
 	}
 
+	bool exponent_is_regular = true;
+	try {
+		exponent.subs(r, subs_options::no_pattern);
+	} catch (pole_error) {
+		exponent_is_regular = false;
+	}
+
+	if (!exponent_is_regular) {
+		ex l = exponent*log(basis);
+		// this == exp(l);
+		ex le = l.series(r, order, options);
+		// Note: expanding exp(l) won't help, since that will attempt
+		// Taylor expansion, and fail (because exponent is "singular")
+		// Still l itself might be expanded in Taylor series.
+		// Examples:
+		// sin(x)/x*log(cos(x))
+		// 1/x*log(1 + x)
+		return exp(le).series(r, order, options);
+		// Note: if l happens to have a Laurent expansion (with
+		// negative powers of (var - point)), expanding exp(le)
+		// will barf (which is The Right Thing).
+	}
+
 	// Is the expression of type something^(-int)?
 	if (!must_expand_basis && !exponent.info(info_flags::negint)
 	 && (!is_a<add>(basis) || !is_a<numeric>(exponent)))
@@ -1072,7 +1169,10 @@ ex power::series(const relational & r, int order, unsigned options) const
 	}
 	const ex& sym = r.lhs();
 	// find existing minimal degree
-	int real_ldegree = basis.expand().ldegree(sym-r.rhs());
+	ex eb = basis.expand();
+	int real_ldegree = 0;
+	if (eb.info(info_flags::rational_function))
+		real_ldegree = eb.ldegree(sym-r.rhs());
 	if (real_ldegree == 0) {
 		int orderloop = 0;
 		do {
@@ -1203,12 +1303,10 @@ ex ex::series(const ex & r, int order, unsigned options) const
 	else
 		throw (std::logic_error("ex::series(): expansion point has unknown type"));
 	
-	try {
-		e = bp->series(rel_, order, options);
-	} catch (std::exception &x) {
-		throw (std::logic_error(std::string("unable to compute series (") + x.what() + ")"));
-	}
+	e = bp->series(rel_, order, options);
 	return e;
 }
 
+GINAC_BIND_UNARCHIVER(pseries);
+
 } // namespace GiNaC