X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=cd3511ad0e43c9c04c3ef2c15af42b0364d6e5c5;hp=f5d33ec79cb15c54820d22f6bbac9fe9e00092d5;hb=f69a24a4fd6caf42ef773d1cef21562a8afa068a;hpb=cbd91ff9206af199888748f0738664d320f456dd

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index f5d33ec7..cd3511ad 100644
--- a/ginac/inifcns_nstdsums.cpp
+++ b/ginac/inifcns_nstdsums.cpp
@@ -47,7 +47,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2007 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
@@ -471,7 +471,13 @@ namespace {
 // performs the actual series summation for multiple polylogarithms
 cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
 {
+	// ensure all x <> 0.
+	for (std::vector<cln::cl_N>::const_iterator it = x.begin(); it != x.end(); ++it) {
+		if ( *it == 0 ) return cln::cl_float(0, cln::float_format(Digits));
+	}
+
 	const int j = s.size();
+	bool flag_accidental_zero = false;
 
 	std::vector<cln::cl_N> t(j);
 	cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
@@ -480,19 +486,13 @@ cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl
 	int q = 0;
 	do {
 		t0buf = t[0];
-		// do it once ...
 		q++;
 		t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
 		for (int k=j-2; k>=0; k--) {
+			flag_accidental_zero = cln::zerop(t[k+1]);
 			t[k] = t[k] + t[k+1] * cln::expt(x[k], q+j-1-k) / cln::expt(cln::cl_I(q+j-1-k), s[k]);
 		}
-		// ... and do it again (to avoid premature drop out due to special arguments)
-		q++;
-		t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
-		for (int k=j-2; k>=0; k--) {
-			t[k] = t[k] + t[k+1] * cln::expt(x[k], q+j-1-k) / cln::expt(cln::cl_I(q+j-1-k), s[k]);
-		}
-	} while (t[0] != t0buf);
+	} while ( (t[0] != t0buf) || flag_accidental_zero );
 
 	return t[0];
 }
