X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=4312377f10cf2d868972db41a56dac127973855d;hp=22e3125198b1ba12712c4b5af77d6c678b1e7481;hb=515171f0bcd42099c266713c3d605cd92cedd2e2;hpb=773570d668a8a86d4c646dfe9d8572746329652c

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index 22e31251..4312377f 100644
--- a/ginac/inifcns_nstdsums.cpp
+++ b/ginac/inifcns_nstdsums.cpp
@@ -47,7 +47,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2010 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
@@ -64,11 +64,6 @@
  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
-#include <sstream>
-#include <stdexcept>
-#include <vector>
-#include <cln/cln.h>
-
 #include "inifcns.h"
 
 #include "add.h"
@@ -84,6 +79,10 @@
 #include "utils.h"
 #include "wildcard.h"
 
+#include <cln/cln.h>
+#include <sstream>
+#include <stdexcept>
+#include <vector>
 
 namespace GiNaC {
 
@@ -210,7 +209,7 @@ void double_Xn()
 			}
 		}
 		// X_n
-		for (int n=2; n<Xn.size(); ++n) {
+		for (size_t n=2; n<Xn.size(); ++n) {
 			for (int i=xninitsize+1; i<=xend; ++i) {
 				if (i & 1) {
 					result = 0; // k == 0
@@ -501,19 +500,6 @@ cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl
 }
 
 
-// converts parameter types and calls multipleLi_do_sum (convenience function for G_numeric)
-cln::cl_N mLi_do_summation(const lst& m, const lst& x)
-{
-	std::vector<int> m_int;
-	std::vector<cln::cl_N> x_cln;
-	for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
-		m_int.push_back(ex_to<numeric>(*itm).to_int());
-		x_cln.push_back(ex_to<numeric>(*itx).to_cl_N());
-	}
-	return multipleLi_do_sum(m_int, x_cln);
-}
-
-
 // forward declaration for Li_eval()
 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf);
 
@@ -583,11 +569,12 @@ ex G_eval(const Gparameter& a, int scale, const exvector& gsyms)
 
 			Gparameter newa;
 			Gparameter::const_iterator it2 = short_a.begin();
-			for (--it2; it2 != it;) {
-				++it2;
+			for (; it2 != it; ++it2) {
 				newa.push_back(*it2);
 			}
+			newa.push_back(*it);
 			newa.push_back(a[0]);
+			it2 = it;
 			++it2;
 			for (; it2 != short_a.end(); ++it2) {
 				newa.push_back(*it2);	
@@ -672,6 +659,8 @@ Gparameter::const_iterator check_parameter_G(const Gparameter& a, int scale,
 			++trailing_zeros;
 		}
 	}
+	if (lastnonzero == a.end())
+		return a.end();
 	return ++lastnonzero;
 }
 
@@ -990,199 +979,228 @@ ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2
 	     + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms);
 }
 
-
 // handles the transformations and the numerical evaluation of G
 // the parameter x, s and y must only contain numerics
