X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=01a57569125cf1675a7dfbd522e21740c6cd3bc0;hp=6f72964e67c886412ebf6e29bd58b817dba331cc;hb=6335b2cfb6f3f9ec339ff683bed415a5879d7909;hpb=695f6ae955ec530cded8f21efd5569df39447f76;ds=sidebyside

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index 6f72964e..01a57569 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-2008 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
@@ -61,7 +61,7 @@
  *
  *  You should have received a copy of the GNU General Public License
  *  along with this program; if not, write to the Free Software
- *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
 #include <sstream>
@@ -320,7 +320,7 @@ cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
 
 
 // forward declaration needed by function Li_projection and C below
-numeric S_num(int n, int p, const numeric& x);
+const cln::cl_N S_num(int n, int p, const cln::cl_N& x);
 
 
 // helper function for classical polylog Li
@@ -371,7 +371,7 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
 		} else {
 			cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
 			for (int j=0; j<n-1; j++) {
-				result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
+				result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
 				                  * cln::expt(cln::log(x), j) / cln::factorial(j);
 			}
 			return result;
@@ -379,15 +379,14 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
 	}
 }
 
-
 // helper function for classical polylog Li
-numeric Lin_numeric(int n, const numeric& x)
+const cln::cl_N Lin_numeric(const int n, const cln::cl_N& x)
 {
 	if (n == 1) {
 		// just a log
-		return -cln::log(1-x.to_cl_N());
+		return -cln::log(1-x);
 	}
-	if (x.is_zero()) {
+	if (zerop(x)) {
 		return 0;
 	}
 	if (x == 1) {
@@ -398,12 +397,11 @@ numeric Lin_numeric(int n, const numeric& x)
 		// [Kol] (2.22)
 		return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
 	}
-	if (abs(x.real()) < 0.4 && abs(abs(x)-1) < 0.01) {
-		cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
-		cln::cl_N result = -cln::expt(cln::log(x_), n-1) * cln::log(1-x_) / cln::factorial(n-1);
+	if (cln::abs(realpart(x)) < 0.4 && cln::abs(cln::abs(x)-1) < 0.01) {
+		cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
 		for (int j=0; j<n-1; j++) {
-			result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x_).to_cl_N())
-				* cln::expt(cln::log(x_), j) / cln::factorial(j);
+			result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
+				* cln::expt(cln::log(x), j) / cln::factorial(j);
 		}
 		return result;
 	}
@@ -411,11 +409,11 @@ numeric Lin_numeric(int n, const numeric& x)
 	// what is the desired float format?
 	// first guess: default format
 	cln::float_format_t prec = cln::default_float_format;
-	const cln::cl_N value = x.to_cl_N();
+	const cln::cl_N value = x;
 	// second guess: the argument's format
-	if (!x.real().is_rational())
+	if (!instanceof(realpart(x), cln::cl_RA_ring))
 		prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
-	else if (!x.imag().is_rational())
+	else if (!instanceof(imagpart(x), cln::cl_RA_ring))
 		prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 	
 	// [Kol] (5.15)
@@ -441,7 +439,7 @@ numeric Lin_numeric(int n, const numeric& x)
 		cln::cl_N add;
 		for (int j=0; j<n-1; j++) {
 			add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
-			            * Lin_numeric(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
+			            * Lin_numeric(n-j,1) * cln::expt(cln::log(-value),j) / cln::factorial(j);
 		}
 		result = result - add;
 		return result;
@@ -471,7 +469,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 +484,18 @@ 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) || cln::zerop(t[0]) || flag_accidental_zero );
 
 	return t[0];
 }
@@ -980,9 +983,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;
@@ -992,10 +994,62 @@ ex G_numeric(const lst& x, const lst& s, const ex& y)
 	if (x.op(x.nops()-1).is_zero()) {
 		need_trafo = true;
 	}
-	if (depth == 1 && !need_trafo) {
+	if (depth == 1 && x.nops() == 2 && !need_trafo) {
 		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 +1125,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;
@@ -1376,12 +1378,18 @@ static ex Li_evalf(const ex& m_, const ex& x_)
 	// classical polylogs
 	if (m_.info(info_flags::posint)) {
 		if (x_.info(info_flags::numeric)) {
-			return Lin_numeric(ex_to<numeric>(m_).to_int(), ex_to<numeric>(x_));
+			int m__ = ex_to<numeric>(m_).to_int();
+			const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
+			const cln::cl_N result = Lin_numeric(m__, x__);
+			return numeric(result);
 		} else {
 			// try to numerically evaluate second argument
 			ex x_val = x_.evalf();
 			if (x_val.info(info_flags::numeric)) {
-				return Lin_numeric(ex_to<numeric>(m_).to_int(), ex_to<numeric>(x_val));
+				int m__ = ex_to<numeric>(m_).to_int();
+				const cln::cl_N x__ = ex_to<numeric>(x_val).to_cl_N();
+				const cln::cl_N result = Lin_numeric(m__, x__);
+				return numeric(result);
 			}
 		}
 	}
@@ -1495,7 +1503,10 @@ static ex Li_eval(const ex& m_, const ex& x_)
 		}
 	}
 	if (m_.info(info_flags::posint) && x_.info(info_flags::numeric) && !x_.info(info_flags::crational)) {
-		return Lin_numeric(ex_to<numeric>(m_).to_int(), ex_to<numeric>(x_));
+		int m__ = ex_to<numeric>(m_).to_int();
+		const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
+		const cln::cl_N result = Lin_numeric(m__, x__);
+		return numeric(result);
 	}
 
 	return Li(m_, x_).hold();
