/** @file inifcns_nstdsums.cpp
*
* Implementation of some special functions that have a representation as nested sums.
+ *
* The functions are:
* classical polylogarithm Li(n,x)
* multiple polylogarithm Li(lst(n_1,...,n_k),lst(x_1,...,x_k))
* multiple zeta value zeta(n) or zeta(lst(n_1,...,n_k))
*
* Some remarks:
+ *
* - All formulae used can be looked up in the following publications:
* [Kol] Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
* [Cra] Fast Evaluation of Multiple Zeta Sums, R.E.Crandall, Math.Comp. 67 (1998), pp. 1163-1172.
- * - Classical polylogarithms (Li) and nielsen's generalized polylogarithms (S) can be numerically
- * evaluated in the whole complex plane.
- * - The calculation of classical polylogarithms is speed up by using Euler-Maclaurin summation (EuMac).
- * - The remaining functions can only be numerically evaluated with arguments lying in the unit sphere
- * at the moment.
- * - The functions have no series expansion. To do it, you have to convert these functions
+ * [ReV] Harmonic Polylogarithms, E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
+ *
+ * - The order of parameters and arguments of H, Li and zeta is defined according to their order in the
+ * nested sums representation.
+ *
+ * - Except for the multiple polylogarithm all functions can be nummerically evaluated with arguments in
+ * the whole complex plane. Multiple polylogarithms evaluate only if each argument x_i is smaller than
+ * one. The parameters for every function (n, p or n_i) must be positive integers.
+ *
+ * - The calculation of classical polylogarithms is speed up by using Bernoulli numbers and
+ * look-up tables. S uses look-up tables as well. The zeta function applies the algorithm in
+ * [Cra] for speed up.
+ *
+ * - The functions have no series expansion as nested sums. To do it, you have to convert these functions
* into the appropriate objects from the nestedsums library, do the expansion and convert the
* result back.
+ *
* - Numerical testing of this implementation has been performed by doing a comparison of results
- * between this software and the commercial M.......... 4.1.
+ * between this software and the commercial M.......... 4.1. Multiple zeta values have been checked
+ * by means of evaluations into simple zeta values. Harmonic polylogarithms have been checked by
+ * comparison to S(n,p,x) for corresponding parameter combinations and by continuity checks
+ * around |x|=1 along with comparisons to corresponding zeta functions.
*
*/
#include "utils.h"
#include "wildcard.h"
-//DEBUG
-#include <iostream>
-using namespace std;
namespace GiNaC {
//////////////////////////////////////////////////////////////////////
+// anonymous namespace for helper functions
namespace {
-// lookup table for Euler-Maclaurin optimization
+
+
+// lookup table for factors built from Bernoulli numbers
// see fill_Xn()
std::vector<std::vector<cln::cl_N> > Xn;
int xnsize = 0;
-}
-// This function calculates the X_n. The X_n are needed for the Euler-Maclaurin summation (EMS) of
-// classical polylogarithms.
-// With EMS the polylogs can be calculated as follows:
+// This function calculates the X_n. The X_n are needed for speed up of classical polylogarithms.
+// With these numbers the polylogs can be calculated as follows:
// Li_p (x) = \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with u = -log(1-x)
// X_0(n) = B_n (Bernoulli numbers)
// X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
// X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
// This results in a slightly more complicated algorithm for the X_n.
// The first index in Xn corresponds to the index of the polylog minus 2.
-// The second index in Xn corresponds to the index from the EMS.
-static void fill_Xn(int n)
+// The second index in Xn corresponds to the index from the actual sum.
+void fill_Xn(int n)
{
// rule of thumb. needs to be improved. TODO
const int initsize = Digits * 3 / 2;
}
-// calculates Li(2,x) without EuMac
-static cln::cl_N Li2_do_sum(const cln::cl_N& x)
+// calculates Li(2,x) without Xn
+cln::cl_N Li2_do_sum(const cln::cl_N& x)
{
cln::cl_N res = x;
cln::cl_N resbuf;
}
-// calculates Li(2,x) with EuMac
-static cln::cl_N Li2_do_sum_EuMac(const cln::cl_N& x)
+// calculates Li(2,x) with Xn
+cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x)
{
std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
cln::cl_N u = -cln::log(1-x);
}
-// calculates Li(n,x), n>2 without EuMac
-static cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
+// calculates Li(n,x), n>2 without Xn
+cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
{
cln::cl_N factor = x;
cln::cl_N res = x;
}
-// calculates Li(n,x), n>2 with EuMac
-static cln::cl_N Lin_do_sum_EuMac(int n, const cln::cl_N& x)
+// calculates Li(n,x), n>2 with Xn
+cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
{
std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
cln::cl_N u = -cln::log(1-x);
// forward declaration needed by function Li_projection and C below
-static numeric S_num(int n, int p, const numeric& x);
+numeric S_num(int n, int p, const numeric& x);
// helper function for classical polylog Li
-static cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
+cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
{
// treat n=2 as special case
if (n == 2) {
return Li2_do_sum(x);
} else {
- return Li2_do_sum_EuMac(x);
+ return Li2_do_sum_Xn(x);
}
} else {
// choose the faster algorithm
if (cln::abs(cln::realpart(x)) > 0.75) {
return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
} else {
- return -Li2_do_sum_EuMac(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
+ return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
}
}
} else {
if (cln::realpart(x) < 0.5) {
// choose the faster algorithm
- // with n>=12 the "normal" summation always wins against EuMac
+ // with n>=12 the "normal" summation always wins against the method with Xn
if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
return Lin_do_sum(n, x);
} else {
- return Lin_do_sum_EuMac(n, x);
+ 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);
// helper function for classical polylog Li
-static numeric Li_num(int n, const numeric& x)
+numeric Li_num(int n, const numeric& x)
{
if (n == 1) {
// just a log
}
+} // end of anonymous namespace
+
+
//////////////////////////////////////////////////////////////////////
//
// Multiple polylogarithm Li
//////////////////////////////////////////////////////////////////////
-static cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
+// anonymous namespace for helper function
+namespace {
+
+
+cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
{
const int j = s.size();
}
+} // end of anonymous namespace
+
+
//////////////////////////////////////////////////////////////////////
//
// Classical polylogarithm and multiple polylogarithm Li
std::vector<int> m;
std::vector<cln::cl_N> x;
- for (int i=ex_to<numeric>(x1.nops()).to_int()-1; i>=0; i--) {
+ for (int i=0; i<ex_to<numeric>(x1.nops()).to_int(); i++) {
m.push_back(ex_to<numeric>(x1.op(i)).to_int());
x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
}
//////////////////////////////////////////////////////////////////////
+// anonymous namespace for helper functions
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;
int ynsize = 0; // number of Yn[]
int ynlength = 100; // initial length of all Yn[i]
-}
// This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
// The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
// representing S_{n,p}(x).
// The calculation of Y_n uses the values from Y_{n-1}.
-static void fill_Yn(int n, const cln::float_format_t& prec)
+void fill_Yn(int n, const cln::float_format_t& prec)
{
const int initsize = ynlength;
//const int initsize = initsize_Yn;
// make Yn longer ...
-static void make_Yn_longer(int newsize, const cln::float_format_t& prec)
+void make_Yn_longer(int newsize, const cln::float_format_t& prec)
{
cln::cl_N one = cln::cl_float(1, prec);
// helper function for S(n,p,x)
// [Kol] (7.2)
-static cln::cl_N C(int n, int p)
+cln::cl_N C(int n, int p)
{
cln::cl_N result;
// helper function for S(n,p,x)
// [Kol] remark to (9.1)
-static cln::cl_N a_k(int k)
+cln::cl_N a_k(int k)
{
cln::cl_N result;
// helper function for S(n,p,x)
// [Kol] remark to (9.1)
-static cln::cl_N b_k(int k)
+cln::cl_N b_k(int k)
{
cln::cl_N result;
// helper function for S(n,p,x)
-static cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
+cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
{
if (p==1) {
return Li_projection(n+1, x, prec);
// helper function for S(n,p,x)
-static cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
+cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
{
// [Kol] (5.3)
if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
// helper function for S(n,p,x)
-static numeric S_num(int n, int p, const numeric& x)
+numeric S_num(int n, int p, const numeric& x)
{
if (x == 1) {
if (n == 1) {
}
+} // end of anonymous namespace
+
+
//////////////////////////////////////////////////////////////////////
//
// Nielsen's generalized polylogarithm S
//////////////////////////////////////////////////////////////////////
-static ex convert_from_RV(const lst& parameterlst, const ex& arg);
+// anonymous namespace for helper functions
+namespace {
+
+// forward declaration
+ex convert_from_RV(const lst& parameterlst, const ex& arg);
-static ex trafo_H_mult(const ex& h1, const ex& h2)
+
+// multiplies an one-dimensional H with another H
+// [ReV] (18)
+ex trafo_H_mult(const ex& h1, const ex& h2)
{
ex res;
ex hshort;
}
-namespace {
+// applies trafo_H_mult recursively on expressions
struct map_trafo_H_mult : public map_function
{
ex operator()(const ex& e)
return e;
}
};
-} // anonymous namespace
-static ex trafo_H_prepend_one(const ex& e, const ex& arg)
+// do integration [ReV] (49)
+// put parameter 1 in front of existing parameters
+ex trafo_H_prepend_one(const ex& e, const ex& arg)
{
ex h;
std::string name;
}
-static ex trafo_H_prepend_zero(const ex& e, const ex& arg)
+// do integration [ReV] (55)
+// put parameter 0 in front of existing parameters
+ex trafo_H_prepend_zero(const ex& e, const ex& arg)
{
ex h;
std::string name;
}
-namespace {
+// do x -> 1-x transformation
struct map_trafo_H_1mx : public map_function
{
ex operator()(const ex& e)
};
+// do x -> 1/x transformation
struct map_trafo_H_1overx : public map_function
{
ex operator()(const ex& e)
};
+// remove trailing zeros from H-parameters
struct map_trafo_H_reduce_trailing_zeros : public map_function
{
ex operator()(const ex& e)
};
+// recursively call convert_from_RV on expression
struct map_trafo_H_convert : public map_function
{
ex operator()(const ex& e)
return e;
}
};
-} // anonymous namespace
-static lst convert_to_RV(const lst& o)
+// translate notation from nested sums to Remiddi/Vermaseren
+lst convert_to_RV(const lst& o)
{
lst res;
for (lst::const_iterator it = o.begin(); it != o.end(); it++) {
}
-static ex convert_from_RV(const lst& parameterlst, const ex& arg)
+// translate notation from Remiddi/Vermaseren to nested sums
+ex convert_from_RV(const lst& parameterlst, const ex& arg)
{
lst newparameterlst;
}
-static cln::cl_N H_do_sum(const std::vector<int>& s, const cln::cl_N& x)
+// do the actual summation.
+cln::cl_N H_do_sum(const std::vector<int>& s, const cln::cl_N& x)
{
const int j = s.size();
}
+} // end of anonymous namespace
+
+
//////////////////////////////////////////////////////////////////////
//
// Harmonic polylogarithm H
return res.subs(xtemp==x2).evalf();
}
+ // since the x->1-x transformation produces a lot of terms, it is only
+ // efficient for argument near one.
if (cln::realpart(x) > 0.95) {
symbol xtemp("xtemp");
map_trafo_H_1mx trafo;
return res.subs(xtemp==x2).evalf();
}
- // no trafo
+ // no trafo -> do summation
int count = x1.nops();
std::vector<int> r(count);
for (int i=0; i<count; i++) {
//////////////////////////////////////////////////////////////////////
+// anonymous namespace for helper functions
namespace {
+
+
+// parameters and data for [Cra] algorithm
const cln::cl_N lambda = cln::cl_N("319/320");
int L1;
int L2;
std::vector<cln::cl_N> crB;
std::vector<std::vector<cln::cl_N> > crG;
std::vector<cln::cl_N> crX;
-}
-static void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
+void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
{
const int size = a.size();
for (int n=0; n<size; n++) {
// [Cra] section 4
-static void initcX(const std::vector<int>& s)
+void initcX(const std::vector<int>& s)
{
const int k = s.size();
// [Cra] section 4
-static cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
+cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
{
cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
cln::cl_N factor = cln::expt(lambda, Sqk);
// [Cra] section 4
-static void calc_f(int maxr)
+void calc_f(int maxr)
{
f_kj.clear();
f_kj.resize(L1);
// [Cra] (3.1)
-static cln::cl_N crandall_Z(const std::vector<int>& s)
+cln::cl_N crandall_Z(const std::vector<int>& s)
{
const int j = s.size();
// [Cra] (2.4)
-static cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
+cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
{
std::vector<int> r = s;
const int j = r.size();
}
-static cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
+cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
{
const int j = r.size();
}
+} // end of anonymous namespace
+
+
//////////////////////////////////////////////////////////////////////
//
// Multiple zeta values zeta
git://www.ginac.de / ginac.git / commitdiff