X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=5215dd8dba570c9d0379b0890b376651f432abd8;hp=3d1f8bcf45b37b721f090bfc962c4634ebe262dd;hb=8bf3668614986630b9a2825d4fde73a0c1f2c31c;hpb=da64e515abf7243bc4c84ca3631470931c4e6691

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index 3d1f8bcf..5215dd8d 100644
--- a/ginac/inifcns_nstdsums.cpp
+++ b/ginac/inifcns_nstdsums.cpp
@@ -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--) {
-			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--) {
+			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]);
 		}
-	} while (t[0] != t0buf);
+	} while ( (t[0] != t0buf) || flag_accidental_zero );
 
 	return t[0];
 }
@@ -1504,9 +1504,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()
 }
 
 
@@ -1988,9 +2016,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()
 }
 
 
@@ -2414,6 +2481,37 @@ ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
 }
 
 
+// do integration [ReV] (49)
+// put parameter 1 in front of existing parameters
+ex trafo_H_prepend_one(const ex& e, const ex& arg)
+{
+	ex h;
+	std::string name;
+	if (is_a<function>(e)) {
+		name = ex_to<function>(e).get_name();
+	}
+	if (name == "H") {
+		h = e;
+	} else {
+		for (int i=0; i<e.nops(); i++) {
+			if (is_a<function>(e.op(i))) {
+				std::string name = ex_to<function>(e.op(i)).get_name();
+				if (name == "H") {
+					h = e.op(i);
+				}
+			}
+		}
+	}
+	if (h != 0) {
+		lst newparameter = ex_to<lst>(h.op(0));
+		newparameter.prepend(1);
+		return e.subs(h == H(newparameter, h.op(1)).hold());
+	} else {
+		return e * H(lst(1),1-arg).hold();
+	}
+}
+
+
 // do integration [ReV] (55)
 // put parameter -1 in front of existing parameters
 ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
@@ -2509,6 +2607,109 @@ ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
 }
 
 
+// do x -> 1-x transformation
+struct map_trafo_H_1mx : public map_function
+{
+	ex operator()(const ex& e)
+	{
+		if (is_a<add>(e) || is_a<mul>(e)) {
+			return e.map(*this);
+		}
+		
+		if (is_a<function>(e)) {
+			std::string name = ex_to<function>(e).get_name();
+			if (name == "H") {
+
+				lst parameter = ex_to<lst>(e.op(0));
+				ex arg = e.op(1);
+
+				// special cases if all parameters are either 0, 1 or -1
+				bool allthesame = true;
+				if (parameter.op(0) == 0) {
+					for (int i=1; i<parameter.nops(); i++) {
+						if (parameter.op(i) != 0) {
+							allthesame = false;
+							break;
+						}
+					}
+					if (allthesame) {
+						lst newparameter;
+						for (int i=parameter.nops(); i>0; i--) {
+							newparameter.append(0);
+						}
+						return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
+					}
+				} else if (parameter.op(0) == -1) {
+					throw std::runtime_error("map_trafo_H_1mx: cannot handle weights equal -1!");
+				} else {
+					for (int i=1; i<parameter.nops(); i++) {
+						if (parameter.op(i) != 1) {
+							allthesame = false;
+							break;
+						}
+					}
+					if (allthesame) {
+						lst newparameter;
+						for (int i=parameter.nops(); i>0; i--) {
+							newparameter.append(1);
+						}
+						return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
+					}
+				}
+
+				lst newparameter = parameter;
+				newparameter.remove_first();
+
+				if (parameter.op(0) == 0) {
+
+					// leading zero
+					ex res = convert_H_to_zeta(parameter);
+					//ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
+					map_trafo_H_1mx recursion;
+					ex buffer = recursion(H(newparameter, arg).hold());
+					if (is_a<add>(buffer)) {
+						for (int i=0; i<buffer.nops(); i++) {
+							res -= trafo_H_prepend_one(buffer.op(i), arg);
+						}
+					} else {
+						res -= trafo_H_prepend_one(buffer, arg);
+					}
+					return res;
+
+				} else {
+
+					// leading one
+					map_trafo_H_1mx recursion;
+					map_trafo_H_mult unify;
+					ex res;
+					int firstzero = 0;
+					while (parameter.op(firstzero) == 1) {
+						firstzero++;
+					}
+					for (int i=firstzero-1; i<parameter.nops()-1; i++) {
+						lst newparameter;
+						int j=0;
+						for (; j<=i; j++) {
+							newparameter.append(parameter[j+1]);
+						}
+						newparameter.append(1);
+						for (; j<parameter.nops()-1; j++) {
+							newparameter.append(parameter[j+1]);
+						}
+						res -= H(newparameter, arg).hold();
+					}
+					return (unify((-H(lst(0), 1-arg).hold() * recursion(H(newparameter, arg).hold())).expand()) +
+							recursion(res)) / firstzero;
+
+				}
+
+			}
+		}
+		return e;
+	}
+};
+
+
 // do x -> 1/x transformation
 struct map_trafo_H_1overx : public map_function
 {
@@ -2822,6 +3023,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
 			return filter(H(x1, xtemp).hold()).subs(xtemp==x2).evalf();
 		}
 		// ... and expand parameter notation