@@ -980,9 +980,8 @@ ex G_numeric(const lst& x, const lst& s, const ex& y)
 	for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
 		if (!(*it).is_zero()) {
 			++depth;
-			if (abs(*it) - y < -pow(10,-Digits+2)) {
+			if (abs(*it) - y < -pow(10,-Digits+1)) {
 				need_trafo = true;
-				break;
 			}
 			if (abs((abs(*it) - y)/y) < 0.01) {
 				need_hoelder = true;
@@ -996,6 +995,58 @@ ex G_numeric(const lst& x, const lst& s, const ex& y)
 		return -Li(x.nops(), y / x.op(x.nops()-1)).evalf();
 	}
 	
+	// do acceleration transformation (hoelder convolution [BBB])
+	if (need_hoelder) {
+		
+		ex result;
+		const int size = x.nops();
+		lst newx;
+		for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
+			newx.append(*it / y);
+		}
+		
+		for (int r=0; r<=size; ++r) {
+			ex buffer = pow(-1, r);
+			ex p = 2;
+			bool adjustp;
+			do {
+				adjustp = false;
+				for (lst::const_iterator it = newx.begin(); it != newx.end(); ++it) {
+					if (*it == 1/p) {
+						p += (3-p)/2; 
+						adjustp = true;
+						continue;
+					}
+				}
+			} while (adjustp);
+			ex q = p / (p-1);
+			lst qlstx;
+			lst qlsts;
+			for (int j=r; j>=1; --j) {
+				qlstx.append(1-newx.op(j-1));
+				if (newx.op(j-1).info(info_flags::real) && newx.op(j-1) > 1 && newx.op(j-1) <= 2) {
+					qlsts.append( s.op(j-1));
+				} else {
+					qlsts.append( -s.op(j-1));
+				}
+			}
+			if (qlstx.nops() > 0) {
+				buffer *= G_numeric(qlstx, qlsts, 1/q);
+			}
+			lst plstx;
+			lst plsts;
+			for (int j=r+1; j<=size; ++j) {
+				plstx.append(newx.op(j-1));
+				plsts.append(s.op(j-1));
+			}
+			if (plstx.nops() > 0) {
+				buffer *= G_numeric(plstx, plsts, 1/p);
+			}
+			result += buffer;
+		}
+		return result;
+	}
+	
 	// convergence transformation
 	if (need_trafo) {
 
@@ -1071,58 +1122,6 @@ ex G_numeric(const lst& x, const lst& s, const ex& y)
 		return result;
 	}
 
-	// do acceleration transformation (hoelder convolution [BBB])
-	if (need_hoelder) {
-		
-		ex result;
-		const int size = x.nops();
-		lst newx;
-		for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
-			newx.append(*it / y);
-		}
-		
-		for (int r=0; r<=size; ++r) {
-			ex buffer = pow(-1, r);
-			ex p = 2;
-			bool adjustp;
-			do {
-				adjustp = false;
-				for (lst::const_iterator it = newx.begin(); it != newx.end(); ++it) {
-					if (*it == 1/p) {
-						p += (3-p)/2; 
-						adjustp = true;
-						continue;
-					}
-				}
-			} while (adjustp);
-			ex q = p / (p-1);
-			lst qlstx;
-			lst qlsts;
-			for (int j=r; j>=1; --j) {
-				qlstx.append(1-newx.op(j-1));
-				if (newx.op(j-1).info(info_flags::real) && newx.op(j-1) > 1 && newx.op(j-1) <= 2) {
-					qlsts.append( s.op(j-1));
-				} else {
-					qlsts.append( -s.op(j-1));
-				}
-			}
-			if (qlstx.nops() > 0) {
-				buffer *= G_numeric(qlstx, qlsts, 1/q);
-			}
-			lst plstx;
-			lst plsts;
-			for (int j=r+1; j<=size; ++j) {
-				plstx.append(newx.op(j-1));
-				plsts.append(s.op(j-1));
-			}
-			if (plstx.nops() > 0) {
-				buffer *= G_numeric(plstx, plsts, 1/p);
-			}
-			result += buffer;
-		}
-		return result;
-	}
-	
 	// do summation
 	lst newx;
 	lst m;
@@ -1504,9 +1503,37 @@ static ex Li_eval(const ex& m_, const ex& x_)
 
 static ex Li_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
 {
-	epvector seq;
-	seq.push_back(expair(Li(m, x), 0));
-	return pseries(rel, seq);
+	if (is_a<lst>(m) || is_a<lst>(x)) {
+		// multiple polylog
+		epvector seq;
+		seq.push_back(expair(Li(m, x), 0));
+		return pseries(rel, seq);
+	}
+	
+	// classical polylog
+	const ex x_pt = x.subs(rel, subs_options::no_pattern);
+	if (m.info(info_flags::numeric) && x_pt.info(info_flags::numeric)) {
+		// First special case: x==0 (derivatives have poles)
+		if (x_pt.is_zero()) {
+			const symbol s;
+			ex ser;
+			// manually construct the primitive expansion
+			for (int i=1; i<order; ++i)
+				ser += pow(s,i) / pow(numeric(i), m);
+			// substitute the argument's series expansion
+			ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
+			// maybe that was terminating, so add a proper order term
+			epvector nseq;
+			nseq.push_back(expair(Order(_ex1), order));
+			ser += pseries(rel, nseq);
+			// reexpanding it will collapse the series again
+			return ser.series(rel, order);
+		}
+		// TODO special cases: x==1 (branch point) and x real, >=1 (branch cut)
+		throw std::runtime_error("Li_series: don't know how to do the series expansion at this point!");
+	}
+	// all other cases should be safe, by now:
+	throw do_taylor();  // caught by function::series()
 }
 
 
@@ -1883,7 +1910,7 @@ numeric S_num(int n, int p, const numeric& x)
 		prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 
 	// [Kol] (5.3)
