]> www.ginac.de Git - ginac.git/blobdiff - ginac/inifcns_nstdsums.cpp
Fix G3_evalf and G3_eval for real and negative x.
[ginac.git] / ginac / inifcns_nstdsums.cpp
index d85044d04bd1f311b2d26ad5ab8ccf0572d1ebe7..c33a4cf8370d22020e302b0c08e1ade04f6a24f1 100644 (file)
@@ -47,7 +47,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2011 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 {
 
@@ -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
@@ -347,7 +346,11 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
                } 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);
                        }
@@ -369,7 +372,8 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
                                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) - S_num(1, n-j-1, 1-x))
                                                  * cln::expt(cln::log(x), j) / cln::factorial(j);
@@ -501,19 +505,6 @@ cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl
 }
 
 
-// converts parameter types and calls multipleLi_do_sum (convenience function for G_numeric)
-cln::cl_N mLi_do_summation(const lst& m, const lst& x)
-{
-       std::vector<int> m_int;
-       std::vector<cln::cl_N> x_cln;
-       for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
-               m_int.push_back(ex_to<numeric>(*itm).to_int());
-               x_cln.push_back(ex_to<numeric>(*itx).to_cl_N());
-       }
-       return multipleLi_do_sum(m_int, x_cln);
-}
-
-
 // forward declaration for Li_eval()
 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf);
 
@@ -583,11 +574,12 @@ ex G_eval(const Gparameter& a, int scale, const exvector& gsyms)
 
                        Gparameter newa;
                        Gparameter::const_iterator it2 = short_a.begin();
-                       for (--it2; it2 != it;) {
-                               ++it2;
+                       for (; it2 != it; ++it2) {
                                newa.push_back(*it2);
                        }
+                       newa.push_back(*it);
                        newa.push_back(a[0]);
+                       it2 = it;
                        ++it2;
                        for (; it2 != short_a.end(); ++it2) {
                                newa.push_back(*it2);   
@@ -672,6 +664,8 @@ Gparameter::const_iterator check_parameter_G(const Gparameter& a, int scale,
                        ++trailing_zeros;
                }
        }
+       if (lastnonzero == a.end())
+               return a.end();
        return ++lastnonzero;
 }
 
@@ -990,199 +984,233 @@ ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2
             + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms);
 }
 
-
 // handles the transformations and the numerical evaluation of G
 // the parameter x, s and y must only contain numerics
