Fix convergence of Li_projection at ±I.
[ginac.git] / ginac / inifcns_nstdsums.cpp
index 88f0fa4ee5f313d3a3c2925294616984782c94c4..67b06bae54f605599aa09d5acbb06352d9f714c5 100644 (file)
@@ -4,12 +4,12 @@
  *
  *  The functions are:
  *    classical polylogarithm              Li(n,x)
- *    multiple polylogarithm               Li(lst(m_1,...,m_k),lst(x_1,...,x_k))
- *                                         G(lst(a_1,...,a_k),y) or G(lst(a_1,...,a_k),lst(s_1,...,s_k),y)
+ *    multiple polylogarithm               Li(lst{m_1,...,m_k},lst{x_1,...,x_k})
+ *                                         G(lst{a_1,...,a_k},y) or G(lst{a_1,...,a_k},lst{s_1,...,s_k},y)
  *    Nielsen's generalized polylogarithm  S(n,p,x)
- *    harmonic polylogarithm               H(m,x) or H(lst(m_1,...,m_k),x)
- *    multiple zeta value                  zeta(m) or zeta(lst(m_1,...,m_k))
- *    alternating Euler sum                zeta(m,s) or zeta(lst(m_1,...,m_k),lst(s_1,...,s_k))
+ *    harmonic polylogarithm               H(m,x) or H(lst{m_1,...,m_k},x)
+ *    multiple zeta value                  zeta(m) or zeta(lst{m_1,...,m_k})
+ *    alternating Euler sum                zeta(m,s) or zeta(lst{m_1,...,m_k},lst{s_1,...,s_k})
  *
  *  Some remarks:
  *
@@ -25,7 +25,7 @@
  *      0, 1 and -1 --- or in compactified --- a string with zeros in front of 1 or -1 is written as a single
  *      number --- notation.
  *
- *    - All functions can be nummerically evaluated with arguments in the whole complex plane. The parameters
+ *    - All functions can be numerically evaluated with arguments in the whole complex plane. The parameters
  *      for Li, zeta and S must be positive integers. If you want to have an alternating Euler sum, you have
  *      to give the signs of the parameters as a second argument s to zeta(m,s) containing 1 and -1.
  *
@@ -47,7 +47,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2016 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
  *  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"
 #include "utils.h"
 #include "wildcard.h"
 
