X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=d85044d04bd1f311b2d26ad5ab8ccf0572d1ebe7;hb=dae52def8dd15914d32c94c2c21e2bd8419b65bb;hp=6f72964e67c886412ebf6e29bd58b817dba331cc;hpb=695f6ae955ec530cded8f21efd5569df39447f76;p=ginac.git

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index 6f72964e..d85044d0 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];
 }
@@ -515,16 +518,12 @@ cln::cl_N mLi_do_summation(const lst& m, const lst& x)
 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf);
 
 
-// holding dummy-symbols for the G/Li transformations
-std::vector<ex> gsyms;
-
-
 // type used by the transformation functions for G
 typedef std::vector<int> Gparameter;
 
 
 // G_eval1-function for G transformations
-ex G_eval1(int a, int scale)
+ex G_eval1(int a, int scale, const exvector& gsyms)
 {
 	if (a != 0) {
 		const ex& scs = gsyms[std::abs(scale)];
@@ -541,7 +540,7 @@ ex G_eval1(int a, int scale)
 
 
 // G_eval-function for G transformations
-ex G_eval(const Gparameter& a, int scale)
+ex G_eval(const Gparameter& a, int scale, const exvector& gsyms)
 {
 	// check for properties of G
 	ex sc = gsyms[std::abs(scale)];
@@ -575,7 +574,7 @@ ex G_eval(const Gparameter& a, int scale)
 		for (; it != a.end(); ++it) {
 			short_a.push_back(*it);
 		}
-		ex result = G_eval1(a.front(), scale) * G_eval(short_a, scale);
+		ex result = G_eval1(a.front(), scale, gsyms) * G_eval(short_a, scale, gsyms);
 		it = short_a.begin();
 		for (int i=1; i<count_ones; ++i) {
 			++it;
@@ -593,14 +592,14 @@ ex G_eval(const Gparameter& a, int scale)
 			for (; it2 != short_a.end(); ++it2) {
 				newa.push_back(*it2);	
 			}
-			result -= G_eval(newa, scale);
+			result -= G_eval(newa, scale, gsyms);
 		}
 		return result / count_ones;
 	}
 
 	// G({1,...,1};y) -> G({1};y)^k / k!
 	if (all_ones && a.size() > 1) {
-		return pow(G_eval1(a.front(),scale), count_ones) / factorial(count_ones);
+		return pow(G_eval1(a.front(),scale, gsyms), count_ones) / factorial(count_ones);
 	}
 
 	// G({0,...,0};y) -> log(y)^k / k!
@@ -693,7 +692,7 @@ Gparameter prepare_pending_integrals(const Gparameter& pending_integrals, int sc
 
 
 // handles trailing zeroes for an otherwise convergent integral
-ex trailing_zeros_G(const Gparameter& a, int scale)
+ex trailing_zeros_G(const Gparameter& a, int scale, const exvector& gsyms)
 {
 	bool convergent;
 	int depth, trailing_zeros;
@@ -705,23 +704,23 @@ ex trailing_zeros_G(const Gparameter& a, int scale)
 	if ((trailing_zeros > 0) && (depth > 0)) {
 		ex result;
 		Gparameter new_a(a.begin(), a.end()-1);
-		result += G_eval1(0, scale) * trailing_zeros_G(new_a, scale);
+		result += G_eval1(0, scale, gsyms) * trailing_zeros_G(new_a, scale, gsyms);
 		for (Gparameter::const_iterator it = a.begin(); it != last; ++it) {
 			Gparameter new_a(a.begin(), it);
 			new_a.push_back(0);
 			new_a.insert(new_a.end(), it, a.end()-1);
-			result -= trailing_zeros_G(new_a, scale);
+			result -= trailing_zeros_G(new_a, scale, gsyms);
 		}
 
 		return result / trailing_zeros;
 	} else {
-		return G_eval(a, scale);
+		return G_eval(a, scale, gsyms);
 	}
 }
 
 
 // G transformation [VSW] (57),(58)
-ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, int scale)
+ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, int scale, const exvector& gsyms)
 {
 	// pendint = ( y1, b1, ..., br )
 	//       a = ( 0, ..., 0, amin )
@@ -752,15 +751,21 @@ ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, i
 			result -= I*Pi;
 		}
 		if (psize) {
-			result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front());
+			result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals),
+				                   pending_integrals.front(),
+						   gsyms);
 		}
 		
 		// G(y2_{-+}; sr)