+		bool has_minus_one = false;
 		lst m;
 		for (lst::const_iterator it = morg.begin(); it != morg.end(); it++) {
 			if (*it > 1) {
@@ -2829,19 +3031,18 @@ static ex H_evalf(const ex& x1, const ex& x2)
 					m.append(0);
 				}
 				m.append(1);
-			} else if (*it < -1) {
+			} else if (*it <= -1) {
 				for (ex count=*it+1; count < 0; count++) {
 					m.append(0);
 				}
 				m.append(-1);
+				has_minus_one = true;
 			} else {
 				m.append(*it);
 			}
 		}
 
-		// since the transformations produce a lot of terms, they are only efficient for
-		// argument near one.
-		// no transformation needed -> do summation
+		// do summation
 		if (cln::abs(x) < 0.95) {
 			lst m_lst;
 			lst s_lst;
@@ -2871,6 +3072,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
 			}
 		}
 
+		symbol xtemp("xtemp");
 		ex res = 1;	
 		
 		// ensure that the realpart of the argument is positive
@@ -2884,14 +3086,8 @@ static ex H_evalf(const ex& x1, const ex& x2)
 			}
 		}
 
-		// choose transformations
-		symbol xtemp("xtemp");
-		if (cln::abs(x-1) < 1.4142) {
-			// x -> (1-x)/(1+x)
-			map_trafo_H_1mxt1px trafo;
-			res *= trafo(H(m, xtemp));
-		} else {
-			// x -> 1/x
+		// x -> 1/x
+		if (cln::abs(x) >= 2.0) {
 			map_trafo_H_1overx trafo;
 			res *= trafo(H(m, xtemp));
 			if (cln::imagpart(x) <= 0) {
@@ -2899,17 +3095,25 @@ static ex H_evalf(const ex& x1, const ex& x2)
 			} else {
 				res = res.subs(H_polesign == I*Pi);
 			}
+			return res.subs(xtemp == numeric(x)).evalf();
+		}
+		
+		// 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
+		if (cln::abs(x-9.53) <= 9.47) {
+			// x -> (1-x)/(1+x)
+			map_trafo_H_1mxt1px trafo;
+			res *= trafo(H(m, xtemp));
+		} else {
+			// x -> 1-x
+			if (has_minus_one) {
+				map_trafo_H_convert_to_Li filter;
+				return filter(H(m, numeric(x)).hold()).evalf();
+			}
+			map_trafo_H_1mx trafo;
+			res *= trafo(H(m, xtemp));
 		}
-
-		// simplify result
-// TODO
-//		map_trafo_H_convert converter;
-//		res = converter(res).expand();
-//		lst ll;
-//		res.find(H(wild(1),wild(2)), ll);
-//		res.find(zeta(wild(1)), ll);
-//		res.find(zeta(wild(1),wild(2)), ll);
-//		res = res.collect(ll);
 
 		return res.subs(xtemp == numeric(x)).evalf();
 	}
@@ -3537,18 +3741,18 @@ static ex zeta1_eval(const ex& m)
 			if (y.is_zero()) {
 				return _ex_1_2;
 			}
-			if (y.is_equal(_num1)) {
+			if (y.is_equal(*_num1_p)) {
 				return zeta(m).hold();
 			}
 			if (y.info(info_flags::posint)) {
 				if (y.info(info_flags::odd)) {
 					return zeta(m).hold();
 				} else {
-					return abs(bernoulli(y)) * pow(Pi, y) * pow(_num2, y-_num1) / factorial(y);
+					return abs(bernoulli(y)) * pow(Pi, y) * pow(*_num2_p, y-(*_num1_p)) / factorial(y);
 				}
 			} else {
 				if (y.info(info_flags::odd)) {
-					return -bernoulli(_num1-y) / (_num1-y);
+					return -bernoulli((*_num1_p)-y) / ((*_num1_p)-y);
 				} else {
 					return _ex0;
 				}