+#include <cln/cln.h>
+#include <sstream>
+#include <stdexcept>
+#include <vector>
 
 namespace GiNaC {
 
@@ -103,7 +102,7 @@ namespace {
 
 // lookup table for factors built from Bernoulli numbers
 // see fill_Xn()
-std::vector<std::vector<cln::cl_N> > Xn;
+std::vector<std::vector<cln::cl_N>> Xn;
 // initial size of Xn that should suffice for 32bit machines (must be even)
 const int xninitsizestep = 26;
 int xninitsize = xninitsizestep;
@@ -125,7 +124,7 @@ void fill_Xn(int n)
        if (n>1) {
                // calculate X_2 and higher (corresponding to Li_4 and higher)
                std::vector<cln::cl_N> buf(xninitsize);
-               std::vector<cln::cl_N>::iterator it = buf.begin();
+               auto it = buf.begin();
                cln::cl_N result;
                *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
                it++;
@@ -150,7 +149,7 @@ void fill_Xn(int n)
        } else if (n==1) {
                // special case to handle the X_0 correct
                std::vector<cln::cl_N> buf(xninitsize);
-               std::vector<cln::cl_N>::iterator it = buf.begin();
+               auto it = buf.begin();
                cln::cl_N result;
                *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
                it++;
@@ -174,7 +173,7 @@ void fill_Xn(int n)
        } else {
                // calculate X_0
                std::vector<cln::cl_N> buf(xninitsize/2);
-               std::vector<cln::cl_N>::iterator it = buf.begin();
+               auto it = buf.begin();
                for (int i=1; i<=xninitsize/2; i++) {
                        *it = bernoulli(i*2).to_cl_N();
                        it++;
@@ -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
@@ -320,7 +319,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
@@ -338,16 +337,20 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
                        // the switching point was empirically determined. the optimal point
                        // depends on hardware, Digits, ... so an approx value is okay.
                        // it solves also the problem with precision due to the u=-log(1-x) transformation
-                       if (cln::abs(cln::realpart(x)) < 0.25) {
-                               
+                       if (cln::abs(x) < 0.25) {
                                return Li2_do_sum(x);
                        } else {
+                               // Li2_do_sum practically doesn't converge near x == ±I
                                return Li2_do_sum_Xn(x);
                        }
                } else {
                        // choose the faster algorithm
                        if (cln::abs(cln::realpart(x)) > 0.75) {
-                               return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
+                               if ( x == 1 ) {
+                                       return cln::zeta(2);
+                               } else {
+                                       return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
+                               }
                        } else {
                                return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
                        }
@@ -363,15 +366,17 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
                if (cln::realpart(x) < 0.5) {
                        // choose the faster algorithm
                        // with n>=12 the "normal" summation always wins against the method with Xn
-                       if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
+                       if ((cln::abs(x) < 0.3) || (n >= 12)) {
                                return Lin_do_sum(n, x);
                        } else {
+                               // Li2_do_sum practically doesn't converge near x == ±I
                                return Lin_do_sum_Xn(n, x);
                        }
                } else {
-                       cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
+                       cln::cl_N result = 0;
+                       if ( x != 1 ) 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 +384,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 +402,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 +414,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 +444,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 +474,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 (const auto & it : x) {
+               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,51 +489,33 @@ 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];
 }
 
 
-// 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);
 
 
-// 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 +532,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)];
@@ -549,9 +540,9 @@ ex G_eval(const Gparameter& a, int scale)
        bool all_zero = true;
        bool all_ones = true;
        int count_ones = 0;
-       for (Gparameter::const_iterator it = a.begin(); it != a.end(); ++it) {
-               if (*it != 0) {
-                       const ex sym = gsyms[std::abs(*it)];
+       for (const auto & it : a) {
+               if (it != 0) {
+                       const ex sym = gsyms[std::abs(it)];
                        newa.append(sym);
                        all_zero = false;
                        if (sym != sc) {
@@ -569,38 +560,25 @@ ex G_eval(const Gparameter& a, int scale)
        // later on in the transformation
        if (newa.nops() > 1 && newa.op(0) == sc && !all_ones && a.front()!=0) {
                // do shuffle
-               Gparameter short_a;
-               Gparameter::const_iterator it = a.begin();
-               ++it;
-               for (; it != a.end(); ++it) {
-                       short_a.push_back(*it);
-               }
-               ex result = G_eval1(a.front(), scale) * G_eval(short_a, scale);
-               it = short_a.begin();
-               for (int i=1; i<count_ones; ++i) {
-                       ++it;
-               }
+               Gparameter short_a(a.begin()+1, a.end());
+               ex result = G_eval1(a.front(), scale, gsyms) * G_eval(short_a, scale, gsyms);
+
+               auto it = short_a.begin();
+               advance(it, count_ones-1);
                for (; it != short_a.end(); ++it) {
 
-                       Gparameter newa;
-                       Gparameter::const_iterator it2 = short_a.begin();
-                       for (--it2; it2 != it;) {
-                               ++it2;
-                               newa.push_back(*it2);
-                       }
+                       Gparameter newa(short_a.begin(), it);
+                       newa.push_back(*it);
                        newa.push_back(a[0]);
-                       ++it2;
-                       for (; it2 != short_a.end(); ++it2) {
-                               newa.push_back(*it2);   
-                       }
-                       result -= G_eval(newa, scale);
+                       newa.insert(newa.end(), it+1, short_a.end());
+                       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!
@@ -613,9 +591,9 @@ ex G_eval(const Gparameter& a, int scale)
        lst x;
        ex argbuf = gsyms[std::abs(scale)];
        ex mval = _ex1;
-       for (Gparameter::const_iterator it=a.begin(); it!=a.end(); ++it) {
-               if (*it != 0) {
-                       const ex& sym = gsyms[std::abs(*it)];
+       for (const auto & it : a) {
+               if (it != 0) {
+                       const ex& sym = gsyms[std::abs(it)];
                        x.append(argbuf / sym);
                        m.append(mval);
                        mval = _ex1;
@@ -651,14 +629,14 @@ Gparameter convert_pending_integrals_G(const Gparameter& pending_integrals)
 // trailing_zeros : number of trailing zeros of a
 // min_it         : iterator of a pointing on the smallest element in a
 Gparameter::const_iterator check_parameter_G(const Gparameter& a, int scale,
-               bool& convergent, int& depth, int& trailing_zeros, Gparameter::const_iterator& min_it)
+                                             bool& convergent, int& depth, int& trailing_zeros, Gparameter::const_iterator& min_it)
 {
        convergent = true;
        depth = 0;
        trailing_zeros = 0;
        min_it = a.end();
-       Gparameter::const_iterator lastnonzero = a.end();
-       for (Gparameter::const_iterator it = a.begin(); it != a.end(); ++it) {
+       auto lastnonzero = a.end();
+       for (auto it = a.begin(); it != a.end(); ++it) {
                if (std::abs(*it) > 0) {
                        ++depth;
                        trailing_zeros = 0;
@@ -673,6 +651,8 @@ Gparameter::const_iterator check_parameter_G(const Gparameter& a, int scale,
                        ++trailing_zeros;
                }
        }
+       if (lastnonzero == a.end())
+               return a.end();
        return ++lastnonzero;
 }
 
@@ -693,7 +673,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 +685,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);
-               for (Gparameter::const_iterator it = a.begin(); it != last; ++it) {
+               result += G_eval1(0, scale, gsyms) * trailing_zeros_G(new_a, scale, gsyms);
+               for (auto 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 +732,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 +758,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 +775,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 +789,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, bool flag_trailing_zeros_only);
 
 
 // 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, bool flag_trailing_zeros_only)
 {
        // main recursion routine
        //
@@ -818,21 +810,22 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
        bool convergent;
        int depth, trailing_zeros;
        Gparameter::const_iterator min_it;
-       Gparameter::const_iterator firstzero = 
-               check_parameter_G(a, scale, convergent, depth, trailing_zeros, min_it);
-       int min_it_pos = min_it - a.begin();
+       auto firstzero = check_parameter_G(a, scale, convergent, depth, trailing_zeros, min_it);
+       int min_it_pos = distance(a.begin(), min_it);
 
        // special case: all a's are zero
        if (depth == 0) {
                ex result;
 
                if (a.size() == 0) {
-                 result = 1;
+                       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,28 +834,30 @@ 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);
-               for (Gparameter::const_iterator it = a.begin(); it != firstzero; ++it) {
+               result += G_eval1(0, scale, gsyms) * G_transform(pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
+               for (auto 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, flag_trailing_zeros_only);
                }
                return result / trailing_zeros;
        }
 
-       // convergence case
-       if (convergent) {
+       // convergence case or flag_trailing_zeros_only
+       if (convergent || flag_trailing_zeros_only) {
                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 +872,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, flag_trailing_zeros_only)*
+                           G_transform(empty, a1, scale, gsyms, flag_trailing_zeros_only);
 
-               result -= shuffle_G(empty,a1,a2,pendint,a,scale);
+               result -= shuffle_G(empty, a1, a2, pendint, a, scale, gsyms, flag_trailing_zeros_only);
                return result;
        }
 
@@ -890,9 +886,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, flag_trailing_zeros_only);
        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 +898,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, flag_trailing_zeros_only);
                int buffer = *changeit;
                *changeit = *min_it;
-               result += G_transform(new_pendint, new_a, scale);
+               result += G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
                *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, flag_trailing_zeros_only);
                *changeit = *min_it;
-               result -= G_transform(new_pendint, new_a, scale);
+               result -= G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
        } 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, flag_trailing_zeros_only);
                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, flag_trailing_zeros_only);
                *changeit = *min_it;
-               result += G_transform(new_pendint, new_a, scale);
+               result += G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
        }
        return result;
 }
@@ -932,27 +933,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, bool flag_trailing_zeros_only)
 {
        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, flag_trailing_zeros_only);
        }
 
        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, flag_trailing_zeros_only);
        }
 
        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, flag_trailing_zeros_only);
        }
 
        Gparameter a1_removed(a1.begin()+1,a1.end());
@@ -964,203 +966,258 @@ 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, flag_trailing_zeros_only)
+            + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms, flag_trailing_zeros_only);
 }
 
-
 // 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+2)) {
-                               need_trafo = true;
-                               break;
+       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) {
+                               qlsts.push_back(1);
+                       } else {
+                               qlsts.push_back(-s[j-1]);
                        }
                }