-ex G_numeric(const lst& x, const lst& s, const ex& y)
+static cln::cl_N
+G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
+         const cln::cl_N& y);
+
+// do acceleration transformation (hoelder convolution [BBB])
+// the parameter x, s and y must only contain numerics
+static cln::cl_N
+G_do_hoelder(std::vector<cln::cl_N> x, /* yes, it's passed by value */
+            const std::vector<int>& s, const cln::cl_N& y)
 {
-       // check for convergence and necessary accelerations
-       bool need_trafo = false;
-       bool need_hoelder = false;
-       int depth = 0;
-       for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
-               if (!(*it).is_zero()) {
-                       ++depth;
-                       if (abs(*it) - y < -pow(10,-Digits+1)) {
-                               need_trafo = true;
+       cln::cl_N result;
+       const std::size_t size = x.size();
+       for (std::size_t i = 0; i < size; ++i)
+               x[i] = x[i]/y;
+
+       for (std::size_t r = 0; r <= size; ++r) {
+               cln::cl_N buffer(1 & r ? -1 : 1);
+               cln::cl_RA p(2);
+               bool adjustp;
+               do {
+                       adjustp = false;
+                       for (std::size_t i = 0; i < size; ++i) {
+                               if (x[i] == cln::cl_RA(1)/p) {
+                                       p = p/2 + cln::cl_RA(3)/2;
+                                       adjustp = true;
+                                       continue;
+                               }
                        }
-                       if (abs((abs(*it) - y)/y) < 0.01) {
-                               need_hoelder = true;
+               } while (adjustp);
+               cln::cl_RA q = p/(p-1);
+               std::vector<cln::cl_N> qlstx;
+               std::vector<int> qlsts;
+               for (std::size_t j = r; j >= 1; --j) {
+                       qlstx.push_back(cln::cl_N(1) - x[j-1]);
+                       if (instanceof(x[j-1], cln::cl_R_ring) &&
+                           realpart(x[j-1]) > 1 && realpart(x[j-1]) <= 2) {
+                               qlsts.push_back(s[j-1]);
+                       } else {
+                               qlsts.push_back(-s[j-1]);
                        }
                }
-       }
-       if (x.op(x.nops()-1).is_zero()) {
-               need_trafo = true;
-       }
-       if (depth == 1 && x.nops() == 2 && !need_trafo) {
-               return -Li(x.nops(), y / x.op(x.nops()-1)).evalf();
-       }
-       
-       // do acceleration transformation (hoelder convolution [BBB])
-       if (need_hoelder) {
-               
-               ex result;
-               const int size = x.nops();
-               lst newx;
-               for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
-                       newx.append(*it / y);
+               if (qlstx.size() > 0) {
+                       buffer = buffer*G_numeric(qlstx, qlsts, 1/q);
                }
-               
-               for (int r=0; r<=size; ++r) {
-                       ex buffer = pow(-1, r);
-                       ex p = 2;
-                       bool adjustp;
-                       do {
-                               adjustp = false;
-                               for (lst::const_iterator it = newx.begin(); it != newx.end(); ++it) {
-                                       if (*it == 1/p) {
-                                               p += (3-p)/2; 
-                                               adjustp = true;
-                                               continue;
-                                       }
-                               }
-                       } while (adjustp);
-                       ex q = p / (p-1);
-                       lst qlstx;
-                       lst qlsts;
-                       for (int j=r; j>=1; --j) {
-                               qlstx.append(1-newx.op(j-1));
-                               if (newx.op(j-1).info(info_flags::real) && newx.op(j-1) > 1 && newx.op(j-1) <= 2) {
-                                       qlsts.append( s.op(j-1));
-                               } else {
-                                       qlsts.append( -s.op(j-1));
-                               }
-                       }
-                       if (qlstx.nops() > 0) {
-                               buffer *= G_numeric(qlstx, qlsts, 1/q);
-                       }
-                       lst plstx;
-                       lst plsts;
-                       for (int j=r+1; j<=size; ++j) {
-                               plstx.append(newx.op(j-1));
-                               plsts.append(s.op(j-1));
-                       }
-                       if (plstx.nops() > 0) {
-                               buffer *= G_numeric(plstx, plsts, 1/p);
-                       }
-                       result += buffer;
+               std::vector<cln::cl_N> plstx;
+               std::vector<int> plsts;
+               for (std::size_t j = r+1; j <= size; ++j) {
+                       plstx.push_back(x[j-1]);
+                       plsts.push_back(s[j-1]);
                }
-               return result;
+               if (plstx.size() > 0) {
+                       buffer = buffer*G_numeric(plstx, plsts, 1/p);
+               }
+               result = result + buffer;
        }
-       
-       // convergence transformation
-       if (need_trafo) {
-
-               // sort (|x|<->position) to determine indices
-               std::multimap<ex,int> sortmap;
-               int size = 0;
-               for (int i=0; i<x.nops(); ++i) {
-                       if (!x[i].is_zero()) {
-                               sortmap.insert(std::pair<ex,int>(abs(x[i]), i));
-                               ++size;
-                       }
-               }
-               // include upper limit (scale)
-               sortmap.insert(std::pair<ex,int>(abs(y), x.nops()));
-
-               // generate missing dummy-symbols
-               int i = 1;
-               // holding dummy-symbols for the G/Li transformations
-               exvector gsyms;
-               gsyms.push_back(symbol("GSYMS_ERROR"));
-               ex lastentry;
-               for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
-                       if (it != sortmap.begin()) {
-                               if (it->second < x.nops()) {
-                                       if (x[it->second] == lastentry) {
-                                               gsyms.push_back(gsyms.back());
-                                               continue;
-                                       }
-                               } else {
-                                       if (y == lastentry) {
-                                               gsyms.push_back(gsyms.back());
-                                               continue;
-                                       }
+       return result;
+}
+
+// convergence transformation, used for numerical evaluation of G function.
+// the parameter x, s and y must only contain numerics
+static cln::cl_N
+G_do_trafo(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
+          const cln::cl_N& y)
+{
+       // sort (|x|<->position) to determine indices
+       typedef std::multimap<cln::cl_R, std::size_t> sortmap_t;
+       sortmap_t sortmap;
+       std::size_t size = 0;
+       for (std::size_t i = 0; i < x.size(); ++i) {
+               if (!zerop(x[i])) {
+                       sortmap.insert(std::make_pair(abs(x[i]), i));
+                       ++size;
+               }
+       }
+       // include upper limit (scale)
+       sortmap.insert(std::make_pair(abs(y), x.size()));
+
+       // generate missing dummy-symbols
+       int i = 1;
+       // holding dummy-symbols for the G/Li transformations
+       exvector gsyms;
+       gsyms.push_back(symbol("GSYMS_ERROR"));
+       cln::cl_N lastentry(0);
+       for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
+               if (it != sortmap.begin()) {
+                       if (it->second < x.size()) {
+                               if (x[it->second] == lastentry) {
+                                       gsyms.push_back(gsyms.back());
+                                       continue;
                                }
-                       }
-                       std::ostringstream os;
-                       os << "a" << i;
-                       gsyms.push_back(symbol(os.str()));
-                       ++i;
-                       if (it->second < x.nops()) {
-                               lastentry = x[it->second];
                        } else {
-                               lastentry = y;
+                               if (y == lastentry) {
+                                       gsyms.push_back(gsyms.back());
+                                       continue;
+                               }
                        }
                }
+               std::ostringstream os;
+               os << "a" << i;
+               gsyms.push_back(symbol(os.str()));
+               ++i;
+               if (it->second < x.size()) {
+                       lastentry = x[it->second];
+               } else {
+                       lastentry = y;
+               }
+       }
 
-               // fill position data according to sorted indices and prepare substitution list
-               Gparameter a(x.nops());
-               lst subslst;
-               int pos = 1;
-               int scale;
-               for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
-                       if (it->second < x.nops()) {
-                               if (s[it->second] > 0) {
-                                       a[it->second] = pos;
-                               } else {
-                                       a[it->second] = -pos;
-                               }
-                               subslst.append(gsyms[pos] == x[it->second]);
+       // fill position data according to sorted indices and prepare substitution list
+       Gparameter a(x.size());
+       exmap subslst;
+       std::size_t pos = 1;
+       int scale = pos;
+       for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
+               if (it->second < x.size()) {
+                       if (s[it->second] > 0) {
+                               a[it->second] = pos;
                        } else {
-                               scale = pos;
-                               subslst.append(gsyms[pos] == y);
+                               a[it->second] = -int(pos);
                        }
-                       ++pos;
+                       subslst[gsyms[pos]] = numeric(x[it->second]);
+               } else {
+                       scale = pos;
+                       subslst[gsyms[pos]] = numeric(y);
                }
+               ++pos;
+       }
 
-               // do transformation
-               Gparameter pendint;
-               ex result = G_transform(pendint, a, scale, gsyms);
-               // replace dummy symbols with their values
-               result = result.eval().expand();
-               result = result.subs(subslst).evalf();
-               
-               return result;
+       // do transformation
+       Gparameter pendint;
+       ex result = G_transform(pendint, a, scale, gsyms);
+       // replace dummy symbols with their values
+       result = result.eval().expand();
+       result = result.subs(subslst).evalf();
+       if (!is_a<numeric>(result))
+               throw std::logic_error("G_do_trafo: G_transform returned non-numeric result");
+       
+       cln::cl_N ret = ex_to<numeric>(result).to_cl_N();
+       return ret;
+}
+
+// handles the transformations and the numerical evaluation of G
+// the parameter x, s and y must only contain numerics
+static cln::cl_N
+G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
+         const cln::cl_N& y)
+{
+       // check for convergence and necessary accelerations
+       bool need_trafo = false;
+       bool need_hoelder = false;
+       std::size_t depth = 0;
+       for (std::size_t i = 0; i < x.size(); ++i) {
+               if (!zerop(x[i])) {
+                       ++depth;
+                       const cln::cl_N x_y = abs(x[i]) - y;
+                       if (instanceof(x_y, cln::cl_R_ring) &&
+                           realpart(x_y) < cln::least_negative_float(cln::float_format(Digits - 2)))
+                               need_trafo = true;
+
+                       if (abs(abs(x[i]/y) - 1) < 0.01)
+                               need_hoelder = true;
+               }
        }
+       if (zerop(x[x.size() - 1]))
+               need_trafo = true;
+
+       if (depth == 1 && x.size() == 2 && !need_trafo)
+               return - Li_projection(2, y/x[1], cln::float_format(Digits));
+       
+       // do acceleration transformation (hoelder convolution [BBB])
+       if (need_hoelder)
+               return G_do_hoelder(x, s, y);
+       
+       // convergence transformation
+       if (need_trafo)
+               return G_do_trafo(x, s, y);
 
        // do summation
-       lst newx;
-       lst m;
+       std::vector<cln::cl_N> newx;
+       newx.reserve(x.size());
+       std::vector<int> m;
+       m.reserve(x.size());
        int mcount = 1;
-       ex sign = 1;
-       ex factor = y;
-       for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
-               if ((*it).is_zero()) {
+       int sign = 1;
+       cln::cl_N factor = y;
+       for (std::size_t i = 0; i < x.size(); ++i) {
+               if (zerop(x[i])) {
                        ++mcount;
                } else {
-                       newx.append(factor / (*it));
-                       factor = *it;
-                       m.append(mcount);
+                       newx.push_back(factor/x[i]);
+                       factor = x[i];
+                       m.push_back(mcount);
                        mcount = 1;
                        sign = -sign;
                }
        }
 
-       return sign * numeric(mLi_do_summation(m, newx));
+       return sign*multipleLi_do_sum(m, newx);
 }
 
 
 ex mLi_numeric(const lst& m, const lst& x)
 {
        // let G_numeric do the transformation
-       lst newx;
-       lst s;
-       ex factor = 1;
+       std::vector<cln::cl_N> newx;
+       newx.reserve(x.nops());
+       std::vector<int> s;
+       s.reserve(x.nops());
+       cln::cl_N factor(1);
        for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
                for (int i = 1; i < *itm; ++i) {
-                       newx.append(0);
-                       s.append(1);
+                       newx.push_back(cln::cl_N(0));
+                       s.push_back(1);
+               }
+               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)));
 }
 
 
@@ -1210,7 +1238,8 @@ static ex G2_evalf(const ex& x_, const ex& y)
        if (x.op(0) == y) {
                return G(x_, y).hold();
        }
-       lst s;
+       std::vector<int> s;
+       s.reserve(x.nops());
        bool all_zero = true;
        for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
                if (!(*it).info(info_flags::numeric)) {
@@ -1220,16 +1249,21 @@ static ex G2_evalf(const ex& x_, const ex& y)
                        all_zero = false;
                }
                if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
-                       s.append(-1);
+                       s.push_back(-1);
                }
                else {
-                       s.append(+1);
+                       s.push_back(1);
                }
        }
        if (all_zero) {
                return pow(log(y), x.nops()) / factorial(x.nops());
        }
-       return G_numeric(x, s, y);
+       std::vector<cln::cl_N> xv;
+       xv.reserve(x.nops());
+       for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+               xv.push_back(ex_to<numeric>(*it).to_cl_N());
+       cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
+       return numeric(result);
 }
 
 
