]> www.ginac.de Git - ginac.git/blobdiff - ginac/inifcns_nstdsums.cpp
Avoid x^0 and Order(x^0) terms together in series expansion.
[ginac.git] / ginac / inifcns_nstdsums.cpp
index 3284ff698d267757e47e6b23b1f79a563e121283..6c17e32e729bedea127eb09573dbb8ff400f2138 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-2009 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
@@ -102,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;
@@ -346,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);
                        }
@@ -368,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);
@@ -798,12 +803,12 @@ 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 exvector& gsyms);
+            const exvector& gsyms, bool flag_trailing_zeros_only);
 
 
 // G transformation [VSW]
 ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale,
-              const exvector& gsyms)
+              const exvector& gsyms, bool flag_trailing_zeros_only)
 {
        // main recursion routine
        //
@@ -843,18 +848,18 @@ 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, gsyms) * G_transform(pendint, new_a, scale, gsyms);
+               result += G_eval1(0, scale, gsyms) * G_transform(pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
                for (Gparameter::const_iterator it = a.begin(); it != firstzero; ++it) {
                        Gparameter new_a(a.begin(), it);
                        new_a.push_back(0);
                        new_a.insert(new_a.end(), it, a.end()-1);
-                       result -= G_transform(pendint, new_a, scale, gsyms);
+                       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(), gsyms)*
@@ -881,10 +886,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, gsyms)*
-                       G_transform(empty, a1, scale, gsyms);
+               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, gsyms);
+               result -= shuffle_G(empty, a1, a2, pendint, a, scale, gsyms, flag_trailing_zeros_only);
                return result;
        }
 
@@ -895,7 +900,7 @@ 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, gsyms);
+       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(), gsyms);
@@ -909,31 +914,31 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale,
                new_pendint.push_back(*changeit);
                result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint),
                                           new_pendint.front(), gsyms)*
-                       G_transform(empty, new_a, scale, 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, gsyms);
+               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(), gsyms)*
-                       G_transform(empty, new_a, scale, gsyms);
+                         G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only);
                *changeit = *min_it;
-               result -= G_transform(new_pendint, new_a, scale, gsyms);
+               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(), gsyms)*
-                       G_transform(empty, new_a, scale, 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(), gsyms)*
-                       G_transform(empty, new_a, scale, gsyms);
+                         G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only);
                *changeit = *min_it;
-               result += G_transform(new_pendint, new_a, scale, gsyms);
+               result += G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
        }
        return result;
 }
@@ -943,27 +948,27 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale,
 // 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 exvector& gsyms) 
+            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, gsyms);
+               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, gsyms);
+               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, gsyms);
+               return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms, flag_trailing_zeros_only);
        }
 
        Gparameter a1_removed(a1.begin()+1,a1.end());
@@ -975,8 +980,8 @@ ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2
        a01.push_back( a1[0] );
        a02.push_back( a2[0] );
 
-       return shuffle_G(a01, a1_removed, a2, pendint, a_old, scale, gsyms)
-            + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms);
+       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
@@ -1015,9 +1020,8 @@ G_do_hoelder(std::vector<cln::cl_N> x, /* yes, it's passed by value */
                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]);
+                       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]);
                        }
@@ -1039,24 +1043,43 @@ G_do_hoelder(std::vector<cln::cl_N> x, /* yes, it's passed by value */
        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, 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)
+           const cln::cl_N& y, bool flag_trailing_zeros_only)
 {
        // sort (|x|<->position) to determine indices
-       typedef std::multimap<cln::cl_R, std::size_t> sortmap_t;
+       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(abs(x[i]), i));
+                       sortmap.insert(std::make_pair(x[i], i));
                        ++size;
                }
        }
        // include upper limit (scale)
-       sortmap.insert(std::make_pair(abs(y), x.size()));
+       sortmap.insert(std::make_pair(y, x.size()));
 
        // generate missing dummy-symbols
        int i = 1;
@@ -1111,9 +1134,9 @@ G_do_trafo(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
 
        // do transformation
        Gparameter pendint;
-       ex result = G_transform(pendint, a, scale, gsyms);
+       ex result = G_transform(pendint, a, scale, gsyms, flag_trailing_zeros_only);
        // replace dummy symbols with their values
-       result = result.eval().expand();
+       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");
@@ -1131,6 +1154,7 @@ G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
        // 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 (std::size_t i = 0; i < x.size(); ++i) {
                if (!zerop(x[i])) {
@@ -1144,19 +1168,21 @@ G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
                                need_hoelder = true;
                }
        }
