X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=6c17e32e729bedea127eb09573dbb8ff400f2138;hp=d85044d04bd1f311b2d26ad5ab8ccf0572d1ebe7;hb=f8c2455fbbd34dbbfb72ac96eec0f45fd453bade;hpb=b8150271ff497c2ef0c2a8748f5f53b4f1bab7c8 diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index d85044d0..6c17e32e 100644 --- a/ginac/inifcns_nstdsums.cpp +++ b/ginac/inifcns_nstdsums.cpp @@ -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-2008 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 @@ -64,11 +64,6 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ -#include -#include -#include -#include - #include "inifcns.h" #include "add.h" @@ -84,6 +79,10 @@ #include "utils.h" #include "wildcard.h" +#include +#include +#include +#include namespace GiNaC { @@ -103,7 +102,7 @@ namespace { // lookup table for factors built from Bernoulli numbers // see fill_Xn() -std::vector > Xn; +std::vector> Xn; // initial size of Xn that should suffice for 32bit machines (must be even) const int xninitsizestep = 26; int xninitsize = xninitsizestep; @@ -210,7 +209,7 @@ void double_Xn() } } // X_n - for (int n=2; n 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& s, const std::vector m_int; - std::vector x_cln; - for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) { - m_int.push_back(ex_to(*itm).to_int()); - x_cln.push_back(ex_to(*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; } @@ -809,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 // @@ -854,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)* @@ -892,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; } @@ -906,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); @@ -920,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; } @@ -954,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()); @@ -986,203 +980,258 @@ ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2 a01.push_back( a1[0] ); a02.push_back( a2[0] ); - return shuffle_G(a01, a1_removed, a2, pendint, a_old, scale, 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 // 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& x, const std::vector& 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 x, /* yes, it's passed by value */ + const std::vector& 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 qlstx; + std::vector qlsts; + for (std::size_t j = r; j >= 1; --j) { + qlstx.push_back(cln::cl_N(1) - x[j-1]); + if (instanceof(x[j-1], cln::cl_R_ring) && realpart(x[j-1]) > 1) { + qlsts.push_back(1); + } else { + qlsts.push_back(-s[j-1]); } } - } - if (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 plstx; + std::vector 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 sortmap; - int size = 0; - for (int i=0; i(abs(x[i]), i)); - ++size; - } - } - // include upper limit (scale) - sortmap.insert(std::pair(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::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; +} + +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& x, const std::vector& s, + const cln::cl_N& y, bool flag_trailing_zeros_only) +{ + // sort (|x|<->position) to determine indices + typedef std::multimap sortmap_t; + sortmap_t sortmap; + std::size_t size = 0; + for (std::size_t i = 0; i < x.size(); ++i) { + if (!zerop(x[i])) { + sortmap.insert(std::make_pair(x[i], i)); + ++size; + } + } + // include upper limit (scale) + sortmap.insert(std::make_pair(y, x.size())); + + // generate missing dummy-symbols + int i = 1; + // holding dummy-symbols for the G/Li transformations + exvector gsyms; + gsyms.push_back(symbol("GSYMS_ERROR")); + cln::cl_N lastentry(0); + for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) { + if (it != sortmap.begin()) { + if (it->second < x.size()) { + if (x[it->second] == lastentry) { + gsyms.push_back(gsyms.back()); + continue; } - } - std::ostringstream os; - os << "a" << i; - gsyms.push_back(symbol(os.str())); - ++i; - if (it->second < x.nops()) { - lastentry = x[it->second]; } else { - lastentry = y; + if (y == lastentry) { + gsyms.push_back(gsyms.back()); + continue; + } } } + std::ostringstream os; + os << "a" << i; + gsyms.push_back(symbol(os.str())); + ++i; + if (it->second < x.size()) { + lastentry = x[it->second]; + } else { + lastentry = y; + } + } - // fill position data according to sorted indices and prepare substitution list - Gparameter a(x.nops()); - lst subslst; - int pos = 1; - int scale; - for (std::multimap::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, flag_trailing_zeros_only); + // replace dummy symbols with their values + result = result.expand(); + result = result.subs(subslst).evalf(); + if (!is_a(result)) + throw std::logic_error("G_do_trafo: G_transform returned non-numeric result"); + + cln::cl_N ret = ex_to(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& x, const std::vector& s, + const cln::cl_N& y) +{ + // check for convergence and necessary accelerations + bool need_trafo = false; + bool need_hoelder = false; + bool have_trailing_zero = false; + std::size_t depth = 0; + for (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.back())) { + have_trailing_zero = true; + need_trafo = true; } + if (depth == 1 && x.size() == 2 && !need_trafo) + return - Li_projection(2, y/x[1], cln::float_format(Digits)); + + // do acceleration transformation (hoelder convolution [BBB]) + if (need_hoelder && !have_trailing_zero) + return G_do_hoelder(x, s, y); + + // convergence transformation + if (need_trafo) + return G_do_trafo(x, s, y, have_trailing_zero); + // do summation - lst newx; - lst m; + std::vector newx; + newx.reserve(x.size()); + std::vector 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 newx; + newx.reserve(x.nops()); + std::vector 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(*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))); } @@ -1200,17 +1249,18 @@ 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(x_) ? ex_to(x_) : lst(x_); + lst x = is_a(x_) ? ex_to(x_) : lst{x_}; if (x.nops() == 0) { return _ex1; } if (x.op(0) == y) { return G(x_, y).hold(); } - lst s; + std::vector 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 +1270,21 @@ static ex G2_evalf(const ex& x_, const ex& y) all_zero = false; } if ( !ex_to(*it).is_real() && ex_to(*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 xv; + xv.reserve(x.nops()); + for (lst::const_iterator it = x.begin(); it != x.end(); ++it) + xv.push_back(ex_to(*it).to_cl_N()); + cln::cl_N result = G_numeric(xv, s, ex_to(y).to_cl_N()); + return numeric(result); } @@ -1237,17 +1292,18 @@ 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(x_) ? ex_to(x_) : lst(x_); + lst x = is_a(x_) ? ex_to(x_) : lst{x_}; if (x.nops() == 0) { return _ex1; } if (x.op(0) == y) { return G(x_, y).hold(); } - lst s; + std::vector 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 +1317,10 @@ static ex G2_eval(const ex& x_, const ex& y) all_zero = false; } if ( !ex_to(*it).is_real() && ex_to(*it).imag() < 0 ) { - s.append(-1); + s.push_back(-1); } else { - s.append(+1); + s.push_back(+1); } } if (all_zero) { @@ -1276,14 +1332,19 @@ static ex G2_eval(const ex& x_, const ex& y) if (crational) { return G(x_, y).hold(); } - return G_numeric(x, s, y); + std::vector xv; + xv.reserve(x.nops()); + for (lst::const_iterator it = x.begin(); it != x.end(); ++it) + xv.push_back(ex_to(*it).to_cl_N()); + cln::cl_N result = G_numeric(xv, s, ex_to(y).to_cl_N()); + return numeric(result); } +// option do_not_evalf_params() removed. unsigned G2_SERIAL::serial = function::register_new(function_options("G", 2). evalf_func(G2_evalf). eval_func(G2_eval). - do_not_evalf_params(). overloaded(2)); //TODO // derivative_func(G2_deriv). @@ -1292,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(x_) ? ex_to(x_) : lst(x_); - lst s = is_a(s_) ? ex_to(s_) : lst(s_); + lst x = is_a(x_) ? ex_to(x_) : lst{x_}; + lst s = is_a(s_) ? ex_to(s_) : lst{s_}; if (x.nops() != s.nops()) { return G(x_, s_, y).hold(); } @@ -1306,7 +1367,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 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 +1381,35 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y) all_zero = false; } if ( ex_to(*itx).is_real() ) { - if ( *its >= 0 ) { - sn.append(+1); - } - else { - sn.append(-1); + if ( ex_to(*itx).is_positive() ) { + if ( *its >= 0 ) { + sn.push_back(1); + } + else { + sn.push_back(-1); + } + } else { + sn.push_back(1); } } else { if ( ex_to(*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 xn; + xn.reserve(x.nops()); + for (lst::const_iterator it = x.begin(); it != x.end(); ++it) + xn.push_back(ex_to(*it).to_cl_N()); + cln::cl_N result = G_numeric(xn, sn, ex_to(y).to_cl_N()); + return numeric(result); } @@ -1346,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(x_) ? ex_to(x_) : lst(x_); - lst s = is_a(s_) ? ex_to(s_) : lst(s_); + lst x = is_a(x_) ? ex_to(x_) : lst{x_}; + lst s = is_a(s_) ? ex_to(s_) : lst{s_}; if (x.nops() != s.nops()) { return G(x_, s_, y).hold(); } @@ -1360,7 +1431,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 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 +1449,23 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y) all_zero = false; } if ( ex_to(*itx).is_real() ) { - if ( *its >= 0 ) { - sn.append(+1); - } - else { - sn.append(-1); + if ( ex_to(*itx).is_positive() ) { + if ( *its >= 0 ) { + sn.push_back(1); + } + else { + sn.push_back(-1); + } + } else { + sn.push_back(1); } } else { if ( ex_to(*itx).imag() > 0 ) { - sn.append(+1); + sn.push_back(1); } else { - sn.append(-1); + sn.push_back(-1); } } } @@ -1402,14 +1478,21 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y) if (crational) { return G(x_, s_, y).hold(); } - return G_numeric(x, sn, y); + std::vector xn; + xn.reserve(x.nops()); + for (lst::const_iterator it = x.begin(); it != x.end(); ++it) + xn.push_back(ex_to(*it).to_cl_N()); + cln::cl_N result = G_numeric(xn, sn, ex_to(y).to_cl_N()); + return numeric(result); } +// option do_not_evalf_params() removed. +// This is safe: in the code above it only matters if s_ > 0 or s_ < 0, +// s_ is allowed to be of floating type. unsigned G3_SERIAL::serial = function::register_new(function_options("G", 3). evalf_func(G3_evalf). eval_func(G3_eval). - do_not_evalf_params(). overloaded(2)); //TODO // derivative_func(G3_deriv). @@ -1517,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; @@ -1569,9 +1662,8 @@ static ex Li_series(const ex& m, const ex& x, const relational& rel, int order, { if (is_a(m) || is_a(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 @@ -1587,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); } @@ -1636,15 +1727,15 @@ static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c) if (is_a(m_)) { m = ex_to(m_); } else { - m = lst(m_); + m = lst{m_}; } lst x; if (is_a(x_)) { x = ex_to(x_); } else { - x = lst(x_); + x = lst{x_}; } - c.s << "\\mbox{Li}_{"; + c.s << "\\mathrm{Li}_{"; lst::const_iterator itm = m.begin(); (*itm).print(c); itm++; @@ -1688,7 +1779,7 @@ namespace { // lookup table for special Euler-Zagier-Sums (used for S_n,p(x)) // see fill_Yn() -std::vector > Yn; +std::vector> Yn; int ynsize = 0; // number of Yn[] int ynlength = 100; // initial length of all Yn[i] @@ -1984,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::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); @@ -2025,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(res).to_cl_N(); + } else { return S_projection(n, p, value, prec); } @@ -2072,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(p).to_int()-1; i>0; i--) { m.append(1); } @@ -2130,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); } @@ -2160,7 +2262,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 +2343,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(e) || is_a(e)) { return e.map(*this); @@ -2270,7 +2372,7 @@ struct map_trafo_H_convert_to_Li : public map_function if (is_a(e.op(0))) { parameter = ex_to(e.op(0)); } else { - parameter = lst(e.op(0)); + parameter = lst{e.op(0)}; } ex arg = e.op(1); @@ -2281,7 +2383,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(e) || is_a(e)) { return e.map(*this); @@ -2309,7 +2411,7 @@ struct map_trafo_H_convert_to_zeta : public map_function if (is_a(e.op(0))) { parameter = ex_to(e.op(0)); } else { - parameter = lst(e.op(0)); + parameter = lst{e.op(0)}; } lst m; @@ -2330,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(e) || is_a(e)) { return e.map(*this); @@ -2342,7 +2444,7 @@ struct map_trafo_H_reduce_trailing_zeros : public map_function if (is_a(e.op(0))) { parameter = ex_to(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) { @@ -2363,7 +2465,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--; } @@ -2422,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++; @@ -2449,12 +2556,12 @@ ex trafo_H_mult(const ex& h1, const ex& h2) if (h2nops > 1) { hlong = ex_to(h2.op(0)); } else { - hlong = h2.op(0).op(0); + hlong = lst{h2.op(0).op(0)}; } } - for (int i=0; i<=hlong.nops(); i++) { + for (std::size_t i=0; i<=hlong.nops(); i++) { lst newparameter; - int j=0; + std::size_t j=0; for (; j(e)) { return e.map(*this); @@ -2482,7 +2589,7 @@ struct map_trafo_H_mult : public map_function ex result = 1; ex firstH; lst Hlst; - for (int pos=0; pos(e.op(pos)) && is_a(e.op(pos).op(0))) { std::string name = ex_to(e.op(pos).op(0)).get_name(); if (name == "H") { @@ -2516,7 +2623,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(e.op(i))) { std::string name = ex_to(e.op(i)).get_name(); if (name == "H") { @@ -2559,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()); } } @@ -2576,7 +2683,7 @@ ex trafo_H_prepend_one(const ex& e, const ex& arg) if (name == "H") { h = e; } else { - for (int i=0; i(e.op(i))) { std::string name = ex_to(e.op(i)).get_name(); if (name == "H") { @@ -2590,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(); } } @@ -2607,7 +2714,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.op(i))) { std::string name = ex_to(e.op(i)).get_name(); if (name == "H") { @@ -2622,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(); } } @@ -2640,7 +2747,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.op(i))) { std::string name = ex_to(e.op(i)).get_name(); if (name == "H") { @@ -2654,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(); } } @@ -2671,7 +2778,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.op(i))) { std::string name = ex_to(e.op(i)).get_name(); if (name == "H") { @@ -2685,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(); } } @@ -2693,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(e) || is_a(e)) { return e.map(*this); @@ -2709,7 +2816,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(buffer)) { - for (int i=0; i 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(e) || is_a(e)) { return e.map(*this); @@ -2810,7 +2917,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(buffer)) { - for (int i=0; i(buffer)) { - for (int i=0; i (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(e) || is_a(e)) { return e.map(*this); @@ -2931,7 +3038,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(buffer)) { - for (int i=0; i(buffer)) { - for (int i=0; i 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 { @@ -3185,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) { @@ -3193,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(); @@ -3209,7 +3316,7 @@ static ex H_eval(const ex& m_, const ex& x) if (is_a(m_)) { m = ex_to(m_); } else { - m = lst(m_); + m = lst{m_}; } if (m.nops() == 0) { return _ex1; @@ -3320,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)); } @@ -3336,7 +3442,7 @@ static ex H_deriv(const ex& m_, const ex& x, unsigned deriv_param) if (is_a(m_)) { m = ex_to(m_); } else { - m = lst(m_); + m = lst{m_}; } ex mb = *m.begin(); if (mb > _ex1) { @@ -3364,9 +3470,9 @@ static void H_print_latex(const ex& m_, const ex& x, const print_context& c) if (is_a(m_)) { m = ex_to(m_); } else { - m = lst(m_); + m = lst{m_}; } - c.s << "\\mbox{H}_{"; + c.s << "\\mathrm{H}_{"; lst::const_iterator itm = m.begin(); (*itm).print(c); itm++; @@ -3397,7 +3503,7 @@ ex convert_H_to_Li(const ex& m, const ex& x) if (is_a(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())); } } @@ -3441,8 +3547,8 @@ static void initcX(std::vector& crX, int Sm = 0; int Smp1 = 0; - std::vector > crG(s.size() - 1, std::vector(L2 + 1)); - for (int m=0; m < s.size() - 1; m++) { + std::vector> crG(s.size() - 1, std::vector(L2 + 1)); + for (int m=0; m < (int)s.size() - 1; m++) { Sm += s[m]; Smp1 = Sm + s[m+1]; for (int i = 0; i <= L2; i++) @@ -3481,12 +3587,12 @@ static cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk, // [Cra] section 4 -static void calc_f(std::vector >& f_kj, +static void calc_f(std::vector>& f_kj, const int maxr, const int L1) { cln::cl_N t0, t1, t2, t3, t4; int i, j, k; - std::vector >::iterator it = f_kj.begin(); + std::vector>::iterator it = f_kj.begin(); cln::cl_F one = cln::cl_float(1, cln::float_format(Digits)); t0 = cln::exp(-lambda); @@ -3510,7 +3616,7 @@ static void calc_f(std::vector >& f_kj, // [Cra] (3.1) static cln::cl_N crandall_Z(const std::vector& s, - const std::vector >& f_kj) + const std::vector>& f_kj) { const int j = s.size(); @@ -3587,7 +3693,7 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector& s) } } - std::vector > f_kj(L1); + std::vector> f_kj(L1); calc_f(f_kj, maxr, L1); const cln::cl_N r0factorial = cln::factorial(r[0]-1); @@ -3668,7 +3774,7 @@ cln::cl_N zeta_do_Hoelder_convolution(const std::vector& 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(m_)) { m = ex_to(m_); } else { - m = lst(m_); + m = lst{m_}; } lst s; if (is_a(s_)) { s = ex_to(s_); } else { - s = lst(s_); + s = lst{s_}; } c.s << "\\zeta("; lst::const_iterator itm = m.begin();