@@ -1247,7 +1281,8 @@ static ex G2_eval(const ex& x_, const ex& y)
        if (x.op(0) == y) {
                return G(x_, y).hold();
        }
-       lst s;
+       std::vector<int> s;
+       s.reserve(x.nops());
        bool all_zero = true;
        bool crational = true;
        for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
@@ -1261,10 +1296,10 @@ static ex G2_eval(const ex& x_, const ex& y)
                        all_zero = false;
                }
                if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
-                       s.append(-1);
+                       s.push_back(-1);
                }
                else {
-                       s.append(+1);
+                       s.push_back(+1);
                }
        }
        if (all_zero) {
@@ -1276,7 +1311,12 @@ static ex G2_eval(const ex& x_, const ex& y)
        if (crational) {
                return G(x_, y).hold();
        }
-       return G_numeric(x, s, y);
+       std::vector<cln::cl_N> xv;
+       xv.reserve(x.nops());
+       for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+               xv.push_back(ex_to<numeric>(*it).to_cl_N());
+       cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
+       return numeric(result);
 }
 
 
@@ -1306,7 +1346,8 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
        if (x.op(0) == y) {
                return G(x_, s_, y).hold();
        }
-       lst sn;
+       std::vector<int> sn;
+       sn.reserve(s.nops());
        bool all_zero = true;
        for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
                if (!(*itx).info(info_flags::numeric)) {
@@ -1319,26 +1360,35 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
                        all_zero = false;
                }
                if ( ex_to<numeric>(*itx).is_real() ) {
-                       if ( *its >= 0 ) {
-                               sn.append(+1);
-                       }
-                       else {
-                               sn.append(-1);
+                       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.append(+1);
+                               sn.push_back(1);
                        }
                        else {
-                               sn.append(-1);
+                               sn.push_back(-1);
                        }
                }
        }
        if (all_zero) {
                return pow(log(y), x.nops()) / factorial(x.nops());
        }