-ex G_numeric(const lst& x, const lst& s, const ex& y)
+static cln::cl_N
+G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
+	  const cln::cl_N& y);
+
+// do acceleration transformation (hoelder convolution [BBB])
+// the parameter x, s and y must only contain numerics
+static cln::cl_N
+G_do_hoelder(std::vector<cln::cl_N> x, /* yes, it's passed by value */
+	     const std::vector<int>& s, const cln::cl_N& y)
 {
-	// check for convergence and necessary accelerations
-	bool need_trafo = false;
-	bool need_hoelder = false;
-	int depth = 0;
-	for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
-		if (!(*it).is_zero()) {
-			++depth;
-			if (abs(*it) - y < -pow(10,-Digits+1)) {
-				need_trafo = true;
+	cln::cl_N result;
+	const std::size_t size = x.size();
+	for (std::size_t i = 0; i < size; ++i)
+		x[i] = x[i]/y;
+
+	for (std::size_t r = 0; r <= size; ++r) {
+		cln::cl_N buffer(1 & r ? -1 : 1);
+		cln::cl_RA p(2);
+		bool adjustp;
+		do {
+			adjustp = false;
+			for (std::size_t i = 0; i < size; ++i) {
+				if (x[i] == cln::cl_RA(1)/p) {
+					p = p/2 + cln::cl_RA(3)/2;
+					adjustp = true;
+					continue;
+				}
 			}
-			if (abs((abs(*it) - y)/y) < 0.01) {
-				need_hoelder = true;
+		} while (adjustp);
+		cln::cl_RA q = p/(p-1);
+		std::vector<cln::cl_N> qlstx;
+		std::vector<int> qlsts;
+		for (std::size_t j = r; j >= 1; --j) {
+			qlstx.push_back(cln::cl_N(1) - x[j-1]);
+			if (instanceof(x[j-1], cln::cl_R_ring) &&
+			    realpart(x[j-1]) > 1 && realpart(x[j-1]) <= 2) {
+				qlsts.push_back(s[j-1]);
+			} else {
+				qlsts.push_back(-s[j-1]);
 			}
 		}
-	}
-	if (x.op(x.nops()-1).is_zero()) {
-		need_trafo = true;
-	}
-	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);
+		if (qlstx.size() > 0) {
+			buffer = buffer*G_numeric(qlstx, qlsts, 1/q);
 		}
-		
-		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;
+		std::vector<cln::cl_N> plstx;
+		std::vector<int> plsts;
+		for (std::size_t j = r+1; j <= size; ++j) {
+			plstx.push_back(x[j-1]);
+			plsts.push_back(s[j-1]);
 		}
-		return result;
+		if (plstx.size() > 0) {
+			buffer = buffer*G_numeric(plstx, plsts, 1/p);
+		}
+		result = result + buffer;
 	}
-	
-	// convergence transformation
-	if (need_trafo) {
-
-		// sort (|x|<->position) to determine indices
-		std::multimap<ex,int> sortmap;
-		int size = 0;
-		for (int i=0; i<x.nops(); ++i) {
-			if (!x[i].is_zero()) {
-				sortmap.insert(std::pair<ex,int>(abs(x[i]), i));
-				++size;
-			}
-		}
-		// include upper limit (scale)
-		sortmap.insert(std::pair<ex,int>(abs(y), x.nops()));
-
-		// generate missing dummy-symbols
-		int i = 1;
-		// 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) {
-			if (it != sortmap.begin()) {
-				if (it->second < x.nops()) {
-					if (x[it->second] == lastentry) {
-						gsyms.push_back(gsyms.back());
-						continue;
-					}
-				} else {
-					if (y == lastentry) {
-						gsyms.push_back(gsyms.back());
-						continue;
-					}
+	return result;
+}
+
+// convergence transformation, used for numerical evaluation of G function.
+// the parameter x, s and y must only contain numerics
+static cln::cl_N
+G_do_trafo(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
+	   const cln::cl_N& y)
+{
+	// sort (|x|<->position) to determine indices
+	typedef std::multimap<cln::cl_R, std::size_t> sortmap_t;
+	sortmap_t sortmap;
+	std::size_t size = 0;
+	for (std::size_t i = 0; i < x.size(); ++i) {
+		if (!zerop(x[i])) {
+			sortmap.insert(std::make_pair(abs(x[i]), i));
+			++size;
+		}
+	}
+	// include upper limit (scale)
+	sortmap.insert(std::make_pair(abs(y), x.size()));
+
+	// generate missing dummy-symbols
+	int i = 1;
+	// holding dummy-symbols for the G/Li transformations
+	exvector gsyms;
+	gsyms.push_back(symbol("GSYMS_ERROR"));
+	cln::cl_N lastentry(0);
+	for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
+		if (it != sortmap.begin()) {
+			if (it->second < x.size()) {
+				if (x[it->second] == lastentry) {
+					gsyms.push_back(gsyms.back());
+					continue;
 				}
-			}
-			std::ostringstream os;
-			os << "a" << i;
-			gsyms.push_back(symbol(os.str()));
-			++i;
-			if (it->second < x.nops()) {
-				lastentry = x[it->second];
 			} else {
-				lastentry = y;
+				if (y == lastentry) {
+					gsyms.push_back(gsyms.back());
+					continue;
+				}
 			}
 		}
+		std::ostringstream os;
+		os << "a" << i;
+		gsyms.push_back(symbol(os.str()));
+		++i;
+		if (it->second < x.size()) {
+			lastentry = x[it->second];
+		} else {
+			lastentry = y;
+		}
+	}
 
-		// fill position data according to sorted indices and prepare substitution list
-		Gparameter a(x.nops());
-		lst subslst;
-		int pos = 1;
-		int scale;
-		for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
-			if (it->second < x.nops()) {
-				if (s[it->second] > 0) {
-					a[it->second] = pos;
-				} else {
-					a[it->second] = -pos;
-				}
-				subslst.append(gsyms[pos] == x[it->second]);
+	// fill position data according to sorted indices and prepare substitution list
+	Gparameter a(x.size());
+	exmap subslst;
+	std::size_t pos = 1;
+	int scale = pos;
+	for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
+		if (it->second < x.size()) {
+			if (s[it->second] > 0) {
+				a[it->second] = pos;
 			} else {
-				scale = pos;
-				subslst.append(gsyms[pos] == y);
+				a[it->second] = -int(pos);
 			}
-			++pos;
+			subslst[gsyms[pos]] = numeric(x[it->second]);
+		} else {
+			scale = pos;
+			subslst[gsyms[pos]] = numeric(y);
 		}
+		++pos;
+	}
 
-		// do transformation
-		Gparameter pendint;
-		ex result = G_transform(pendint, a, scale, gsyms);
-		// replace dummy symbols with their values
-		result = result.eval().expand();
-		result = result.subs(subslst).evalf();
-		
-		return result;
+	// do transformation
+	Gparameter pendint;
+	ex result = G_transform(pendint, a, scale, gsyms);
+	// replace dummy symbols with their values
+	result = result.eval().expand();
+	result = result.subs(subslst).evalf();
+	if (!is_a<numeric>(result))
+		throw std::logic_error("G_do_trafo: G_transform returned non-numeric result");
+	
+	cln::cl_N ret = ex_to<numeric>(result).to_cl_N();
+	return ret;
+}
+
+// handles the transformations and the numerical evaluation of G
+// the parameter x, s and y must only contain numerics
+static cln::cl_N
+G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
+	  const cln::cl_N& y)
+{
+	// check for convergence and necessary accelerations
+	bool need_trafo = false;
+	bool need_hoelder = false;
+	std::size_t depth = 0;
+	for (std::size_t i = 0; i < x.size(); ++i) {
+		if (!zerop(x[i])) {
+			++depth;
+			const cln::cl_N x_y = abs(x[i]) - y;
+			if (instanceof(x_y, cln::cl_R_ring) &&
+			    realpart(x_y) < cln::least_negative_float(cln::float_format(Digits - 2)))
+				need_trafo = true;
+
+			if (abs(abs(x[i]/y) - 1) < 0.01)
+				need_hoelder = true;
+		}
 	}
+	if (zerop(x[x.size() - 1]))
+		need_trafo = true;
+
+	if (depth == 1 && x.size() == 2 && !need_trafo)
+		return - Li_projection(2, y/x[1], cln::float_format(Digits));
+	
+	// do acceleration transformation (hoelder convolution [BBB])
+	if (need_hoelder)
+		return G_do_hoelder(x, s, y);
+	
+	// convergence transformation
+	if (need_trafo)
+		return G_do_trafo(x, s, y);
 
 	// do summation
-	lst newx;
-	lst m;
+	std::vector<cln::cl_N> newx;
+	newx.reserve(x.size());
+	std::vector<int> m;
+	m.reserve(x.size());
 	int mcount = 1;
-	ex sign = 1;
-	ex factor = y;
-	for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
-		if ((*it).is_zero()) {
+	int sign = 1;
+	cln::cl_N factor = y;
+	for (std::size_t i = 0; i < x.size(); ++i) {
+		if (zerop(x[i])) {
 			++mcount;
 		} else {
-			newx.append(factor / (*it));
-			factor = *it;
-			m.append(mcount);
+			newx.push_back(factor/x[i]);
+			factor = x[i];
+			m.push_back(mcount);
 			mcount = 1;
 			sign = -sign;
 		}
 	}
 
-	return sign * numeric(mLi_do_summation(m, newx));
+	return sign*multipleLi_do_sum(m, newx);
 }
 
 
 ex mLi_numeric(const lst& m, const lst& x)
 {
 	// let G_numeric do the transformation
-	lst newx;
-	lst s;
-	ex factor = 1;
+	std::vector<cln::cl_N> newx;
+	newx.reserve(x.nops());
+	std::vector<int> s;
+	s.reserve(x.nops());
+	cln::cl_N factor(1);
 	for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
 		for (int i = 1; i < *itm; ++i) {
-			newx.append(0);
-			s.append(1);
+			newx.push_back(cln::cl_N(0));
+			s.push_back(1);
 		}
-		newx.append(factor / *itx);
-		factor /= *itx;
-		s.append(1);
+		const cln::cl_N xi = ex_to<numeric>(*itx).to_cl_N();
+		newx.push_back(factor/xi);
+		factor = factor/xi;
+		s.push_back(1);
 	}
-	return pow(-1, m.nops()) * G_numeric(newx, s, _ex1);
+	return numeric(cln::cl_N(1 & m.nops() ? - 1 : 1)*G_numeric(newx, s, cln::cl_N(1)));
 }
 
 
@@ -1210,7 +1228,8 @@ static ex G2_evalf(const ex& x_, const ex& y)
 	if (x.op(0) == y) {
 		return G(x_, y).hold();
 	}
-	lst s;
+	std::vector<int> s;
+	s.reserve(x.nops());
 	bool all_zero = true;
 	for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
 		if (!(*it).info(info_flags::numeric)) {
@@ -1220,16 +1239,21 @@ static ex G2_evalf(const ex& x_, const ex& y)
 			all_zero = false;
 		}
 		if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
-			s.append(-1);
+			s.push_back(-1);
 		}
 		else {
-			s.append(+1);
+			s.push_back(1);
 		}
 	}
 	if (all_zero) {
 		return pow(log(y), x.nops()) / factorial(x.nops());
 	}
