X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpseries.cpp;h=d799e62b132a576b23430298a2237b09b0b5eef7;hp=48dfa1ebbf7d98337954ccac7ef116f9f7e1fef2;hb=846e62e85c483f74d1e253c7841050ea11d70451;hpb=b301f03a61bc9f72b27940ca7fe1f8d0b343a4e2

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index 48dfa1eb..d799e62b 100644
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -4,7 +4,7 @@
  *  methods for series expansion. */
 
 /*
- *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2004 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,7 +21,7 @@
  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  */
 
-#include <iostream>
+#include <numeric>
 #include <stdexcept>
 
 #include "pseries.h"
@@ -33,13 +33,17 @@
 #include "relational.h"
 #include "operators.h"
 #include "symbol.h"
-#include "print.h"
 #include "archive.h"
 #include "utils.h"
 
 namespace GiNaC {
 
-GINAC_IMPLEMENT_REGISTERED_CLASS(pseries, basic)
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(pseries, basic,
+  print_func<print_context>(&pseries::do_print).
+  print_func<print_latex>(&pseries::do_print_latex).
+  print_func<print_tree>(&pseries::do_print_tree).
+  print_func<print_python>(&pseries::do_print_python).
+  print_func<print_python_repr>(&pseries::do_print_python_repr))
 
 
 /*
@@ -108,106 +112,113 @@ DEFAULT_UNARCHIVE(pseries)
 // functions overriding virtual functions from base classes
 //////////
 
-void pseries::print(const print_context & c, unsigned level) const
+void pseries::print_series(const print_context & c, const char *openbrace, const char *closebrace, const char *mul_sym, const char *pow_sym, unsigned level) const
 {
-	if (is_a<print_tree>(c)) {
+	if (precedence() <= level)
+		c.s << '(';
+		
+	// objects of type pseries must not have any zero entries, so the
+	// trivial (zero) pseries needs a special treatment here:
+	if (seq.empty())
+		c.s << '0';
 
-		c.s << std::string(level, ' ') << class_name()
-		    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
-		    << std::endl;
-		unsigned delta_indent = static_cast<const print_tree &>(c).delta_indent;
-		size_t num = seq.size();
-		for (size_t i=0; i<num; ++i) {
-			seq[i].rest.print(c, level + delta_indent);
-			seq[i].coeff.print(c, level + delta_indent);
-			c.s << std::string(level + delta_indent, ' ') << "-----" << std::endl;
-		}
-		var.print(c, level + delta_indent);
-		point.print(c, level + delta_indent);
+	epvector::const_iterator i = seq.begin(), end = seq.end();
+	while (i != end) {
 
-	} else if (is_a<print_python_repr>(c)) {
-		c.s << class_name() << "(relational(";
-		var.print(c);
-		c.s << ',';
-		point.print(c);
-		c.s << "),[";
-		size_t num = seq.size();
-		for (size_t i=0; i<num; ++i) {
-			if (i)
-				c.s << ',';
-			c.s << '(';
-			seq[i].rest.print(c);
-			c.s << ',';
-			seq[i].coeff.print(c);
-			c.s << ')';
-		}
-		c.s << "])";
-	} else {
+		// print a sign, if needed
+		if (i != seq.begin())
+			c.s << '+';
 
-		if (precedence() <= level)
-			c.s << "(";
-		
-		std::string par_open = is_a<print_latex>(c) ? "{(" : "(";
-		std::string par_close = is_a<print_latex>(c) ? ")}" : ")";
-		
-		// objects of type pseries must not have any zero entries, so the
-		// trivial (zero) pseries needs a special treatment here:
-		if (seq.empty())
-			c.s << '0';
-		epvector::const_iterator i = seq.begin(), end = seq.end();
-		while (i != end) {
-			// print a sign, if needed
-			if (i != seq.begin())
-				c.s << '+';
-			if (!is_order_function(i->rest)) {
-				// print 'rest', i.e. the expansion coefficient
-				if (i->rest.info(info_flags::numeric) &&
-					i->rest.info(info_flags::positive)) {
-					i->rest.print(c);
-				} else {
-					c.s << par_open;
-					i->rest.print(c);
-					c.s << par_close;
-				}
-				// print 'coeff', something like (x-1)^42
-				if (!i->coeff.is_zero()) {
-					if (is_a<print_latex>(c))
-						c.s << ' ';
-					else
-						c.s << '*';
-					if (!point.is_zero()) {
-						c.s << par_open;
-						(var-point).print(c);
-						c.s << par_close;
+		if (!is_order_function(i->rest)) {
+
+			// print 'rest', i.e. the expansion coefficient
+			if (i->rest.info(info_flags::numeric) &&
+				i->rest.info(info_flags::positive)) {
+				i->rest.print(c);
+			} else {
+				c.s << openbrace << '(';
+				i->rest.print(c);
+				c.s << ')' << closebrace;
+			}
+
+			// print 'coeff', something like (x-1)^42
+			if (!i->coeff.is_zero()) {
+				c.s << mul_sym;
+				if (!point.is_zero()) {
+					c.s << openbrace << '(';
+					(var-point).print(c);
+					c.s << ')' << closebrace;
+				} else
+					var.print(c);
+				if (i->coeff.compare(_ex1)) {
+					c.s << pow_sym;
+					c.s << openbrace;
+					if (i->coeff.info(info_flags::negative)) {
+						c.s << '(';
+						i->coeff.print(c);
+						c.s << ')';
 					} else
-						var.print(c);
-					if (i->coeff.compare(_ex1)) {
-						if (is_a<print_python>(c))
-							c.s << "**";
-						else
-							c.s << '^';
-						if (i->coeff.info(info_flags::negative)) {
-							c.s << par_open;
-							i->coeff.print(c);
-							c.s << par_close;
-						} else {
-							if (is_a<print_latex>(c)) {
-								c.s << '{';
-								i->coeff.print(c);
-								c.s << '}';
-							} else
-								i->coeff.print(c);
-						}
-					}
+						i->coeff.print(c);
+					c.s << closebrace;
 				}
-			} else
-				Order(power(var-point,i->coeff)).print(c);
-			++i;
-		}
+			}
+		} else
+			Order(power(var-point,i->coeff)).print(c);
+		++i;
+	}
 
-		if (precedence() <= level)
-			c.s << ")";
+	if (precedence() <= level)
+		c.s << ')';
+}
+
+void pseries::do_print(const print_context & c, unsigned level) const
+{
+	print_series(c, "", "", "*", "^", level);
+}
+
+void pseries::do_print_latex(const print_latex & c, unsigned level) const
+{
+	print_series(c, "{", "}", " ", "^", level);
+}
+
+void pseries::do_print_python(const print_python & c, unsigned level) const
+{
+	print_series(c, "", "", "*", "**", level);
+}
+
+void pseries::do_print_tree(const print_tree & c, unsigned level) const
+{
+	c.s << std::string(level, ' ') << class_name() << " @" << this
+	    << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
+	    << std::endl;
+	size_t num = seq.size();
+	for (size_t i=0; i<num; ++i) {
+		seq[i].rest.print(c, level + c.delta_indent);
+		seq[i].coeff.print(c, level + c.delta_indent);
+		c.s << std::string(level + c.delta_indent, ' ') << "-----" << std::endl;
+	}
+	var.print(c, level + c.delta_indent);
+	point.print(c, level + c.delta_indent);
+}
+
+void pseries::do_print_python_repr(const print_python_repr & c, unsigned level) const
+{
+	c.s << class_name() << "(relational(";
+	var.print(c);
+	c.s << ',';
+	point.print(c);
+	c.s << "),[";
+	size_t num = seq.size();
+	for (size_t i=0; i<num; ++i) {
+		if (i)
+			c.s << ',';
+		c.s << '(';
+		seq[i].rest.print(c);
+		c.s << ',';
+		seq[i].coeff.print(c);
+		c.s << ')';
 	}
+	c.s << "])";
 }
 
 int pseries::compare_same_type(const basic & other) const
@@ -396,6 +407,23 @@ ex pseries::evalf(int level) const
 	return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
 }
 
+ex pseries::conjugate() const
+{
+	epvector * newseq = conjugateepvector(seq);
+	ex newvar = var.conjugate();
+	ex newpoint = point.conjugate();
+
+	if (!newseq	&& are_ex_trivially_equal(newvar, var) && are_ex_trivially_equal(point, newpoint)) {
+		return *this;
+	}
+
+	ex result = (new pseries(newvar==newpoint, newseq ? *newseq : seq))->setflag(status_flags::dynallocated);
+	if (newseq) {
+		delete newseq;
+	}
+	return result;
+}
+
 ex pseries::subs(const exmap & m, unsigned options) const
 {
 	// If expansion variable is being substituted, convert the series to a
@@ -501,14 +529,23 @@ bool pseries::is_terminating() const
 ex basic::series(const relational & r, int order, unsigned options) const
 {
 	epvector seq;
+	const symbol &s = ex_to<symbol>(r.lhs());
+
+	// default for order-values that make no sense for Taylor expansion
+	if ((order <= 0) && this->has(s)) {
+		seq.push_back(expair(Order(_ex1), order));
+		return pseries(r, seq);
+	}
+
+	// do Taylor expansion
 	numeric fac = 1;
 	ex deriv = *this;
 	ex coeff = deriv.subs(r, subs_options::no_pattern);
-	const symbol &s = ex_to<symbol>(r.lhs());
-	
-	if (!coeff.is_zero())
+
+	if (!coeff.is_zero()) {
 		seq.push_back(expair(coeff, _ex0));
-	
+	}
+
 	int n;
 	for (n=1; n<order; ++n) {
 		fac = fac.mul(n);
@@ -737,11 +774,39 @@ ex mul::series(const relational & r, int order, unsigned options) const
 {
 	pseries acc; // Series accumulator
 
-	// Multiply with remaining terms
+	GINAC_ASSERT(is_a<symbol>(r.lhs()));
+	const ex& sym = r.lhs();
+		
+	// holds ldegrees of the series of individual factors
+	std::vector<int> ldegrees;
+
+	// find minimal degrees
 	const epvector::const_iterator itbeg = seq.begin();
 	const epvector::const_iterator itend = seq.end();
 	for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
-		ex op = recombine_pair_to_ex(*it).series(r, order, options);
+		
+		ex buf = recombine_pair_to_ex(*it);
+		
+		int real_ldegree = buf.expand().ldegree(sym-r.rhs());
+		if (real_ldegree == 0) {
+			int orderloop = 0;
+			do {
+				orderloop++;
+				real_ldegree = buf.series(r, orderloop, options).ldegree(sym);
+			} while (real_ldegree == orderloop);
+		}
+
+		ldegrees.push_back(real_ldegree);
+	}
+
+	int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
+		
+	// Multiply with remaining terms
+	std::vector<int>::const_iterator itd = ldegrees.begin();
+	for (epvector::const_iterator it=itbeg; it!=itend; ++it, ++itd) {
+
+		// do series expansion with adjusted order
+		ex op = recombine_pair_to_ex(*it).series(r, order-degsum+(*itd), options);
 
 		// Series multiplication
 		if (it==itbeg)
@@ -749,6 +814,7 @@ ex mul::series(const relational & r, int order, unsigned options) const
 		else
 			acc = ex_to<pseries>(acc.mul_series(ex_to<pseries>(op)));
 	}
+
 	return acc.mul_const(ex_to<numeric>(overall_coeff));
 }
 
@@ -795,6 +861,9 @@ ex pseries::power_const(const numeric &p, int deg) const
 	if (!(p*ldeg).is_integer())
 		throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
 
+	// adjust number of coefficients
+	deg = deg - p.to_int()*ldeg;
+	
 	// O(x^n)^(-m) is undefined
 	if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
 		throw pole_error("pseries::power_const(): division by zero",1);
@@ -803,7 +872,6 @@ ex pseries::power_const(const numeric &p, int deg) const
 	exvector co;
 	co.reserve(deg);
 	co.push_back(power(coeff(var, ldeg), p));
-	bool all_sums_zero = true;
 	for (int i=1; i<deg; ++i) {
 		ex sum = _ex0;
 		for (int j=1; j<=i; ++j) {
@@ -814,8 +882,6 @@ ex pseries::power_const(const numeric &p, int deg) const
 			} else
 				sum += (p * j - (i - j)) * co[i - j] * c;
 		}
-		if (!sum.is_zero())
-			all_sums_zero = false;
 		co.push_back(sum / coeff(var, ldeg) / i);
 	}
 	
@@ -830,8 +896,9 @@ ex pseries::power_const(const numeric &p, int deg) const
 			break;
 		}
 	}
-	if (!higher_order && !all_sums_zero)
+	if (!higher_order)
 		new_seq.push_back(expair(Order(_ex1), p * ldeg + deg));
+
 	return pseries(relational(var,point), new_seq);
 }
 
@@ -865,11 +932,11 @@ ex power::series(const relational & r, int order, unsigned options) const
 	} catch (pole_error) {
 		must_expand_basis = true;
 	}
-		
+
 	// Is the expression of type something^(-int)?
 	if (!must_expand_basis && !exponent.info(info_flags::negint))
 		return basic::series(r, order, options);
-		
+
 	// Is the expression of type 0^something?
 	if (!must_expand_basis && !basis.subs(r, subs_options::no_pattern).is_zero())
 		return basic::series(r, order, options);
@@ -885,8 +952,32 @@ ex power::series(const relational & r, int order, unsigned options) const
 	}
 
 	// No, expand basis into series
-	ex e = basis.series(r, order, options);
-	return ex_to<pseries>(e).power_const(ex_to<numeric>(exponent), order);
+
+	int intexp = ex_to<numeric>(exponent).to_int();
+	const ex& sym = r.lhs();
+	// find existing minimal degree
+	int real_ldegree = basis.expand().ldegree(sym-r.rhs());
+	if (real_ldegree == 0) {
+		int orderloop = 0;
+		do {
+			orderloop++;
+			real_ldegree = basis.series(r, orderloop, options).ldegree(sym);
+		} while (real_ldegree == orderloop);
+	}
+
+	ex e = basis.series(r, order + real_ldegree*(1-intexp), options);
+	
+	ex result;
+	try {
+		result = ex_to<pseries>(e).power_const(intexp, order);
+	}
+	catch (pole_error) {
+		epvector ser;
+		ser.push_back(expair(Order(_ex1), order));
+		result = pseries(r, ser);
+	}
+
+	return result;
 }