+               if (qlstx.size() > 0) {
+                       buffer = buffer*G_numeric(qlstx, qlsts, 1/q);
+               }
+               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]);
+               }
+               if (plstx.size() > 0) {
+                       buffer = buffer*G_numeric(plstx, plsts, 1/p);
+               }
+               result = result + buffer;
        }
-       if (x.op(x.nops()-1).is_zero()) {
-               need_trafo = true;
-       }
-       if (depth == 1 && !need_trafo) {
-               return -Li(x.nops(), y / x.op(x.nops()-1)).evalf();
+       return result;
+}
+
+class less_object_for_cl_N
+{
+public:
+       bool operator() (const cln::cl_N & a, const cln::cl_N & b) const
+       {
+               // absolute value?
+               if (abs(a) != abs(b))
+                       return (abs(a) < abs(b)) ? true : false;
+
+               // complex phase?
+               if (phase(a) != phase(b))
+                       return (phase(a) < phase(b)) ? true : false;
+
+               // equal, therefore "less" is not true
+               return false;
        }
-       
-       // 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;
-               gsyms.clear();
-               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;
-                                       }
+};
+
+
+// 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, bool flag_trailing_zeros_only)
+{
+       // sort (|x|<->position) to determine indices
+       typedef std::multimap<cln::cl_N, std::size_t, less_object_for_cl_N> 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(x[i], i));
+                       ++size;
+               }
+       }
+       // include upper limit (scale)
+       sortmap.insert(std::make_pair(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);
                }
-
-               // do transformation
-               Gparameter pendint;
-               ex result = G_transform(pendint, a, scale);
-               // replace dummy symbols with their values
-               result = result.eval().expand();
-               result = result.subs(subslst).evalf();
-               
-               return result;
+               ++pos;
        }
 
-       // 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;
+       // do transformation
+       Gparameter pendint;
+       ex result = G_transform(pendint, a, scale, gsyms, flag_trailing_zeros_only);
+       // replace dummy symbols with their values
+       result = result.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;
+       bool have_trailing_zero = false;
+       std::size_t depth = 0;
+       for (auto & xi : x) {
+               if (!zerop(xi)) {
+                       ++depth;
+                       const cln::cl_N x_y = abs(xi) - 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(xi/y) - 1) < 0.01)
+                               need_hoelder = true;
                }
-               return result;
        }
+       if (zerop(x.back())) {
+               have_trailing_zero = true;
+               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 && !have_trailing_zero)
+               return G_do_hoelder(x, s, y);
+       
+       // convergence transformation
+       if (need_trafo)
+               return G_do_trafo(x, s, y, have_trailing_zero);
+
        // 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 (auto & xi : x) {
+               if (zerop(xi)) {
                        ++mcount;
                } else {
-                       newx.append(factor / (*it));
-                       factor = *it;
-                       m.append(mcount);
+                       newx.push_back(factor/xi);
+                       factor = xi;
+                       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;
-       for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
+       std::vector<cln::cl_N> newx;
+       newx.reserve(x.nops());
+       std::vector<int> s;
+       s.reserve(x.nops());
+       cln::cl_N factor(1);
+       for (auto 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);
+               }
+               const cln::cl_N xi = ex_to<numeric>(*itx).to_cl_N();
+               factor = factor/xi;
+               newx.push_back(factor);
+               if ( !instanceof(factor, cln::cl_R_ring) && imagpart(factor) < 0 ) {
+                       s.push_back(-1);
+               }
+               else {
+                       s.push_back(1);
                }
-               newx.append(factor / *itx);
-               factor /= *itx;
-               s.append(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)));
 }
 
 
