+// 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:
+// 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)
+// The calculation of Xn depends on X0 and X{n-1}.
+// 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)
+{
+ // rule of thumb. needs to be improved. TODO
+ const int initsize = Digits * 3 / 2;
+
+ if (n>1) {
+ // calculate X_2 and higher (corresponding to Li_4 and higher)
+ std::vector<cln::cl_N> buf(initsize);
+ std::vector<cln::cl_N>::iterator it = buf.begin();
+ cln::cl_N result;
+ *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
+ it++;
+ for (int i=2; i<=initsize; i++) {
+ if (i&1) {
+ result = 0; // k == 0
+ } else {
+ result = Xn[0][i/2-1]; // k == 0
+ }
+ for (int k=1; k<i-1; k++) {
+ if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
+ result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
+ }
+ }
+ result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
+ result = result + Xn[n-1][i-1] / (i+1); // k == i
+
+ *it = result;
+ it++;
+ }
+ Xn.push_back(buf);
+ } else if (n==1) {
+ // special case to handle the X_0 correct
+ std::vector<cln::cl_N> buf(initsize);
+ std::vector<cln::cl_N>::iterator it = buf.begin();
+ cln::cl_N result;
+ *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
+ it++;
+ *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
+ it++;
+ for (int i=3; i<=initsize; i++) {
+ if (i & 1) {
+ result = -Xn[0][(i-3)/2]/2;
+ *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
+ it++;
+ } else {
+ result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
+ for (int k=1; k<i/2; k++) {
+ result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
+ }
+ *it = result;
+ it++;
+ }
+ }
+ Xn.push_back(buf);
+ } else {
+ // calculate X_0
+ std::vector<cln::cl_N> buf(initsize/2);
+ std::vector<cln::cl_N>::iterator it = buf.begin();
+ for (int i=1; i<=initsize/2; i++) {
+ *it = bernoulli(i*2).to_cl_N();
+ it++;
+ }
+ Xn.push_back(buf);
+ }
+
+ xnsize++;
+}
+
+
+// This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
+// The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
+// representing S_{n,p}(x).
+// The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
+// equivalent Z-sum.
+// 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)
+{
+ // TODO -> get rid of the magic number
+ const int initsize = ynlength;
+ //const int initsize = initsize_Yn;
+ cln::cl_N one = cln::cl_float(1, prec);
+
+ if (n) {
+ std::vector<cln::cl_N> buf(initsize);
+ std::vector<cln::cl_N>::iterator it = buf.begin();
+ std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
+ *it = (*itprev) / cln::cl_N(n+1) * one;
+ it++;
+ itprev++;
+ // sums with an index smaller than the depth are zero and need not to be calculated.
+ // calculation starts with depth, which is n+2)
+ for (int i=n+2; i<=initsize+n; i++) {
+ *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
+ it++;
+ itprev++;
+ }
+ Yn.push_back(buf);
+ } else {
+ std::vector<cln::cl_N> buf(initsize);
+ std::vector<cln::cl_N>::iterator it = buf.begin();
+ *it = 1 * one;
+ it++;
+ for (int i=2; i<=initsize; i++) {
+ *it = *(it-1) + 1 / cln::cl_N(i) * one;
+ it++;
+ }
+ Yn.push_back(buf);
+ }
+ ynsize++;
+}
+
+
+// make Yn longer ...
+static void make_Yn_longer(int newsize, const cln::float_format_t& prec)
+{
+
+ cln::cl_N one = cln::cl_float(1, prec);
+
+ Yn[0].resize(newsize);
+ std::vector<cln::cl_N>::iterator it = Yn[0].begin();
+ it += ynlength;
+ for (int i=ynlength+1; i<=newsize; i++) {
+ *it = *(it-1) + 1 / cln::cl_N(i) * one;
+ it++;
+ }
+
+ for (int n=1; n<ynsize; n++) {
+ Yn[n].resize(newsize);
+ std::vector<cln::cl_N>::iterator it = Yn[n].begin();
+ std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
+ it += ynlength;
+ itprev += ynlength;
+ for (int i=ynlength+n+1; i<=newsize+n; i++) {
+ *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
+ it++;
+ itprev++;
+ }
+ }
+
+ ynlength = newsize;
+}
+
+
+// calculates Li(2,x) without EuMac
+static cln::cl_N Li2_series(const cln::cl_N& x)
git://www.ginac.de / ginac.git / commitdiff