-       if (zerop(x[x.size() - 1]))
+       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)
+       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);
+               return G_do_trafo(x, s, y, have_trailing_zero);
 
        // do summation
        std::vector<cln::cl_N> newx;
@@ -1196,9 +1222,14 @@ ex mLi_numeric(const lst& m, const lst& x)
                        s.push_back(1);
                }
                const cln::cl_N xi = ex_to<numeric>(*itx).to_cl_N();
-               newx.push_back(factor/xi);
                factor = factor/xi;
-               s.push_back(1);
+               newx.push_back(factor);
+               if ( !instanceof(factor, cln::cl_R_ring) && imagpart(factor) < 0 ) {
+                       s.push_back(-1);
+               }
+               else {
+                       s.push_back(1);
+               }
        }
        return numeric(cln::cl_N(1 & m.nops() ? - 1 : 1)*G_numeric(newx, s, cln::cl_N(1)));
 }
@@ -1218,10 +1249,10 @@ 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;
        }
@@ -1261,10 +1292,10 @@ 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;
        }
@@ -1310,10 +1341,10 @@ static ex G2_eval(const ex& x_, const ex& y)
 }
 
 
+// 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).
@@ -1322,11 +1353,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();
        }
@@ -1350,12 +1381,16 @@ 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 ) {
+                       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 {
-                               sn.push_back(-1);
-                       }
                }
                else {
                        if ( ex_to<numeric>(*itx).imag() > 0 ) {
@@ -1382,11 +1417,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();
        }
@@ -1414,12 +1449,16 @@ 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 ) {
+                       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 {
-                               sn.push_back(-1);
-                       }
                }
                else {
                        if ( ex_to<numeric>(*itx).imag() > 0 ) {
@@ -1448,10 +1487,12 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
 }
 
 
+// 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).
@@ -1559,7 +1600,17 @@ static ex Li_eval(const ex& m_, const ex& x_)
                                }
                        }
                        if (is_zeta) {
-                               return zeta(m_,x_);
+                               lst newx;
+                               for (lst::const_iterator itx = x.begin(); itx != x.end(); ++itx) {
+                                       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;
@@ -1611,9 +1662,8 @@ static ex Li_series(const ex& m, const ex& x, const relational& rel, int order,
 {
        if (is_a<lst>(m) || is_a<lst>(x)) {
                // multiple polylog
-               epvector seq;
-               seq.push_back(expair(Li(m, x), 0));
-               return pseries(rel, seq);
+               epvector seq { expair(Li(m, x), 0) };
+               return pseries(rel, std::move(seq));
        }
        
        // classical polylog
@@ -1629,9 +1679,8 @@ static ex Li_series(const ex& m, const ex& x, const relational& rel, int order,
                        // substitute the argument's series expansion
                        ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
                        // maybe that was terminating, so add a proper order term
-                       epvector nseq;
-                       nseq.push_back(expair(Order(_ex1), order));
-                       ser += pseries(rel, nseq);
+                       epvector nseq { expair(Order(_ex1), order) };
+                       ser += pseries(rel, std::move(nseq));
                        // reexpanding it will collapse the series again
                        return ser.series(rel, order);
                }
@@ -1678,13 +1727,13 @@ 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 << "\\mathrm{Li}_{";
        lst::const_iterator itm = m.begin();
@@ -1730,7 +1779,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]
 
@@ -2026,7 +2075,9 @@ const cln::cl_N S_num(int n, int p, const cln::cl_N& x)
                prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 
        // [Kol] (5.3)
-       if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95))) {
+       // 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);
@@ -2067,6 +2118,16 @@ const cln::cl_N S_num(int n, int p, const cln::cl_N& 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);
        }
@@ -2114,7 +2175,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);
                        }
@@ -2172,9 +2233,8 @@ static ex S_series(const ex& n, const ex& p, const ex& x, const relational& rel,
                        // substitute the argument's series expansion
                        ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
                        // maybe that was terminating, so add a proper order term
-                       epvector nseq;
-                       nseq.push_back(expair(Order(_ex1), order));
-                       ser += pseries(rel, nseq);
+                       epvector nseq { expair(Order(_ex1), order) };
+                       ser += pseries(rel, std::move(nseq));
                        // reexpanding it will collapse the series again
                        return ser.series(rel, order);
                }