-       return G_numeric(x, sn, y);
+       std::vector<cln::cl_N> xn;
+       xn.reserve(x.nops());
+       for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+               xn.push_back(ex_to<numeric>(*it).to_cl_N());
+       cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
+       return numeric(result);
 }
 
 
@@ -1360,7 +1410,8 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
        if (x.op(0) == y) {
                return G(x_, s_, y).hold();
        }
-       lst sn;
+       std::vector<int> sn;
+       sn.reserve(s.nops());
        bool all_zero = true;
        bool crational = true;
        for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
@@ -1377,19 +1428,23 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
                        all_zero = false;
                }
                if ( ex_to<numeric>(*itx).is_real() ) {
-                       if ( *its >= 0 ) {
-                               sn.append(+1);
-                       }
-                       else {
-                               sn.append(-1);
+                       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.append(+1);
+                               sn.push_back(1);
                        }
                        else {
-                               sn.append(-1);
+                               sn.push_back(-1);
                        }
                }
        }
@@ -1402,7 +1457,12 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
        if (crational) {
                return G(x_, s_, y).hold();
        }
-       return G_numeric(x, sn, y);
+       std::vector<cln::cl_N> xn;
+       xn.reserve(x.nops());
+       for (lst::const_iterator it = x.begin(); it != x.end(); ++it)
+               xn.push_back(ex_to<numeric>(*it).to_cl_N());
+       cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
+       return numeric(result);
 }
 
 