-		result += trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front());
+		result += trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals),
+			                   new_pending_integrals.front(),
+					   gsyms);
 		
 		// G(0; sr)
 		new_pending_integrals.back() = 0;
-		result -= trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front());
+		result -= trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals),
+			                   new_pending_integrals.front(),
+					   gsyms);
 
 		return result;
 	}
@@ -772,14 +777,16 @@ ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, i
 	//term zeta_m
 	result -= zeta(a.size());
 	if (psize) {
-		result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front());
+		result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals),
+			                   pending_integrals.front(),
+					   gsyms);
 	}
 	
 	// term int_0^sr dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
 	//    = int_0^sr dt/t G_{m-1}( t_{+-}; y2 )
 	Gparameter new_a(a.begin()+1, a.end());
 	new_pending_integrals.push_back(0);
-	result -= depth_one_trafo_G(new_pending_integrals, new_a, scale);
+	result -= depth_one_trafo_G(new_pending_integrals, new_a, scale, gsyms);
 	
 	// term int_0^y2 dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
 	//    = int_0^y2 dt/t G_{m-1}( t_{+-}; y2 )
@@ -787,10 +794,12 @@ ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, i
 	new_pending_integrals_2.push_back(scale);
 	new_pending_integrals_2.push_back(0);
 	if (psize) {
-		result += trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front())
-		          * depth_one_trafo_G(new_pending_integrals_2, new_a, scale);
+		result += trailing_zeros_G(convert_pending_integrals_G(pending_integrals),
+			                   pending_integrals.front(),
+					   gsyms)
+		          * depth_one_trafo_G(new_pending_integrals_2, new_a, scale, gsyms);
 	} else {
-		result += depth_one_trafo_G(new_pending_integrals_2, new_a, scale);
+		result += depth_one_trafo_G(new_pending_integrals_2, new_a, scale, gsyms);
 	}
 
 	return result;
@@ -799,11 +808,13 @@ ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, i
 
 // forward declaration
 ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2,
-	     const Gparameter& pendint, const Gparameter& a_old, int scale);
+	     const Gparameter& pendint, const Gparameter& a_old, int scale,
+	     const exvector& gsyms);
 
 
 // G transformation [VSW]
-ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
+ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale,
+	       const exvector& gsyms)
 {
 	// main recursion routine
 	//
@@ -829,10 +840,12 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
 		if (a.size() == 0) {
 		  result = 1;
 		} else {
-		  result = G_eval(a, scale);
+		  result = G_eval(a, scale, gsyms);
 		}
 		if (pendint.size() > 0) {
-		  result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front());
+		  result *= trailing_zeros_G(convert_pending_integrals_G(pendint),
+			                     pendint.front(),
+					     gsyms);
 		} 
 		return result;
 	}
@@ -841,12 +854,12 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
 	if (trailing_zeros > 0) {
 		ex result;
 		Gparameter new_a(a.begin(), a.end()-1);
-		result += G_eval1(0, scale) * G_transform(pendint, new_a, scale);
+		result += G_eval1(0, scale, gsyms) * G_transform(pendint, new_a, scale, gsyms);
 		for (Gparameter::const_iterator it = a.begin(); it != firstzero; ++it) {
 			Gparameter new_a(a.begin(), it);
 			new_a.push_back(0);
 			new_a.insert(new_a.end(), it, a.end()-1);
-			result -= G_transform(pendint, new_a, scale);
+			result -= G_transform(pendint, new_a, scale, gsyms);
 		}
 		return result / trailing_zeros;
 	}
@@ -854,15 +867,17 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
 	// convergence case
 	if (convergent) {
 		if (pendint.size() > 0) {
-			return G_eval(convert_pending_integrals_G(pendint), pendint.front()) * G_eval(a, scale);
+			return G_eval(convert_pending_integrals_G(pendint),
+				      pendint.front(), gsyms)*
+				G_eval(a, scale, gsyms);
 		} else {
-			return G_eval(a, scale);
+			return G_eval(a, scale, gsyms);
 		}
 	}
 
 	// call basic transformation for depth equal one
 	if (depth == 1) {
-		return depth_one_trafo_G(pendint, a, scale);
+		return depth_one_trafo_G(pendint, a, scale, gsyms);
 	}
 
 	// do recursion
@@ -877,9 +892,10 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
 		Gparameter a1(a.begin(),min_it+1);
 		Gparameter a2(min_it+1,a.end());
 