-	if ((cln::realpart(value) < -0.5) || (n == 0)) {
+	if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95))) {
 
 		cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
 		                   * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
@@ -1988,9 +2015,48 @@ static ex S_eval(const ex& n, const ex& p, const ex& x)
 
 static ex S_series(const ex& n, const ex& p, const ex& x, const relational& rel, int order, unsigned options)
 {
-	epvector seq;
-	seq.push_back(expair(S(n, p, x), 0));
-	return pseries(rel, seq);
+	if (p == _ex1) {
+		return Li(n+1, x).series(rel, order, options);
+	}
+
+	const ex x_pt = x.subs(rel, subs_options::no_pattern);
+	if (n.info(info_flags::posint) && p.info(info_flags::posint) && x_pt.info(info_flags::numeric)) {
+		// First special case: x==0 (derivatives have poles)
+		if (x_pt.is_zero()) {
+			const symbol s;
+			ex ser;
+			// manually construct the primitive expansion
+			// subsum = Euler-Zagier-Sum is needed
+			// dirty hack (slow ...) calculation of subsum:
+			std::vector<ex> presubsum, subsum;
+			subsum.push_back(0);
+			for (int i=1; i<order-1; ++i) {
+				subsum.push_back(subsum[i-1] + numeric(1, i));
+			}
+			for (int depth=2; depth<p; ++depth) {
+				presubsum = subsum;
+				for (int i=1; i<order-1; ++i) {
+					subsum[i] = subsum[i-1] + numeric(1, i) * presubsum[i-1];
+				}
+			}
+				
+			for (int i=1; i<order; ++i) {
+				ser += pow(s,i) / pow(numeric(i), n+1) * subsum[i-1];
+			}
+			// substitute the argument's series expansion
+			ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
+			// maybe that was terminating, so add a proper order term
+			epvector nseq;
+			nseq.push_back(expair(Order(_ex1), order));
+			ser += pseries(rel, nseq);
+			// reexpanding it will collapse the series again
+			return ser.series(rel, order);
+		}
+		// TODO special cases: x==1 (branch point) and x real, >=1 (branch cut)
+		throw std::runtime_error("S_series: don't know how to do the series expansion at this point!");
+	}
+	// all other cases should be safe, by now:
+	throw do_taylor();  // caught by function::series()
 }
 
 
@@ -2568,7 +2634,7 @@ struct map_trafo_H_1mx : public map_function
 					if (allthesame) {
 						lst newparameter;
 						for (int i=parameter.nops(); i>0; i--) {
-							newparameter.append(0);
+							newparameter.append(1);
 						}
 						return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
 					}
@@ -2584,7 +2650,7 @@ struct map_trafo_H_1mx : public map_function
 					if (allthesame) {
 						lst newparameter;
 						for (int i=parameter.nops(); i>0; i--) {
-							newparameter.append(1);
+							newparameter.append(0);
 						}
 						return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
 					}
@@ -2614,7 +2680,7 @@ struct map_trafo_H_1mx : public map_function
 					// leading one
 					map_trafo_H_1mx recursion;
 					map_trafo_H_mult unify;
-					ex res;
+					ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
 					int firstzero = 0;
 					while (parameter.op(firstzero) == 1) {
 						firstzero++;
@@ -2631,11 +2697,9 @@ struct map_trafo_H_1mx : public map_function
 						}
 						res -= H(newparameter, arg).hold();
 					}
-					return (unify((-H(lst(0), 1-arg).hold() * recursion(H(newparameter, arg).hold())).expand()) +
-							recursion(res)) / firstzero;
-
+					res = recursion(res).expand() / firstzero;
+					return unify(res);
 				}
-
 			}
 		}
 		return e;
@@ -3033,7 +3097,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
 		
 		// check transformations for 0.95 <= |x| < 2.0
 		
-		// |(1-x)/(1+x)| < 0.9 -> circular area with center=9,53+0i and radius=9.47
+		// |(1-x)/(1+x)| < 0.9 -> circular area with center=9.53+0i and radius=9.47
 		if (cln::abs(x-9.53) <= 9.47) {
 			// x -> (1-x)/(1+x)
 			map_trafo_H_1mxt1px trafo;