-	return G_numeric(x, s, y);
+	std::vector<cln::cl_N> xv;
+	xv.reserve(x.nops());
+	for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+		xv.push_back(ex_to<numeric>(*it).to_cl_N());
+	cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
+	return numeric(result);
 }
 
 
@@ -1247,7 +1271,8 @@ static ex G2_eval(const ex& x_, const ex& y)
 	if (x.op(0) == y) {
 		return G(x_, y).hold();
 	}
-	lst s;
+	std::vector<int> s;
+	s.reserve(x.nops());
 	bool all_zero = true;
 	bool crational = true;
 	for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
@@ -1261,10 +1286,10 @@ static ex G2_eval(const ex& x_, const ex& y)
 			all_zero = false;
 		}
 		if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
-			s.append(-1);
+			s.push_back(-1);
 		}
 		else {
-			s.append(+1);
+			s.push_back(+1);
 		}
 	}
 	if (all_zero) {
@@ -1276,7 +1301,12 @@ static ex G2_eval(const ex& x_, const ex& y)
 	if (crational) {
 		return G(x_, y).hold();
 	}
-	return G_numeric(x, s, y);
+	std::vector<cln::cl_N> xv;
+	xv.reserve(x.nops());
+	for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+		xv.push_back(ex_to<numeric>(*it).to_cl_N());
+	cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
+	return numeric(result);
 }
 
 