@@ -1504,9 +1515,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()
 }
 
 
@@ -1688,10 +1727,10 @@ cln::cl_N C(int n, int p)
 			if (k == 0) {
 				if (n & 1) {
 					if (j & 1) {
-						result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
+						result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
 					}
 					else {
-						result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
+						result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
 					}
 				}
 			}
@@ -1699,23 +1738,23 @@ cln::cl_N C(int n, int p)
 				if (k & 1) {
 					if (j & 1) {
 						result = result + cln::factorial(n+k-1)
-						                  * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+						                  * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
 						                  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 					}
 					else {
 						result = result - cln::factorial(n+k-1)
-						                  * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+						                  * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
 						                  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 					}
 				}
 				else {
 					if (j & 1) {
-						result = result - cln::factorial(n+k-1) * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+						result = result - cln::factorial(n+k-1) * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
 						                  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 					}
 					else {
 						result = result + cln::factorial(n+k-1)
-						                  * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+						                  * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
 						                  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 					}
 				}
@@ -1777,10 +1816,20 @@ cln::cl_N b_k(int k)
 // helper function for S(n,p,x)
 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 {
+	static cln::float_format_t oldprec = cln::default_float_format;
+
 	if (p==1) {
 		return Li_projection(n+1, x, prec);
 	}
-	
+
+	// precision has changed, we need to clear lookup table Yn
+	if ( oldprec != prec ) {
+		Yn.clear();
+		ynsize = 0;
+		ynlength = 100;
+		oldprec = prec;
+	}
+		
 	// check if precalculated values are sufficient
 	if (p > ynsize+1) {
 		for (int i=ynsize; i<p-1; i++) {
@@ -1828,7 +1877,7 @@ cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format
 				res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
 				              * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
 			}
-			result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
+			result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
 		}
 
 		return result;
@@ -1839,7 +1888,7 @@ cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format
 
 
 // helper function for S(n,p,x)
-numeric S_num(int n, int p, const numeric& x)
+const cln::cl_N S_num(int n, int p, const cln::cl_N& x)
 {
 	if (x == 1) {
 		if (n == 1) {
@@ -1875,15 +1924,15 @@ numeric S_num(int n, int p, const numeric& x)
 	// what is the desired float format?
 	// first guess: default format
 	cln::float_format_t prec = cln::default_float_format;
-	const cln::cl_N value = x.to_cl_N();
+	const cln::cl_N value = x;
 	// second guess: the argument's format
-	if (!x.real().is_rational())
+	if (!instanceof(realpart(value), cln::cl_RA_ring))
 		prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
-	else if (!x.imag().is_rational())
+	else if (!instanceof(imagpart(value), cln::cl_RA_ring))
 		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);
@@ -1892,9 +1941,9 @@ numeric S_num(int n, int p, const numeric& x)
 			cln::cl_N res2;
 			for (int r=0; r<p; r++) {
 				res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
-				              * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
+				              * S_num(p-r,n-s,1-value) / cln::factorial(r);
 			}
-			result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
+			result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
 		}
 
 		return result;
@@ -1909,7 +1958,7 @@ numeric S_num(int n, int p, const numeric& x)
 			for (int r=0; r<=s; r++) {
 				result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
 				                  / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
-				                  * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
+				                  * S_num(n+s-r,p-s,cln::recip(value));
 			}
 		}
 		result = result * cln::expt(cln::cl_I(-1),n);
@@ -1945,12 +1994,18 @@ numeric S_num(int n, int p, const numeric& x)
 static ex S_evalf(const ex& n, const ex& p, const ex& x)
 {
 	if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
+		const int n_ = ex_to<numeric>(n).to_int();
+		const int p_ = ex_to<numeric>(p).to_int();
 		if (is_a<numeric>(x)) {
-			return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
+			const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
+			const cln::cl_N result = S_num(n_, p_, x_);
+			return numeric(result);
 		} else {
 			ex x_val = x.evalf();
 			if (is_a<numeric>(x_val)) {
-				return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x_val));
+				const cln::cl_N x_val_ = ex_to<numeric>(x_val).to_cl_N();
+				const cln::cl_N result = S_num(n_, p_, x_val_);
+				return numeric(result);
 			}
 		}
 	}
@@ -1975,7 +2030,11 @@ static ex S_eval(const ex& n, const ex& p, const ex& x)
 			return Li(n+1, x);
 		}
 		if (x.info(info_flags::numeric) && (!x.info(info_flags::crational))) {
-			return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
+			int n_ = ex_to<numeric>(n).to_int();
+			int p_ = ex_to<numeric>(p).to_int();
+			const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
+			const cln::cl_N result = S_num(n_, p_, x_);
+			return numeric(result);
 		}
 	}
 	if (n.is_zero()) {
@@ -1988,9 +2047,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 +2512,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 +2638,107 @@ 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(1);
+						}
+						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(0);
+						}
+						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 = H(lst(1), arg).hold() * H(newparameter, arg).hold();
+					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();
+					}
+					res = recursion(res).expand() / firstzero;
+					return unify(res);
+				}
+			}
+		}
+		return e;
+	}
+};
+
+
 // do x -> 1/x transformation
 struct map_trafo_H_1overx : public map_function
 {
@@ -2822,6 +3052,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 +3060,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 +3101,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 +3115,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 +3124,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 +3770,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;
 				}