-		ex result = G_transform(pendint,a2,scale)*G_transform(empty,a1,scale);
+		ex result = G_transform(pendint, a2, scale, gsyms)*
+			G_transform(empty, a1, scale, gsyms);
 
-		result -= shuffle_G(empty,a1,a2,pendint,a,scale);
+		result -= shuffle_G(empty, a1, a2, pendint, a, scale, gsyms);
 		return result;
 	}
 
@@ -890,9 +906,10 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
 	Gparameter new_pendint = prepare_pending_integrals(pendint, a[min_it_pos]);
 	Gparameter new_a = a;
 	new_a[min_it_pos] = 0;
-	ex result = G_transform(empty, new_a, scale);
+	ex result = G_transform(empty, new_a, scale, gsyms);
 	if (pendint.size() > 0) {
-		result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front());
+		result *= trailing_zeros_G(convert_pending_integrals_G(pendint),
+			                   pendint.front(), gsyms);
 	}
 
 	// other terms
@@ -901,29 +918,33 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
 	if (changeit != new_a.begin()) {
 		// smallest in the middle
 		new_pendint.push_back(*changeit);
-		result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
-			* G_transform(empty, new_a, scale);
+		result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint),
+			                   new_pendint.front(), gsyms)*
+			G_transform(empty, new_a, scale, gsyms);
 		int buffer = *changeit;
 		*changeit = *min_it;
-		result += G_transform(new_pendint, new_a, scale);
+		result += G_transform(new_pendint, new_a, scale, gsyms);
 		*changeit = buffer;
 		new_pendint.pop_back();
 		--changeit;
 		new_pendint.push_back(*changeit);
-		result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
-			* G_transform(empty, new_a, scale);
+		result += trailing_zeros_G(convert_pending_integrals_G(new_pendint),
+			                   new_pendint.front(), gsyms)*
+			G_transform(empty, new_a, scale, gsyms);
 		*changeit = *min_it;
-		result -= G_transform(new_pendint, new_a, scale);
+		result -= G_transform(new_pendint, new_a, scale, gsyms);
 	} else {
 		// smallest at the front
 		new_pendint.push_back(scale);
-		result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
-			* G_transform(empty, new_a, scale);
+		result += trailing_zeros_G(convert_pending_integrals_G(new_pendint),
+			                   new_pendint.front(), gsyms)*
+			G_transform(empty, new_a, scale, gsyms);
 		new_pendint.back() =  *changeit;
-		result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
-			* G_transform(empty, new_a, scale);
+		result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint),
+			                   new_pendint.front(), gsyms)*
+			G_transform(empty, new_a, scale, gsyms);
 		*changeit = *min_it;
-		result += G_transform(new_pendint, new_a, scale);
+		result += G_transform(new_pendint, new_a, scale, gsyms);
 	}
 	return result;
 }
@@ -932,27 +953,28 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
 // shuffles the two parameter list a1 and a2 and calls G_transform for every term except
 // for the one that is equal to a_old
 ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2,
-	     const Gparameter& pendint, const Gparameter& a_old, int scale) 
+	     const Gparameter& pendint, const Gparameter& a_old, int scale,
+	     const exvector& gsyms) 
 {
 	if (a1.size()==0 && a2.size()==0) {
 		// veto the one configuration we don't want
 		if ( a0 == a_old ) return 0;
 
-		return G_transform(pendint,a0,scale);
+		return G_transform(pendint, a0, scale, gsyms);
 	}
 
 	if (a2.size()==0) {
 		Gparameter empty;
 		Gparameter aa0 = a0;
 		aa0.insert(aa0.end(),a1.begin(),a1.end());
-		return shuffle_G(aa0,empty,empty,pendint,a_old,scale);
+		return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms);
 	}
 
 	if (a1.size()==0) {
 		Gparameter empty;
 		Gparameter aa0 = a0;
 		aa0.insert(aa0.end(),a2.begin(),a2.end());
-		return shuffle_G(aa0,empty,empty,pendint,a_old,scale);
+		return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms);
 	}
 
 	Gparameter a1_removed(a1.begin()+1,a1.end());
@@ -964,8 +986,8 @@ ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2
 	a01.push_back( a1[0] );
 	a02.push_back( a2[0] );
 