@@ -1306,7 +1336,8 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
 	if (x.op(0) == y) {
 		return G(x_, s_, y).hold();
 	}
-	lst sn;
+	std::vector<int> sn;
+	sn.reserve(s.nops());
 	bool all_zero = true;
 	for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
 		if (!(*itx).info(info_flags::numeric)) {
@@ -1320,25 +1351,30 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
 		}
 		if ( ex_to<numeric>(*itx).is_real() ) {
 			if ( *its >= 0 ) {
-				sn.append(+1);
+				sn.push_back(1);
 			}
 			else {
-				sn.append(-1);
+				sn.push_back(-1);
 			}
 		}
 		else {
 			if ( ex_to<numeric>(*itx).imag() > 0 ) {
-				sn.append(+1);
+				sn.push_back(1);
 			}
 			else {
-				sn.append(-1);
+				sn.push_back(-1);
 			}
 		}
 	}
 	if (all_zero) {
 		return pow(log(y), x.nops()) / factorial(x.nops());
 	}
-	return G_numeric(x, sn, y);
+	std::vector<cln::cl_N> xn;
+	xn.reserve(x.nops());
+	for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+		xn.push_back(ex_to<numeric>(*it).to_cl_N());
+	cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
+	return numeric(result);
 }
 
 
@@ -1360,7 +1396,8 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
 	if (x.op(0) == y) {
 		return G(x_, s_, y).hold();
 	}
