X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=bf5ce093165cb7b30d15f410f7f72e5f23582d2d;hb=bb6b3d82cdf9e7ff4ecac89c47e63024e39ec96b;hp=5c5411b5700af4bf5cae92e6ecd58c5be4b6f991;hpb=cfc58a216ac4a3b18e4e89f7256c2b87ae16b1e4;p=ginac.git

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index 5c5411b5..bf5ce093 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
@@ -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;
@@ -488,11 +486,16 @@ cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl
 		t0buf = t[0];
 		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]);
+		}
+		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) || flag_accidental_zero );
+	} 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();
@@ -1716,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);
 					}
 				}
 			}
@@ -1727,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));
 					}
 				}
@@ -1805,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++) {
@@ -1856,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;
@@ -1867,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) {
@@ -1903,11 +1976,11 @@ 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)
@@ -1920,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;
@@ -1937,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);
@@ -1973,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);
 			}
 		}
 	}
@@ -2003,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()) {
@@ -3098,7 +3181,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;