@@ -1644,7 +1704,7 @@ static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c)
        } else {
                x = lst(x_);
        }
-       c.s << "\\mbox{Li}_{";
+       c.s << "\\mathrm{Li}_{";
        lst::const_iterator itm = m.begin();
        (*itm).print(c);
        itm++;
@@ -2160,7 +2220,7 @@ static ex S_deriv(const ex& n, const ex& p, const ex& x, unsigned deriv_param)
 
 static void S_print_latex(const ex& n, const ex& p, const ex& x, const print_context& c)
 {
-       c.s << "\\mbox{S}_{";
+       c.s << "\\mathrm{S}_{";
        n.print(c);
        c.s << ",";
        p.print(c);
@@ -2241,7 +2301,7 @@ bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
                }
        }
        if (has_negative_parameters) {
-               for (int i=0; i<m.nops(); i++) {
+               for (std::size_t i=0; i<m.nops(); i++) {
                        if (m.op(i) < 0) {
                                m.let_op(i) = -m.op(i);
                                s.append(-1);
@@ -2281,7 +2341,7 @@ struct map_trafo_H_convert_to_Li : public map_function
                                        s.let_op(0) = s.op(0) * arg;
                                        return pf * Li(m, s).hold();
                                } else {
-                                       for (int i=0; i<m.nops(); i++) {
+                                       for (std::size_t i=0; i<m.nops(); i++) {
                                                s.append(1);
                                        }
                                        s.let_op(0) = s.op(0) * arg;
@@ -2363,7 +2423,7 @@ struct map_trafo_H_reduce_trailing_zeros : public map_function
                                        
                                        //
                                        parameter.remove_last();
-                                       int lastentry = parameter.nops();
+                                       std::size_t lastentry = parameter.nops();
                                        while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
                                                lastentry--;
                                        }
@@ -2452,9 +2512,9 @@ ex trafo_H_mult(const ex& h1, const ex& h2)
                        hlong = h2.op(0).op(0);
                }
        }
-       for (int i=0; i<=hlong.nops(); i++) {
+       for (std::size_t i=0; i<=hlong.nops(); i++) {
                lst newparameter;
-               int j=0;
+               std::size_t j=0;
                for (; j<i; j++) {
                        newparameter.append(hlong[j]);
                }
@@ -2482,7 +2542,7 @@ struct map_trafo_H_mult : public map_function
                        ex result = 1;
                        ex firstH;
                        lst Hlst;
-                       for (int pos=0; pos<e.nops(); pos++) {
+                       for (std::size_t pos=0; pos<e.nops(); pos++) {
                                if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
                                        std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
                                        if (name == "H") {
@@ -2516,7 +2576,7 @@ struct map_trafo_H_mult : public map_function
                        if (Hlst.nops() > 0) {
                                ex buffer = trafo_H_mult(firstH, Hlst.op(0));
                                result *= buffer;
-                               for (int i=1; i<Hlst.nops(); i++) {
+                               for (std::size_t i=1; i<Hlst.nops(); i++) {
                                        result *= Hlst.op(i);
                                }
                                result = result.expand();
@@ -2544,7 +2604,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") {
@@ -2559,7 +2619,7 @@ 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());
        }
 }
 
@@ -2576,7 +2636,7 @@ ex trafo_H_prepend_one(const ex& e, const ex& arg)
        if (name == "H") {
                h = e;
        } else {
-               for (int i=0; i<e.nops(); i++) {
+               for (std::size_t i=0; i<e.nops(); i++) {
                        if (is_a<function>(e.op(i))) {
                                std::string name = ex_to<function>(e.op(i)).get_name();
                                if (name == "H") {
@@ -2590,7 +2650,7 @@ ex trafo_H_prepend_one(const ex& e, const ex& arg)
                newparameter.prepend(1);
                return e.subs(h == H(newparameter, h.op(1)).hold());
        } else {
-               return e * H(lst(1),1-arg).hold();
+               return e * H(lst(ex(1)),1-arg).hold();
        }
 }
 
@@ -2607,7 +2667,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") {
@@ -2622,8 +2682,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();
        }
 }
 
@@ -2640,7 +2700,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") {
@@ -2654,7 +2714,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();
        }
 }
 
@@ -2671,7 +2731,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") {
@@ -2685,7 +2745,7 @@ 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();
        }
 }
 
@@ -2709,7 +2769,7 @@ struct map_trafo_H_1mx : public map_function
                                // special cases if all parameters are either 0, 1 or -1
                                bool allthesame = true;
                                if (parameter.op(0) == 0) {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != 0) {
                                                        allthesame = false;
                                                        break;
@@ -2725,7 +2785,7 @@ struct map_trafo_H_1mx : public map_function
                                } else if (parameter.op(0) == -1) {
                                        throw std::runtime_error("map_trafo_H_1mx: cannot handle weights equal -1!");
                                } else {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != 1) {
                                                        allthesame = false;
                                                        break;
@@ -2751,7 +2811,7 @@ struct map_trafo_H_1mx : public map_function
                                        map_trafo_H_1mx recursion;
                                        ex buffer = recursion(H(newparameter, arg).hold());
                                        if (is_a<add>(buffer)) {
-                                               for (int i=0; i<buffer.nops(); i++) {
+                                               for (std::size_t i = 0; i < buffer.nops(); i++) {
                                                        res -= trafo_H_prepend_one(buffer.op(i), arg);
                                                }
                                        } else {
@@ -2764,14 +2824,14 @@ struct map_trafo_H_1mx : public map_function
                                        // leading one
                                        map_trafo_H_1mx recursion;
                                        map_trafo_H_mult unify;
-                                       ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
-                                       int firstzero = 0;
+                                       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]);
                                                }
@@ -2810,7 +2870,7 @@ struct map_trafo_H_1overx : public map_function
                                // special cases if all parameters are either 0, 1 or -1
                                bool allthesame = true;
                                if (parameter.op(0) == 0) {
-                                       for (int i=1; i<parameter.nops(); i++) {
+                                       for (std::size_t i = 1; i < parameter.nops(); i++) {
                                                if (parameter.op(i) != 0) {
                                                        allthesame = false;
                                                        break;
@@ -2820,7 +2880,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;
@@ -2828,11 +2888,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;
@@ -2840,7 +2900,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());
                                        }
                                }
@@ -2855,7 +2915,7 @@ struct map_trafo_H_1overx : public map_function
                                        map_trafo_H_1overx recursion;
                                        ex buffer = recursion(H(newparameter, arg).hold());
                                        if (is_a<add>(buffer)) {
-                                               for (int i=0; i<buffer.nops(); i++) {
+                                               for (std::size_t i = 0; i < buffer.nops(); i++) {
                                                        res += trafo_H_1tx_prepend_zero(buffer.op(i), arg);
                                                }
                                        } else {
@@ -2870,7 +2930,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 {
@@ -2883,14 +2943,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]);
                                                }
@@ -2931,7 +2991,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;
@@ -2939,11 +2999,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;
@@ -2951,11 +3011,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;
@@ -2963,7 +3023,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());
                                        }
                                }
@@ -2978,7 +3038,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                        map_trafo_H_1mxt1px recursion;
                                        ex buffer = recursion(H(newparameter, arg).hold());
                                        if (is_a<add>(buffer)) {
-                                               for (int i=0; i<buffer.nops(); i++) {
+                                               for (std::size_t i = 0; i < buffer.nops(); i++) {
                                                        res -= trafo_H_1mxt1px_prepend_one(buffer.op(i), arg) + trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
                                                }
                                        } else {
@@ -2993,7 +3053,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 {
@@ -3006,14 +3066,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]);
                                                }
@@ -3087,7 +3147,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                        }
                }
 