-	lst sn;
+	std::vector<int> sn;
+	sn.reserve(s.nops());
 	bool all_zero = true;
 	bool crational = true;
 	for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
@@ -1378,18 +1415,18 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
 		}
 		if ( ex_to<numeric>(*itx).is_real() ) {
 			if ( *its >= 0 ) {
-				sn.append(+1);
+				sn.push_back(1);
 			}
 			else {
-				sn.append(-1);
+				sn.push_back(-1);
 			}
 		}
 		else {
 			if ( ex_to<numeric>(*itx).imag() > 0 ) {
-				sn.append(+1);
+				sn.push_back(1);
 			}
 			else {
-				sn.append(-1);
+				sn.push_back(-1);
 			}
 		}
 	}
@@ -1402,7 +1439,12 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
 	if (crational) {
 		return G(x_, s_, y).hold();
 	}
-	return G_numeric(x, sn, y);
+	std::vector<cln::cl_N> xn;
+	xn.reserve(x.nops());
+	for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+		xn.push_back(ex_to<numeric>(*it).to_cl_N());
+	cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
+	return numeric(result);
 }
 
 
@@ -1644,7 +1686,7 @@ static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c)
 	} else {
 		x = lst(x_);
 	}
-	c.s << "\\mbox{Li}_{";
+	c.s << "\\mathrm{Li}_{";
 	lst::const_iterator itm = m.begin();
 	(*itm).print(c);
 	itm++;