-	return shuffle_G(a01,a1_removed,a2,pendint,a_old,scale)
-	     + shuffle_G(a02,a1,a2_removed,pendint,a_old,scale);
+	return shuffle_G(a01, a1_removed, a2, pendint, a_old, scale, gsyms)
+	     + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms);
 }
 
 
@@ -980,9 +1002,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 +1013,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) {
 
@@ -1013,7 +1086,8 @@ ex G_numeric(const lst& x, const lst& s, const ex& y)
 
 		// generate missing dummy-symbols
 		int i = 1;
-		gsyms.clear();
+		// holding dummy-symbols for the G/Li transformations
+		exvector gsyms;
 		gsyms.push_back(symbol("GSYMS_ERROR"));
 		ex lastentry;
 		for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
@@ -1063,7 +1137,7 @@ ex G_numeric(const lst& x, const lst& s, const ex& y)
 
 		// do transformation
 		Gparameter pendint;
-		ex result = G_transform(pendint, a, scale);
+		ex result = G_transform(pendint, a, scale, gsyms);
 		// replace dummy symbols with their values
 		result = result.eval().expand();
 		result = result.subs(subslst).evalf();
@@ -1071,58 +1145,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;
@@ -1197,7 +1219,12 @@ static ex G2_evalf(const ex& x_, const ex& y)
 		if (*it != _ex0) {
 			all_zero = false;
 		}
-		s.append(+1);
+		if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
+			s.append(-1);
+		}
+		else {
+			s.append(+1);
+		}
 	}
 	if (all_zero) {
 		return pow(log(y), x.nops()) / factorial(x.nops());
@@ -1233,7 +1260,12 @@ static ex G2_eval(const ex& x_, const ex& y)
 		if (*it != _ex0) {
 			all_zero = false;
 		}
-		s.append(+1);
+		if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
+			s.append(-1);
+		}
+		else {
+			s.append(+1);
+		}
 	}
 	if (all_zero) {
 		return pow(log(y), x.nops()) / factorial(x.nops());
@@ -1286,10 +1318,21 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
 		if (*itx != _ex0) {
 			all_zero = false;
 		}
-		if (*its >= 0) {
-			sn.append(+1);
-		} else {
-			sn.append(-1);
+		if ( ex_to<numeric>(*itx).is_real() ) {
+			if ( *its >= 0 ) {
+				sn.append(+1);
+			}
+			else {
+				sn.append(-1);
+			}
+		}
+		else {
+			if ( ex_to<numeric>(*itx).imag() > 0 ) {
+				sn.append(+1);
+			}
+			else {
+				sn.append(-1);
+			}
 		}
 	}
 	if (all_zero) {
@@ -1333,10 +1376,21 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
 		if (*itx != _ex0) {
 			all_zero = false;
 		}
-		if (*its >= 0) {
-			sn.append(+1);
-		} else {
-			sn.append(-1);
+		if ( ex_to<numeric>(*itx).is_real() ) {
+			if ( *its >= 0 ) {
+				sn.append(+1);
+			}
+			else {
+				sn.append(-1);
+			}
+		}
+		else {
+			if ( ex_to<numeric>(*itx).imag() > 0 ) {
+				sn.append(+1);
+			}
+			else {
+				sn.append(-1);
+			}
 		}
 	}
 	if (all_zero) {
@@ -1376,12 +1430,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 +1555,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 +1567,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 +1779,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 +1790,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 +1868,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 +1929,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 +1940,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 +1976,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 +1993,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 +2010,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 +2046,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 +2082,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 +2099,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 +2564,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 +2690,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 +3104,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 +3112,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 +3153,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 +3167,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 +3176,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();
 	}
@@ -3132,13 +3417,6 @@ namespace {
 
 // parameters and data for [Cra] algorithm
 const cln::cl_N lambda = cln::cl_N("319/320");
-int L1;
-int L2;
-std::vector<std::vector<cln::cl_N> > f_kj;
-std::vector<cln::cl_N> crB;
-std::vector<std::vector<cln::cl_N> > crG;
-std::vector<cln::cl_N> crX;
-
 
 void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
 {
@@ -3153,44 +3431,39 @@ void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln
 
 
 // [Cra] section 4
-void initcX(const std::vector<int>& s)
+static void initcX(std::vector<cln::cl_N>& crX,
+	           const std::vector<int>& s,
+		   const int L2)
 {
-	const int k = s.size();
-
-	crX.clear();
-	crG.clear();
-	crB.clear();
-
-	for (int i=0; i<=L2; i++) {
-		crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
-	}
+	std::vector<cln::cl_N> crB(L2 + 1);
+	for (int i=0; i<=L2; i++)
+		crB[i] = bernoulli(i).to_cl_N() / cln::factorial(i);
 
 	int Sm = 0;
 	int Smp1 = 0;
-	for (int m=0; m<k-1; m++) {
-		std::vector<cln::cl_N> crGbuf;
-		Sm = Sm + s[m];
+	std::vector<std::vector<cln::cl_N> > crG(s.size() - 1, std::vector<cln::cl_N>(L2 + 1));
+	for (int m=0; m < s.size() - 1; m++) {
+		Sm += s[m];
 		Smp1 = Sm + s[m+1];
-		for (int i=0; i<=L2; i++) {
-			crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
-		}
-		crG.push_back(crGbuf);
+		for (int i = 0; i <= L2; i++)
+			crG[m][i] = cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2);
 	}
 
 	crX = crB;
 
-	for (int m=0; m<k-1; m++) {
-		std::vector<cln::cl_N> Xbuf;
-		for (int i=0; i<=L2; i++) {
-			Xbuf.push_back(crX[i] * crG[m][i]);
-		}
+	for (std::size_t m = 0; m < s.size() - 1; m++) {
+		std::vector<cln::cl_N> Xbuf(L2 + 1);
+		for (int i = 0; i <= L2; i++)
+			Xbuf[i] = crX[i] * crG[m][i];
+
 		halfcyclic_convolute(Xbuf, crB, crX);
 	}
 }
 
 
 // [Cra] section 4
-cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
+static cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk,
+	                         const std::vector<cln::cl_N>& crX)
 {
 	cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
 	cln::cl_N factor = cln::expt(lambda, Sqk);
@@ -3208,11 +3481,9 @@ cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
 
 
 // [Cra] section 4
-void calc_f(int maxr)
+static void calc_f(std::vector<std::vector<cln::cl_N> >& f_kj,
+	           const int maxr, const int L1)
 {
-	f_kj.clear();
-	f_kj.resize(L1);
-	
 	cln::cl_N t0, t1, t2, t3, t4;
 	int i, j, k;
 	std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
@@ -3238,7 +3509,8 @@ void calc_f(int maxr)
 
 
 // [Cra] (3.1)
-cln::cl_N crandall_Z(const std::vector<int>& s)
+static cln::cl_N crandall_Z(const std::vector<int>& s,
+	                    const std::vector<std::vector<cln::cl_N> >& f_kj)
 {
 	const int j = s.size();
 
@@ -3279,6 +3551,8 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
 	std::vector<int> r = s;
 	const int j = r.size();
 
+	std::size_t L1;
+
 	// decide on maximal size of f_kj for crandall_Z
 	if (Digits < 50) {
 		L1 = 150;
@@ -3286,6 +3560,7 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
 		L1 = Digits * 3 + j*2;
 	}
 
+	std::size_t L2;
 	// decide on maximal size of crX for crandall_Y
 	if (Digits < 38) {
 		L2 = 63;
@@ -3312,7 +3587,8 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
 		}
 	}
 
-	calc_f(maxr);
+	std::vector<std::vector<cln::cl_N> > f_kj(L1);
+	calc_f(f_kj, maxr, L1);
 
 	const cln::cl_N r0factorial = cln::factorial(r[0]-1);
 
@@ -3326,12 +3602,13 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
 		Srun -= skp1buf;
 		r.pop_back();
 
-		initcX(r);
+		std::vector<cln::cl_N> crX;
+		initcX(crX, r, L2);
 		
 		for (int q=0; q<skp1buf; q++) {
 			
-			cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
-			cln::cl_N pp2 = crandall_Z(rz);
+			cln::cl_N pp1 = crandall_Y_loop(Srun+q-k, crX);
+			cln::cl_N pp2 = crandall_Z(rz, f_kj);
 
 			rz.front()--;
 			
@@ -3345,9 +3622,11 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
 	}
 	rz.insert(rz.begin(), r.back());
 
-	initcX(rz);
+	std::vector<cln::cl_N> crX;
+	initcX(crX, rz, L2);
 
-	res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
+	res = (res + crandall_Y_loop(S-j, crX)) / r0factorial
+		+ crandall_Z(rz, f_kj);
 
 	return res;
 }
@@ -3537,18 +3816,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;
 				}