-               for (int i=0; i<x1.nops(); i++) {
+               for (std::size_t i = 0; i < x1.nops(); i++) {
                        if (!x1.op(i).info(info_flags::integer)) {
                                return H(x1, x2).hold();
                        }
@@ -3159,7 +3219,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                // ensure that the realpart of the argument is positive
                if (cln::realpart(x) < 0) {
                        x = -x;
-                       for (int i=0; i<m.nops(); i++) {
+                       for (std::size_t i = 0; i < m.nops(); i++) {
                                if (m.op(i) != 0) {
                                        m.let_op(i) = -m.op(i);
                                        res *= -1;
@@ -3366,7 +3426,7 @@ static void H_print_latex(const ex& m_, const ex& x, const print_context& c)
        } else {
                m = lst(m_);
        }
-       c.s << "\\mbox{H}_{";
+       c.s << "\\mathrm{H}_{";
        lst::const_iterator itm = m.begin();
        (*itm).print(c);
        itm++;
@@ -3442,7 +3502,7 @@ static void initcX(std::vector<cln::cl_N>& crX,
        int Sm = 0;
        int Smp1 = 0;
        std::vector<std::vector<cln::cl_N> > crG(s.size() - 1, std::vector<cln::cl_N>(L2 + 1));
-       for (int m=0; m < s.size() - 1; m++) {
+       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++)
@@ -3668,7 +3728,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];