@@ -2300,7 +2360,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);
@@ -2312,7 +2372,7 @@ struct map_trafo_H_convert_to_Li : 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);
 
@@ -2339,7 +2399,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);
@@ -2351,7 +2411,7 @@ struct map_trafo_H_convert_to_zeta : 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)};
                                }
 
                                lst m;
@@ -2372,7 +2432,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);
@@ -2384,7 +2444,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) {
@@ -2464,7 +2524,12 @@ lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf)
        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++;
@@ -2491,7 +2556,7 @@ 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 (std::size_t i=0; i<=hlong.nops(); i++) {
@@ -2513,7 +2578,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);
@@ -2601,7 +2666,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());
        }
 }
 
@@ -2632,7 +2697,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();
        }
 }
 
@@ -2664,8 +2729,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();
        }
 }
 
@@ -2696,7 +2761,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();
        }
 }
 
@@ -2727,7 +2792,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();
        }
 }
 
@@ -2735,7 +2800,7 @@ ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
 // do x -> 1-x transformation
 struct map_trafo_H_1mx : public map_function
 {
-       ex operator()(const ex& e)
+       ex operator()(const ex& e) override
        {
                if (is_a<add>(e) || is_a<mul>(e)) {
                        return e.map(*this);
@@ -2806,7 +2871,7 @@ 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();
+                                       ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
                                        std::size_t firstzero = 0;
                                        while (parameter.op(firstzero) == 1) {
                                                firstzero++;
@@ -2836,7 +2901,7 @@ struct map_trafo_H_1mx : public map_function
 // 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);
@@ -2870,7 +2935,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(), 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 {
@@ -2882,7 +2947,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());
                                        }
                                }
@@ -2925,7 +2990,7 @@ 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();
+                                       ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
                                        std::size_t firstzero = 0;
                                        while (parameter.op(firstzero) == 1) {
                                                firstzero++;
@@ -2957,7 +3022,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);
@@ -2981,7 +3046,7 @@ 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) {
@@ -2993,7 +3058,7 @@ 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 {
@@ -3005,7 +3070,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());
                                        }
                                }
@@ -3048,7 +3113,7 @@ 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();
+                                       ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
                                        std::size_t firstzero = 0;
                                        while (parameter.op(firstzero) == 1) {
                                                firstzero++;
@@ -3212,7 +3277,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                // 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 {
@@ -3227,7 +3292,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                if (cln::abs(x-9.53) <= 9.47) {
                        // x -> (1-x)/(1+x)
                        map_trafo_H_1mxt1px trafo;
-                       res *= trafo(H(m, xtemp));
+                       res *= trafo(H(m, xtemp).hold());
                } else {
                        // x -> 1-x
                        if (has_minus_one) {
@@ -3235,7 +3300,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                                return filter(H(m, numeric(x)).hold()).evalf();
                        }
                        map_trafo_H_1mx trafo;
-                       res *= trafo(H(m, xtemp));
+                       res *= trafo(H(m, xtemp).hold());
                }
 
                return res.subs(xtemp == numeric(x)).evalf();
@@ -3251,7 +3316,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;
@@ -3362,9 +3427,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));
 }
 
 
@@ -3378,7 +3442,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) {
@@ -3406,7 +3470,7 @@ 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 << "\\mathrm{H}_{";
        lst::const_iterator itm = m.begin();
@@ -3439,7 +3503,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()));
        }
 }
 
@@ -3483,7 +3547,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));
+       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];
@@ -3523,12 +3587,12 @@ static cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk,
 
 
 // [Cra] section 4
-static void calc_f(std::vector<std::vector<cln::cl_N> >& f_kj,
+static void calc_f(std::vector<std::vector<cln::cl_N>>& f_kj,
                   const int maxr, const int 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();
+       std::vector<std::vector<cln::cl_N>>::iterator it = f_kj.begin();
        cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
        
        t0 = cln::exp(-lambda);
@@ -3552,7 +3616,7 @@ static void calc_f(std::vector<std::vector<cln::cl_N> >& f_kj,
 
 // [Cra] (3.1)
 static cln::cl_N crandall_Z(const std::vector<int>& s,
-                           const std::vector<std::vector<cln::cl_N> >& f_kj)
+                           const std::vector<std::vector<cln::cl_N>>& f_kj)
 {
        const int j = s.size();
 
@@ -3629,7 +3693,7 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
                }
        }
 
-       std::vector<std::vector<cln::cl_N> > f_kj(L1);
+       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);
@@ -4018,13 +4082,13 @@ 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();