@@ -1178,31 +1235,42 @@ ex mLi_numeric(const lst& m, const lst& x)
 
 static ex G2_evalf(const ex& x_, const ex& y)
 {
-       if (!y.info(info_flags::positive)) {
+       if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
                return G(x_, y).hold();
        }
-       lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
+       lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
        if (x.nops() == 0) {
                return _ex1;
        }
        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)) {
+       for (const auto & it : x) {
+               if (!it.info(info_flags::numeric)) {
                        return G(x_, y).hold();
                }
-               if (*it != _ex0) {
+               if (it != _ex0) {
                        all_zero = false;
                }
-               s.append(+1);
+               if ( !ex_to<numeric>(it).is_real() && ex_to<numeric>(it).imag() < 0 ) {
+                       s.push_back(-1);
+               }
+               else {
+                       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 (const auto & it : x)
+               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);
 }
 
 
@@ -1210,30 +1278,36 @@ static ex G2_eval(const ex& x_, const ex& y)
 {
        //TODO eval to MZV or H or S or Lin
 
-       if (!y.info(info_flags::positive)) {
+       if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
                return G(x_, y).hold();
        }
-       lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
+       lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
        if (x.nops() == 0) {
                return _ex1;
        }
        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) {
-               if (!(*it).info(info_flags::numeric)) {
+       for (const auto & it : x) {
+               if (!it.info(info_flags::numeric)) {
                        return G(x_, y).hold();
                }
-               if (!(*it).info(info_flags::crational)) {
+               if (!it.info(info_flags::crational)) {
                        crational = false;
                }
-               if (*it != _ex0) {
+               if (it != _ex0) {
                        all_zero = false;
                }
-               s.append(+1);
+               if ( !ex_to<numeric>(it).is_real() && ex_to<numeric>(it).imag() < 0 ) {
+                       s.push_back(-1);
+               }
+               else {
+                       s.push_back(+1);
+               }
        }
        if (all_zero) {
                return pow(log(y), x.nops()) / factorial(x.nops());
@@ -1244,14 +1318,19 @@ 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 (const auto & it : x)
+               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);
 }
 
 
+// option do_not_evalf_params() removed.
 unsigned G2_SERIAL::serial = function::register_new(function_options("G", 2).
                                 evalf_func(G2_evalf).
                                 eval_func(G2_eval).
-                                do_not_evalf_params().
                                 overloaded(2));
 //TODO
 //                                derivative_func(G2_deriv).
@@ -1260,11 +1339,11 @@ unsigned G2_SERIAL::serial = function::register_new(function_options("G", 2).
 
 static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
 {
-       if (!y.info(info_flags::positive)) {
+       if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
                return G(x_, s_, y).hold();
        }
-       lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
-       lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst(s_);
+       lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
+       lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst{s_};
        if (x.nops() != s.nops()) {
                return G(x_, s_, y).hold();
        }
@@ -1274,9 +1353,10 @@ 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) {
+       for (auto itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
                if (!(*itx).info(info_flags::numeric)) {
                        return G(x_, y).hold();
                }
@@ -1286,16 +1366,36 @@ 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 ( ex_to<numeric>(*itx).is_positive() ) {
+                               if ( *its >= 0 ) {
+                                       sn.push_back(1);
+                               }
+                               else {
+                                       sn.push_back(-1);
+                               }
+                       } else {
+                               sn.push_back(1);
+                       }
+               }
+               else {
+                       if ( ex_to<numeric>(*itx).imag() > 0 ) {
+                               sn.push_back(1);
+                       }
+                       else {
+                               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 (const auto & it : x)
+               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);
 }
 
 
@@ -1303,11 +1403,11 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
 {
        //TODO eval to MZV or H or S or Lin
 
-       if (!y.info(info_flags::positive)) {
+       if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
                return G(x_, s_, y).hold();
        }
-       lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
-       lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst(s_);
+       lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
+       lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst{s_};
        if (x.nops() != s.nops()) {
                return G(x_, s_, y).hold();
        }
@@ -1317,10 +1417,11 @@ 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) {
+       for (auto itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
                if (!(*itx).info(info_flags::numeric)) {
                        return G(x_, s_, y).hold();
                }
@@ -1333,10 +1434,25 @@ 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 ( ex_to<numeric>(*itx).is_positive() ) {
+                               if ( *its >= 0 ) {
+                                       sn.push_back(1);
+                               }
+                               else {
+                                       sn.push_back(-1);
+                               }
+                       } else {
+                               sn.push_back(1);
+                       }
+               }
+               else {
+                       if ( ex_to<numeric>(*itx).imag() > 0 ) {
+                               sn.push_back(1);
+                       }
+                       else {
+                               sn.push_back(-1);
+                       }
                }
        }
        if (all_zero) {
@@ -1348,14 +1464,21 @@ 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 (const auto & it : x)
+               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);
 }
 
 
+// option do_not_evalf_params() removed.
+// This is safe: in the code above it only matters if s_ > 0 or s_ < 0,
+// s_ is allowed to be of floating type.
 unsigned G3_SERIAL::serial = function::register_new(function_options("G", 3).
                                 evalf_func(G3_evalf).
                                 eval_func(G3_eval).
-                                do_not_evalf_params().
                                 overloaded(2));
 //TODO
 //                                derivative_func(G3_deriv).
@@ -1376,12 +1499,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);
                        }
                }
        }
@@ -1400,7 +1529,7 @@ static ex Li_evalf(const ex& m_, const ex& x_)
                        return Li(m_,x_).hold();
                }
 
-               for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
+               for (auto itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
                        if (!(*itm).info(info_flags::posint)) {
                                return Li(m_, x_).hold();
                        }
@@ -1436,7 +1565,7 @@ static ex Li_eval(const ex& m_, const ex& x_)
                        bool is_zeta = true;
                        bool do_evalf = true;
                        bool crational = true;
-                       for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
+                       for (auto itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
                                if (!(*itm).info(info_flags::posint)) {
                                        return Li(m_,x_).hold();
                                }
@@ -1457,7 +1586,17 @@ static ex Li_eval(const ex& m_, const ex& x_)
                                }
                        }
                        if (is_zeta) {
-                               return zeta(m_,x_);
+                               lst newx;
+                               for (const auto & itx : x) {
+                                       GINAC_ASSERT((itx == _ex1) || (itx == _ex_1));
+                                       // XXX: 1 + 0.0*I is considered equal to 1. However
+                                       // the former is a not automatically converted
+                                       // to a real number. Do the conversion explicitly
+                                       // to avoid the "numeric::operator>(): complex inequality"
+                                       // exception (and similar problems).
+                                       newx.append(itx != _ex_1 ? _ex1 : _ex_1);
+                               }
+                               return zeta(m_, newx);
                        }
                        if (is_H) {
                                ex prefactor;
@@ -1495,7 +1634,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 +1646,35 @@ 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 { expair(Li(m, x), 0) };
+               return pseries(rel, std::move(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 { expair(Order(_ex1), order) };
+                       ser += pseries(rel, std::move(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()
 }
 
 
@@ -1545,16 +1713,16 @@ static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c)
        if (is_a<lst>(m_)) {
                m = ex_to<lst>(m_);
        } else {
-               m = lst(m_);
+               m = lst{m_};
        }
        lst x;
        if (is_a<lst>(x_)) {
                x = ex_to<lst>(x_);
        } else {
-               x = lst(x_);
+               x = lst{x_};
        }
-       c.s << "\\mbox{Li}_{";
-       lst::const_iterator itm = m.begin();
+       c.s << "\\mathrm{Li}_{";
+       auto itm = m.begin();
        (*itm).print(c);
        itm++;
        for (; itm != m.end(); itm++) {
@@ -1562,7 +1730,7 @@ static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c)
                (*itm).print(c);
        }
        c.s << "}(";
-       lst::const_iterator itx = x.begin();
+       auto itx = x.begin();
        (*itx).print(c);
        itx++;
        for (; itx != x.end(); itx++) {
@@ -1597,7 +1765,7 @@ namespace {
 
 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
 // see fill_Yn()
-std::vector<std::vector<cln::cl_N> > Yn;
+std::vector<std::vector<cln::cl_N>> Yn;
 int ynsize = 0; // number of Yn[]
 int ynlength = 100; // initial length of all Yn[i]
 
@@ -1618,8 +1786,8 @@ void fill_Yn(int n, const cln::float_format_t& prec)
 
        if (n) {
                std::vector<cln::cl_N> buf(initsize);
-               std::vector<cln::cl_N>::iterator it = buf.begin();
-               std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
+               auto it = buf.begin();
+               auto itprev = Yn[n-1].begin();
                *it = (*itprev) / cln::cl_N(n+1) * one;
                it++;
                itprev++;
@@ -1633,7 +1801,7 @@ void fill_Yn(int n, const cln::float_format_t& prec)
                Yn.push_back(buf);
        } else {
                std::vector<cln::cl_N> buf(initsize);
-               std::vector<cln::cl_N>::iterator it = buf.begin();
+               auto it = buf.begin();
                *it = 1 * one;
                it++;
                for (int i=2; i<=initsize; i++) {
@@ -1653,7 +1821,7 @@ void make_Yn_longer(int newsize, const cln::float_format_t& prec)
        cln::cl_N one = cln::cl_float(1, prec);
 
        Yn[0].resize(newsize);
-       std::vector<cln::cl_N>::iterator it = Yn[0].begin();
+       auto it = Yn[0].begin();
        it += ynlength;
        for (int i=ynlength+1; i<=newsize; i++) {
                *it = *(it-1) + 1 / cln::cl_N(i) * one;
@@ -1662,8 +1830,8 @@ void make_Yn_longer(int newsize, const cln::float_format_t& prec)
 
        for (int n=1; n<ynsize; n++) {
                Yn[n].resize(newsize);
-               std::vector<cln::cl_N>::iterator it = Yn[n].begin();
-               std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
+               auto it = Yn[n].begin();
+               auto itprev = Yn[n-1].begin();
                it += ynlength;
                itprev += ynlength;
                for (int i=ynlength+n+1; i<=newsize+n; i++) {
@@ -1688,10 +1856,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 +1867,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 +1945,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 +2006,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 +2017,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 +2053,17 @@ 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)) {
+       // the condition abs(1-value)>1 avoids an infinite recursion in the region abs(value)<=1 && abs(value)>0.95 && abs(1-value)<=1 && abs(1-value)>0.95
+       // we don't care here about abs(value)<1 && real(value)>0.5, this will be taken care of in S_projection
+       if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95) && (cln::abs(1-value) > 1) )) {
 
                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 +2072,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 +2089,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);
@@ -1924,6 +2104,16 @@ numeric S_num(int n, int p, const numeric& x)
 
                return result;
        }
+
+       if ((cln::abs(value) > 0.95) && (cln::abs(value-9.53) < 9.47)) {
+               lst m;
+               m.append(n+1);
+               for (int s=0; s<p-1; s++)
+                       m.append(1);
+
+               ex res = H(m,numeric(value)).evalf();
+               return ex_to<numeric>(res).to_cl_N();
+       }
        else {
                return S_projection(n, p, value, prec);
        }
@@ -1945,12 +2135,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);
                        }
                }
        }
@@ -1965,7 +2161,7 @@ static ex S_eval(const ex& n, const ex& p, const ex& x)
                        return _ex0;
                }
                if (x == 1) {
-                       lst m(n+1);
+                       lst m{n+1};
                        for (int i=ex_to<numeric>(p).to_int()-1; i>0; i--) {
                                m.append(1);
                        }
@@ -1975,7 +2171,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 +2188,47 @@ 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 { expair(Order(_ex1), order) };
+                       ser += pseries(rel, std::move(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()
 }
 
 
@@ -2010,7 +2248,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);
@@ -2041,7 +2279,7 @@ REGISTER_FUNCTION(S,
 // anonymous namespace for helper functions
 namespace {
 
-       
+
 // regulates the pole (used by 1/x-transformation)
 symbol H_polesign("IMSIGN");
 
@@ -2053,19 +2291,19 @@ bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
 {
        // expand parameter list
        lst mexp;
-       for (lst::const_iterator it = l.begin(); it != l.end(); it++) {
-               if (*it > 1) {
-                       for (ex count=*it-1; count > 0; count--) {
+       for (const auto & it : l) {
+               if (it > 1) {
+                       for (ex count=it-1; count > 0; count--) {
                                mexp.append(0);
                        }
                        mexp.append(1);
-               } else if (*it < -1) {
-                       for (ex count=*it+1; count < 0; count++) {
+               } else if (it < -1) {
+                       for (ex count=it+1; count < 0; count++) {
                                mexp.append(0);
                        }
                        mexp.append(-1);
                } else {
-                       mexp.append(*it);
+                       mexp.append(it);
                }
        }
        
@@ -2073,25 +2311,25 @@ bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
        pf = 1;
        bool has_negative_parameters = false;
        ex acc = 1;
-       for (lst::const_iterator it = mexp.begin(); it != mexp.end(); it++) {
-               if (*it == 0) {
+       for (const auto & it : mexp) {
+               if (it == 0) {
                        acc++;
                        continue;
                }
-               if (*it > 0) {
-                       m.append((*it+acc-1) * signum);
+               if (it > 0) {
+                       m.append((it+acc-1) * signum);
                } else {
-                       m.append((*it-acc+1) * signum);
+                       m.append((it-acc+1) * signum);
                }
                acc = 1;
-               signum = *it;
-               pf *= *it;
+               signum = it;
+               pf *= it;
                if (pf < 0) {
                        has_negative_parameters = true;
                }
        }
        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);
@@ -2108,7 +2346,7 @@ bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
 // recursivly transforms H to corresponding multiple polylogarithms
 struct map_trafo_H_convert_to_Li : public map_function
 {
-       ex operator()(const ex& e)
+       ex operator()(const ex& e) override
        {
                if (is_a<add>(e) || is_a<mul>(e)) {
                        return e.map(*this);
@@ -2118,9 +2356,9 @@ struct map_trafo_H_convert_to_Li : public map_function
                        if (name == "H") {
                                lst parameter;
                                if (is_a<lst>(e.op(0))) {
-                                               parameter = ex_to<lst>(e.op(0));
+                                       parameter = ex_to<lst>(e.op(0));
                                } else {
-                                       parameter = lst(e.op(0));
+                                       parameter = lst{e.op(0)};
                                }
                                ex arg = e.op(1);
 
@@ -2131,7 +2369,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;
@@ -2147,7 +2385,7 @@ struct map_trafo_H_convert_to_Li : public map_function
 // recursivly transforms H to corresponding zetas
 struct map_trafo_H_convert_to_zeta : public map_function
 {
-       ex operator()(const ex& e)
+       ex operator()(const ex& e) override
        {
                if (is_a<add>(e) || is_a<mul>(e)) {
                        return e.map(*this);
@@ -2157,9 +2395,9 @@ struct map_trafo_H_convert_to_zeta : public map_function
                        if (name == "H") {
                                lst parameter;
                                if (is_a<lst>(e.op(0))) {
-                                               parameter = ex_to<lst>(e.op(0));
+                                       parameter = ex_to<lst>(e.op(0));
                                } else {
-                                       parameter = lst(e.op(0));
+                                       parameter = lst{e.op(0)};
                                }
 
                                lst m;
@@ -2180,7 +2418,7 @@ struct map_trafo_H_convert_to_zeta : public map_function
 // remove trailing zeros from H-parameters
 struct map_trafo_H_reduce_trailing_zeros : public map_function
 {
-       ex operator()(const ex& e)
+       ex operator()(const ex& e) override
        {
                if (is_a<add>(e) || is_a<mul>(e)) {
                        return e.map(*this);
@@ -2192,7 +2430,7 @@ struct map_trafo_H_reduce_trailing_zeros : public map_function
                                if (is_a<lst>(e.op(0))) {
                                        parameter = ex_to<lst>(e.op(0));
                                } else {
-                                       parameter = lst(e.op(0));
+                                       parameter = lst{e.op(0)};
                                }
                                ex arg = e.op(1);
                                if (parameter.op(parameter.nops()-1) == 0) {
@@ -2203,7 +2441,7 @@ struct map_trafo_H_reduce_trailing_zeros : public map_function
                                        }
                                        
                                        //
-                                       lst::const_iterator it = parameter.begin();
+                                       auto it = parameter.begin();
                                        while ((it != parameter.end()) && (*it == 0)) {
                                                it++;
                                        }
@@ -2213,7 +2451,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--;
                                        }
@@ -2265,14 +2503,19 @@ ex convert_H_to_zeta(const lst& m)
 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf)
 {
        lst res;
-       lst::const_iterator itm = m.begin();
-       lst::const_iterator itx = ++x.begin();
+       auto itm = m.begin();
+       auto itx = ++x.begin();
        int signum = 1;
        pf = _ex1;
        res.append(*itm);
        itm++;
        while (itx != x.end()) {
-               signum *= (*itx > 0) ? 1 : -1;
+               GINAC_ASSERT((*itx == _ex1) || (*itx == _ex_1));
+               // XXX: 1 + 0.0*I is considered equal to 1. However the former
+               // is not automatically converted to a real number.
+               // Do the conversion explicitly to avoid the
+               // "numeric::operator>(): complex inequality" exception.
+               signum *= (*itx != _ex_1) ? 1 : -1;
                pf *= signum;
                res.append((*itm) * signum);
                itm++;
@@ -2299,12 +2542,12 @@ ex trafo_H_mult(const ex& h1, const ex& h2)
                if (h2nops > 1) {
                        hlong = ex_to<lst>(h2.op(0));
                } else {
-                       hlong = h2.op(0).op(0);
+                       hlong = lst{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]);
                }
@@ -2321,7 +2564,7 @@ ex trafo_H_mult(const ex& h1, const ex& h2)
 // applies trafo_H_mult recursively on expressions
 struct map_trafo_H_mult : public map_function
 {
-       ex operator()(const ex& e)
+       ex operator()(const ex& e) override
        {
                if (is_a<add>(e)) {
                        return e.map(*this);
@@ -2332,7 +2575,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") {
@@ -2366,7 +2609,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();
@@ -2394,7 +2637,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") {
@@ -2409,7 +2652,38 @@ ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
                ex addzeta = convert_H_to_zeta(newparameter);
                return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
        } else {
-               return e * (-H(lst(0),1/arg).hold());
+               return e * (-H(lst{ex(0)},1/arg).hold());
+       }
+}
+
+
+// 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 (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") {
+                                       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{ex(1)},1-arg).hold();
        }
 }
 
@@ -2426,7 +2700,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") {
@@ -2441,8 +2715,8 @@ ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
                ex addzeta = convert_H_to_zeta(newparameter);
                return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
        } else {
-               ex addzeta = convert_H_to_zeta(lst(-1));
-               return (e * (addzeta - H(lst(-1),1/arg).hold())).expand();
+               ex addzeta = convert_H_to_zeta(lst{ex(-1)});
+               return (e * (addzeta - H(lst{ex(-1)},1/arg).hold())).expand();
        }
 }
 
@@ -2459,7 +2733,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") {
@@ -2473,7 +2747,7 @@ ex trafo_H_1mxt1px_prepend_minusone(const ex& e, const ex& arg)
                newparameter.prepend(-1);
                return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
        } else {
-               return (e * H(lst(-1),(1-arg)/(1+arg)).hold()).expand();
+               return (e * H(lst{ex(-1)},(1-arg)/(1+arg)).hold()).expand();
        }
 }
 
@@ -2490,7 +2764,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") {
@@ -2504,15 +2778,116 @@ ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
                newparameter.prepend(1);
                return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
        } else {
-               return (e * H(lst(1),(1-arg)/(1+arg)).hold()).expand();
+               return (e * H(lst{ex(1)},(1-arg)/(1+arg)).hold()).expand();
        }
 }
 
 
+// do x -> 1-x transformation
+struct map_trafo_H_1mx : public map_function
+{
+       ex operator()(const ex& e) override
+       {
+               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 (std::size_t 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 (std::size_t 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 (std::size_t 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{ex(1)}, arg).hold() * H(newparameter, arg).hold();
+                                       std::size_t firstzero = 0;
+                                       while (parameter.op(firstzero) == 1) {
+                                               firstzero++;
+                                       }
+                                       for (std::size_t i = firstzero-1; i < parameter.nops()-1; i++) {
+                                               lst newparameter;
+                                               std::size_t 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
 {
-       ex operator()(const ex& e)
+       ex operator()(const ex& e) override
        {
                if (is_a<add>(e) || is_a<mul>(e)) {
                        return e.map(*this);
@@ -2528,7 +2903,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;
@@ -2538,7 +2913,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;
@@ -2546,11 +2921,11 @@ struct map_trafo_H_1overx : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(H(lst(-1),1/arg).hold() - H(lst(0),1/arg).hold(), parameter.nops())
+                                               return unify((pow(H(lst{ex(-1)},1/arg).hold() - H(lst{ex(0)},1/arg).hold(), parameter.nops())
                                                       / 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;
@@ -2558,7 +2933,7 @@ struct map_trafo_H_1overx : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() + H_polesign, parameter.nops())
+                                               return unify((pow(H(lst{ex(1)},1/arg).hold() + H(lst{ex(0)},1/arg).hold() + H_polesign, parameter.nops())
                                                       / factorial(parameter.nops())).expand());
                                        }
                                }
@@ -2573,7 +2948,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 {
@@ -2588,7 +2963,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 {
@@ -2601,14 +2976,14 @@ struct map_trafo_H_1overx : public map_function
                                        // leading one
                                        map_trafo_H_1overx recursion;
                                        map_trafo_H_mult unify;
-                                       ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
-                                       int firstzero = 0;
+                                       ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
+                                       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]);
                                                }
@@ -2633,7 +3008,7 @@ struct map_trafo_H_1overx : public map_function
 // do x -> (1-x)/(1+x) transformation
 struct map_trafo_H_1mxt1px : public map_function
 {
-       ex operator()(const ex& e)
+       ex operator()(const ex& e) override
        {
                if (is_a<add>(e) || is_a<mul>(e)) {
                        return e.map(*this);
@@ -2649,7 +3024,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;
@@ -2657,11 +3032,11 @@ struct map_trafo_H_1mxt1px : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(-H(lst(1),(1-arg)/(1+arg)).hold() - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
+                                               return unify((pow(-H(lst{ex(1)},(1-arg)/(1+arg)).hold() - H(lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.nops())
                                                       / 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;
@@ -2669,11 +3044,11 @@ struct map_trafo_H_1mxt1px : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(log(2) - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
+                                               return unify((pow(log(2) - H(lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.nops())
                                                       / 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;
@@ -2681,7 +3056,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(-log(2) - H(lst(0),(1-arg)/(1+arg)).hold() + H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
+                                               return unify((pow(-log(2) - H(lst{ex(0)},(1-arg)/(1+arg)).hold() + H(lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.nops())
                                                       / factorial(parameter.nops())).expand());
                                        }
                                }
@@ -2696,7 +3071,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 {
@@ -2711,7 +3086,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 {
@@ -2724,14 +3099,14 @@ struct map_trafo_H_1mxt1px : public map_function
                                        // leading one
                                        map_trafo_H_1mxt1px recursion;
                                        map_trafo_H_mult unify;
-                                       ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
-                                       int firstzero = 0;
+                                       ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
+                                       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]);
                                                }
@@ -2805,7 +3180,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();
                        }
@@ -2822,33 +3197,27 @@ 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) {
-                               for (ex count=*it-1; count > 0; count--) {
+               for (const auto & it : morg) {
+                       if (it > 1) {
+                               for (ex count=it-1; count > 0; count--) {
                                        m.append(0);
                                }
                                m.append(1);
-                       } else if (*it < -1) {
-                               for (ex count=*it+1; count < 0; count++) {
+                       } 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);
+                               m.append(it);
                        }
                }
 
-               // check for the applicability of transformations
-               // 
-               // first condition: since the transformations produce a lot of terms,
-               //   they are only efficient for argument near the boundary |x| = 1, then no
-               //   transformation is needed
-               // second condition: veto for region around +-I to avoid endless recursion
-               //   with the (1-x)/(1+x) transformation. 1.198 is sqrt(1.4142) is the
-               //   boundary of the problematic transformation.
-               //   
-               if (cln::abs(x) < 0.95 || (cln::abs(x) < 1 && cln::abs(x-1) >= 1.198)) {
+               // do summation
+               if (cln::abs(x) < 0.95) {
                        lst m_lst;
                        lst s_lst;
                        ex pf;
@@ -2856,7 +3225,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                                // negative parameters -> s_lst is filled
                                std::vector<int> m_int;
                                std::vector<cln::cl_N> x_cln;
-                               for (lst::const_iterator it_int = m_lst.begin(), it_cln = s_lst.begin(); 
+                               for (auto it_int = m_lst.begin(), it_cln = s_lst.begin();
                                     it_int != m_lst.end(); it_int++, it_cln++) {
                                        m_int.push_back(ex_to<numeric>(*it_int).to_int());
                                        x_cln.push_back(ex_to<numeric>(*it_cln).to_cl_N());
@@ -2870,19 +3239,20 @@ static ex H_evalf(const ex& x1, const ex& x2)
                                        return Li(m_lst.op(0), x2).evalf();
                                }
                                std::vector<int> m_int;
-                               for (lst::const_iterator it = m_lst.begin(); it != m_lst.end(); it++) {
-                                       m_int.push_back(ex_to<numeric>(*it).to_int());
+                               for (const auto & it : m_lst) {
+                                       m_int.push_back(ex_to<numeric>(it).to_int());
                                }
                                return numeric(H_do_sum(m_int, x));
                        }
                }
 
+               symbol xtemp("xtemp");
                ex res = 1;     
                
                // 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;
@@ -2890,32 +3260,34 @@ 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));
+                       res *= trafo(H(m, xtemp).hold());
                        if (cln::imagpart(x) <= 0) {
                                res = res.subs(H_polesign == -I*Pi);
                        } 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).hold());
+               } 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).hold());
                }
-
-               // 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();
        }
@@ -2930,7 +3302,7 @@ static ex H_eval(const ex& m_, const ex& x)
        if (is_a<lst>(m_)) {
                m = ex_to<lst>(m_);
        } else {
-               m = lst(m_);
+               m = lst{m_};
        }
        if (m.nops() == 0) {
                return _ex1;
@@ -2959,8 +3331,8 @@ static ex H_eval(const ex& m_, const ex& x)
                pos1 = *m.begin();
                p = _ex1;
        }
-       for (lst::const_iterator it = ++m.begin(); it != m.end(); it++) {
-               if ((*it).info(info_flags::integer)) {
+       for (auto it = ++m.begin(); it != m.end(); it++) {
+               if (it->info(info_flags::integer)) {
                        if (step == 0) {
                                if (*it > _ex1) {
                                        if (pos1 == _ex0) {
@@ -3041,9 +3413,8 @@ static ex H_eval(const ex& m_, const ex& x)
 
 static ex H_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
 {
-       epvector seq;
-       seq.push_back(expair(H(m, x), 0));
-       return pseries(rel, seq);
+       epvector seq { expair(H(m, x), 0) };
+       return pseries(rel, std::move(seq));
 }
 
 
@@ -3057,7 +3428,7 @@ static ex H_deriv(const ex& m_, const ex& x, unsigned deriv_param)
        if (is_a<lst>(m_)) {
                m = ex_to<lst>(m_);
        } else {
-               m = lst(m_);
+               m = lst{m_};
        }
        ex mb = *m.begin();
        if (mb > _ex1) {
@@ -3085,10 +3456,10 @@ static void H_print_latex(const ex& m_, const ex& x, const print_context& c)
        if (is_a<lst>(m_)) {
                m = ex_to<lst>(m_);
        } else {
-               m = lst(m_);
+               m = lst{m_};
        }
-       c.s << "\\mbox{H}_{";
-       lst::const_iterator itm = m.begin();
+       c.s << "\\mathrm{H}_{";
+       auto itm = m.begin();
        (*itm).print(c);
        itm++;
        for (; itm != m.end(); itm++) {
@@ -3118,7 +3489,7 @@ ex convert_H_to_Li(const ex& m, const ex& x)
        if (is_a<lst>(m)) {
                return filter2(filter(H(m, x).hold()));
        } else {
-               return filter2(filter(H(lst(m), x).hold()));
+               return filter2(filter(H(lst{m}, x).hold()));
        }
 }
 
@@ -3138,13 +3509,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)
 {
@@ -3159,44 +3523,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);
@@ -3214,14 +3573,12 @@ 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();
+       auto it = f_kj.begin();
        cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
        
        t0 = cln::exp(-lambda);
@@ -3244,7 +3601,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();
 
@@ -3285,6 +3643,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;
@@ -3292,6 +3652,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;
@@ -3318,7 +3679,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);
 
@@ -3332,12 +3694,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()--;
                        
@@ -3351,9 +3714,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;
 }
@@ -3395,7 +3760,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];
@@ -3490,8 +3855,8 @@ static ex zeta1_evalf(const ex& x)
                std::vector<int> r(count);
 
                // check parameters and convert them
-               lst::const_iterator it1 = xlst.begin();
-               std::vector<int>::iterator it2 = r.begin();
+               auto it1 = xlst.begin();
+               auto it2 = r.begin();
                do {
                        if (!(*it1).info(info_flags::posint)) {
                                return zeta(x).hold();
@@ -3586,7 +3951,7 @@ static void zeta1_print_latex(const ex& m_, const print_context& c)
        c.s << "\\zeta(";
        if (is_a<lst>(m_)) {
                const lst& m = ex_to<lst>(m_);
-               lst::const_iterator it = m.begin();
+               auto it = m.begin();
                (*it).print(c);
                it++;
                for (; it != m.end(); it++) {
@@ -3630,10 +3995,10 @@ static ex zeta2_evalf(const ex& x, const ex& s)
                std::vector<int> si(count);
 
                // check parameters and convert them
-               lst::const_iterator it_xread = xlst.begin();
-               lst::const_iterator it_sread = slst.begin();
-               std::vector<int>::iterator it_xwrite = xi.begin();
-               std::vector<int>::iterator it_swrite = si.begin();
+               auto it_xread = xlst.begin();
+               auto it_sread = slst.begin();
+               auto it_xwrite = xi.begin();
+               auto it_swrite = si.begin();
                do {
                        if (!(*it_xread).info(info_flags::posint)) {
                                return zeta(x, s).hold();
@@ -3667,8 +4032,8 @@ static ex zeta2_eval(const ex& m, const ex& s_)
 {
        if (is_exactly_a<lst>(s_)) {
                const lst& s = ex_to<lst>(s_);
-               for (lst::const_iterator it = s.begin(); it != s.end(); it++) {
-                       if ((*it).info(info_flags::positive)) {
+               for (const auto & it : s) {
+                       if (it.info(info_flags::positive)) {
                                continue;
                        }
                        return zeta(m, s_).hold();
@@ -3703,17 +4068,17 @@ static void zeta2_print_latex(const ex& m_, const ex& s_, const print_context& c
        if (is_a<lst>(m_)) {
                m = ex_to<lst>(m_);
        } else {
-               m = lst(m_);
+               m = lst{m_};
        }
        lst s;
        if (is_a<lst>(s_)) {
                s = ex_to<lst>(s_);
        } else {
-               s = lst(s_);
+               s = lst{s_};
        }
        c.s << "\\zeta(";
-       lst::const_iterator itm = m.begin();
-       lst::const_iterator its = s.begin();
+       auto itm = m.begin();
+       auto its = s.begin();
        if (*its < 0) {
                c.s << "\\overline{";
                (*itm).print(c);