@@ -2160,7 +2202,7 @@ static ex S_deriv(const ex& n, const ex& p, const ex& x, unsigned deriv_param)
 
 static void S_print_latex(const ex& n, const ex& p, const ex& x, const print_context& c)
 {
-	c.s << "\\mbox{S}_{";
+	c.s << "\\mathrm{S}_{";
 	n.print(c);
 	c.s << ",";
 	p.print(c);
@@ -2241,7 +2283,7 @@ bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
 		}
 	}
 	if (has_negative_parameters) {
-		for (int i=0; i<m.nops(); i++) {
+		for (std::size_t i=0; i<m.nops(); i++) {
 			if (m.op(i) < 0) {
 				m.let_op(i) = -m.op(i);
 				s.append(-1);
@@ -2281,7 +2323,7 @@ struct map_trafo_H_convert_to_Li : public map_function
 					s.let_op(0) = s.op(0) * arg;
 					return pf * Li(m, s).hold();
 				} else {
-					for (int i=0; i<m.nops(); i++) {
+					for (std::size_t i=0; i<m.nops(); i++) {
 						s.append(1);
 					}
 					s.let_op(0) = s.op(0) * arg;
@@ -2363,7 +2405,7 @@ struct map_trafo_H_reduce_trailing_zeros : public map_function
 					
 					//
 					parameter.remove_last();
-					int lastentry = parameter.nops();
+					std::size_t lastentry = parameter.nops();
 					while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
 						lastentry--;
 					}
@@ -2452,9 +2494,9 @@ ex trafo_H_mult(const ex& h1, const ex& h2)
 			hlong = h2.op(0).op(0);
 		}
 	}
-	for (int i=0; i<=hlong.nops(); i++) {
+	for (std::size_t i=0; i<=hlong.nops(); i++) {
 		lst newparameter;
-		int j=0;
+		std::size_t j=0;
 		for (; j<i; j++) {
 			newparameter.append(hlong[j]);
 		}
@@ -2482,7 +2524,7 @@ struct map_trafo_H_mult : public map_function
 			ex result = 1;
 			ex firstH;
 			lst Hlst;
-			for (int pos=0; pos<e.nops(); pos++) {
+			for (std::size_t pos=0; pos<e.nops(); pos++) {
 				if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
 					std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
 					if (name == "H") {
@@ -2516,7 +2558,7 @@ struct map_trafo_H_mult : public map_function
 			if (Hlst.nops() > 0) {
 				ex buffer = trafo_H_mult(firstH, Hlst.op(0));
 				result *= buffer;
-				for (int i=1; i<Hlst.nops(); i++) {
+				for (std::size_t i=1; i<Hlst.nops(); i++) {
 					result *= Hlst.op(i);
 				}
 				result = result.expand();
@@ -2544,7 +2586,7 @@ ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
 	if (name == "H") {
 		h = e;
 	} else {
-		for (int i=0; i<e.nops(); i++) {
+		for (std::size_t 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") {
@@ -2576,7 +2618,7 @@ ex trafo_H_prepend_one(const ex& e, const ex& arg)
 	if (name == "H") {
 		h = e;
 	} else {
-		for (int i=0; i<e.nops(); i++) {
+		for (std::size_t 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") {
@@ -2607,7 +2649,7 @@ ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
 	if (name == "H") {
 		h = e;
 	} else {
-		for (int i=0; i<e.nops(); i++) {
+		for (std::size_t 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") {
@@ -2640,7 +2682,7 @@ ex trafo_H_1mxt1px_prepend_minusone(const ex& e, const ex& arg)
 	if (name == "H") {
 		h = e;
 	} else {
-		for (int i=0; i<e.nops(); i++) {
+		for (std::size_t 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") {
@@ -2671,7 +2713,7 @@ ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
 	if (name == "H") {
 		h = e;
 	} else {
-		for (int i=0; i<e.nops(); i++) {
+		for (std::size_t 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") {
@@ -2709,7 +2751,7 @@ struct map_trafo_H_1mx : public map_function
 				// 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++) {
+					for (std::size_t i = 1; i < parameter.nops(); i++) {
 						if (parameter.op(i) != 0) {
 							allthesame = false;
 							break;
@@ -2725,7 +2767,7 @@ struct map_trafo_H_1mx : public map_function
 				} 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++) {
+					for (std::size_t i = 1; i < parameter.nops(); i++) {
 						if (parameter.op(i) != 1) {
 							allthesame = false;
 							break;
@@ -2751,7 +2793,7 @@ struct map_trafo_H_1mx : public map_function
 					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++) {
+						for (std::size_t i = 0; i < buffer.nops(); i++) {
 							res -= trafo_H_prepend_one(buffer.op(i), arg);
 						}
 					} else {
@@ -2765,13 +2807,13 @@ struct map_trafo_H_1mx : public map_function
 					map_trafo_H_1mx recursion;
 					map_trafo_H_mult unify;
 					ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
-					int firstzero = 0;
+					std::size_t firstzero = 0;
 					while (parameter.op(firstzero) == 1) {
 						firstzero++;
 					}
-					for (int i=firstzero-1; i<parameter.nops()-1; i++) {
+					for (std::size_t i = firstzero-1; i < parameter.nops()-1; i++) {
 						lst newparameter;
-						int j=0;
+						std::size_t j=0;
 						for (; j<=i; j++) {
 							newparameter.append(parameter[j+1]);
 						}
@@ -2810,7 +2852,7 @@ struct map_trafo_H_1overx : public map_function
 				// 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++) {
+					for (std::size_t i = 1; i < parameter.nops(); i++) {
 						if (parameter.op(i) != 0) {
 							allthesame = false;
 							break;
@@ -2820,7 +2862,7 @@ struct map_trafo_H_1overx : public map_function
 						return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
 					}
 				} else if (parameter.op(0) == -1) {
-					for (int i=1; i<parameter.nops(); i++) {
+					for (std::size_t i = 1; i < parameter.nops(); i++) {
 						if (parameter.op(i) != -1) {
 							allthesame = false;
 							break;
@@ -2832,7 +2874,7 @@ struct map_trafo_H_1overx : public map_function
 						       / factorial(parameter.nops())).expand());
 					}
 				} else {
-					for (int i=1; i<parameter.nops(); i++) {
+					for (std::size_t i = 1; i < parameter.nops(); i++) {
 						if (parameter.op(i) != 1) {
 							allthesame = false;
 							break;
@@ -2855,7 +2897,7 @@ struct map_trafo_H_1overx : public map_function
 					map_trafo_H_1overx recursion;
 					ex buffer = recursion(H(newparameter, arg).hold());
 					if (is_a<add>(buffer)) {
-						for (int i=0; i<buffer.nops(); i++) {
+						for (std::size_t i = 0; i < buffer.nops(); i++) {
 							res += trafo_H_1tx_prepend_zero(buffer.op(i), arg);
 						}
 					} else {
@@ -2870,7 +2912,7 @@ struct map_trafo_H_1overx : public map_function
 					map_trafo_H_1overx recursion;
 					ex buffer = recursion(H(newparameter, arg).hold());
 					if (is_a<add>(buffer)) {
-						for (int i=0; i<buffer.nops(); i++) {
+						for (std::size_t i = 0; i < buffer.nops(); i++) {
 							res += trafo_H_1tx_prepend_zero(buffer.op(i), arg) - trafo_H_1tx_prepend_minusone(buffer.op(i), arg);
 						}
 					} else {
@@ -2884,13 +2926,13 @@ struct map_trafo_H_1overx : public map_function
 					map_trafo_H_1overx recursion;
 					map_trafo_H_mult unify;
 					ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
-					int firstzero = 0;
+					std::size_t firstzero = 0;
 					while (parameter.op(firstzero) == 1) {
 						firstzero++;
 					}
-					for (int i=firstzero-1; i<parameter.nops()-1; i++) {
+					for (std::size_t i = firstzero-1; i < parameter.nops() - 1; i++) {
 						lst newparameter;
-						int j=0;
+						std::size_t j = 0;
 						for (; j<=i; j++) {
 							newparameter.append(parameter[j+1]);
 						}
@@ -2931,7 +2973,7 @@ struct map_trafo_H_1mxt1px : public map_function
 				// 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++) {
+					for (std::size_t i = 1; i < parameter.nops(); i++) {
 						if (parameter.op(i) != 0) {
 							allthesame = false;
 							break;
@@ -2943,7 +2985,7 @@ struct map_trafo_H_1mxt1px : public map_function
 						       / factorial(parameter.nops())).expand());
 					}
 				} else if (parameter.op(0) == -1) {
-					for (int i=1; i<parameter.nops(); i++) {
+					for (std::size_t i = 1; i < parameter.nops(); i++) {
 						if (parameter.op(i) != -1) {
 							allthesame = false;
 							break;
@@ -2955,7 +2997,7 @@ struct map_trafo_H_1mxt1px : public map_function
 						       / factorial(parameter.nops())).expand());
 					}
 				} else {
-					for (int i=1; i<parameter.nops(); i++) {
+					for (std::size_t i = 1; i < parameter.nops(); i++) {
 						if (parameter.op(i) != 1) {
 							allthesame = false;
 							break;
@@ -2978,7 +3020,7 @@ struct map_trafo_H_1mxt1px : public map_function
 					map_trafo_H_1mxt1px recursion;
 					ex buffer = recursion(H(newparameter, arg).hold());
 					if (is_a<add>(buffer)) {
-						for (int i=0; i<buffer.nops(); i++) {
+						for (std::size_t i = 0; i < buffer.nops(); i++) {
 							res -= trafo_H_1mxt1px_prepend_one(buffer.op(i), arg) + trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
 						}
 					} else {
@@ -2993,7 +3035,7 @@ struct map_trafo_H_1mxt1px : public map_function
 					map_trafo_H_1mxt1px recursion;
 					ex buffer = recursion(H(newparameter, arg).hold());
 					if (is_a<add>(buffer)) {
-						for (int i=0; i<buffer.nops(); i++) {
+						for (std::size_t i = 0; i < buffer.nops(); i++) {
 							res -= trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
 						}
 					} else {
@@ -3007,13 +3049,13 @@ struct map_trafo_H_1mxt1px : public map_function
 					map_trafo_H_1mxt1px recursion;
 					map_trafo_H_mult unify;
 					ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
-					int firstzero = 0;
+					std::size_t firstzero = 0;
 					while (parameter.op(firstzero) == 1) {
 						firstzero++;
 					}
-					for (int i=firstzero-1; i<parameter.nops()-1; i++) {
+					for (std::size_t i = firstzero - 1; i < parameter.nops() - 1; i++) {
 						lst newparameter;
-						int j=0;
+						std::size_t j=0;
 						for (; j<=i; j++) {
 							newparameter.append(parameter[j+1]);
 						}
@@ -3087,7 +3129,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
 			}
 		}
 
-		for (int i=0; i<x1.nops(); i++) {
+		for (std::size_t i = 0; i < x1.nops(); i++) {
 			if (!x1.op(i).info(info_flags::integer)) {
 				return H(x1, x2).hold();
 			}
@@ -3159,7 +3201,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
 		// ensure that the realpart of the argument is positive
 		if (cln::realpart(x) < 0) {
 			x = -x;
-			for (int i=0; i<m.nops(); i++) {
+			for (std::size_t i = 0; i < m.nops(); i++) {
 				if (m.op(i) != 0) {
 					m.let_op(i) = -m.op(i);
 					res *= -1;
@@ -3366,7 +3408,7 @@ static void H_print_latex(const ex& m_, const ex& x, const print_context& c)
 	} else {
 		m = lst(m_);
 	}
-	c.s << "\\mbox{H}_{";
+	c.s << "\\mathrm{H}_{";
 	lst::const_iterator itm = m.begin();
 	(*itm).print(c);
 	itm++;
@@ -3417,11 +3459,6 @@ namespace {
 
 // parameters and data for [Cra] algorithm
 const cln::cl_N lambda = cln::cl_N("319/320");
-int L2;
-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)
 {
@@ -3436,44 +3473,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 < (int)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);
@@ -3570,6 +3602,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;
@@ -3611,11 +3644,12 @@ 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 pp1 = crandall_Y_loop(Srun+q-k, crX);
 			cln::cl_N pp2 = crandall_Z(rz, f_kj);
 
 			rz.front()--;
@@ -3630,9 +3664,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, f_kj);
+	res = (res + crandall_Y_loop(S-j, crX)) / r0factorial
+		+ crandall_Z(rz, f_kj);
 
 	return res;
 }
@@ -3674,7 +3710,7 @@ cln::cl_N zeta_do_Hoelder_convolution(const std::vector<int>& m_, const std::vec
 	s_p[0] = s_p[0] * cln::cl_N("1/2");
 	// convert notations
 	int sig = 1;
-	for (int i=0; i<s_.size(); i++) {
+	for (std::size_t i = 0; i < s_.size(); i++) {
 		if (s_[i] < 0) {
 			sig = -sig;
 			s_p[i] = -s_p[i];