From 9dfecd4a3e529be5e8310ec08b1a8e6687860f6e Mon Sep 17 00:00:00 2001 From: Jens Vollinga Date: Mon, 22 Sep 2003 21:26:04 +0000 Subject: [PATCH] Classical polylog now uses Euler-MacLaurin summation. Some speed improvements for S. --- ginac/inifcns_nstdsums.cpp | 263 ++++++++++++++++++++++++++++++++----- 1 file changed, 227 insertions(+), 36 deletions(-) diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index f0cd6ab9..c094c50b 100644 --- a/ginac/inifcns_nstdsums.cpp +++ b/ginac/inifcns_nstdsums.cpp @@ -15,6 +15,7 @@ * - Classical polylogarithms (Li) and nielsen's generalized polylogarithms (S) can be numerically * evaluated in the whole complex plane except for S(n,p,-1) when p is not unit (no formula yet * to tackle these points). And of course, there is still room for speed optimizations ;-). + * - The calculation of classical polylogarithms is speed up by using Euler-MacLaurin summation. * - The remaining functions can only be numerically evaluated with arguments lying in the unit sphere * at the moment. Sorry. The evaluation especially for mZeta is very slow ... better not use it * right now. @@ -59,35 +60,221 @@ namespace GiNaC { +// lookup table for Euler-MacLaurin optimization +// see fill_Xn() +std::vector > Xn; +int xnsize = 0; + + +// lookup table for Euler-Zagier-Sums (used for S_n,p(x)) +// see fill_Yn() +std::vector > Yn; +int ynsize = 0; +//TODO EVIL MAGIC NUMBER !!! but first the transformations for S have to improve ... +const int initsize_Yn = 2000; + + ////////////////////// // helper functions // ////////////////////// -// helper function for classical polylog Li +// 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 buf(initsize); + std::vector::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 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 buf(initsize); + std::vector::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 buf(initsize/2); + std::vector::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) +{ + // TODO -> get rid of the magic number + const int initsize = initsize_Yn; + + if (n) { + std::vector buf(initsize); + std::vector::iterator it = buf.begin(); + std::vector::iterator itprev = Yn[n-1].begin(); + *it = (*itprev) / cln::cl_N(n+1); + 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); + it++; + itprev++; + } + Yn.push_back(buf); + } else { + std::vector buf(initsize); + std::vector::iterator it = buf.begin(); + *it = 1; + it++; + for (int i=2; i<=initsize; i++) { + *it = *(it-1) + 1 / cln::cl_N(i); + it++; + } + Yn.push_back(buf); + } + ynsize++; +} + + static cln::cl_N Li_series(int n, const cln::cl_N& x, const cln::float_format_t& prec) { - // Note: argument must be in the unit circle - cln::cl_N aug, acc; - cln::cl_N num = cln::complex(cln::cl_float(1, prec), 0); - cln::cl_N den = 0; - int i = 1; - do { - num = num * x; - cln::cl_R ii = i; - den = cln::expt(ii, n); - i++; - aug = num / den; - acc = acc + aug; - } while (acc != acc+aug); - return acc; + // check if precalculated values are sufficient + if (n > xnsize+1) { + for (int i=xnsize; i::const_iterator it = Xn[0].begin(); + cln::cl_N u = -cln::log(cln::complex(cln::cl_float(1, prec), 0)-x); + cln::cl_N factor = u; + cln::cl_N res = u - u*u/4; + cln::cl_N resbuf; + for (int i=1; true; i++) { + resbuf = res; + factor = factor * u*u / (2*i * (2*i+1)); + res = res + (*it) * factor; + it++; // should we check it? or rely on initsize? ... + if (cln::zerop(res-resbuf)) + { + break; + } + } + return res; + } else { + // Li_3 and higher + std::vector::const_iterator it = Xn[n-2].begin(); + cln::cl_N u = -cln::log(cln::complex(cln::cl_float(1, prec), 0)-x); + cln::cl_N factor = u; + cln::cl_N res = u; + cln::cl_N resbuf; + for (int i=1; true; i++) { + resbuf = res; + factor = factor * u / (i+1); + res = res + (*it) * factor; + it++; // should we check it? or rely on initsize? ... + if (cln::zerop(res-resbuf)) + { + // should not be needed. +// if (!cln::zerop(*it)) { + break; +// } + } + } + return res; + } } +// forward declaration needed by function C below +static 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) { - return Li_series(n, x, prec); + if (cln::realpart(x) < 0.5) { + return Li_series(n, x, prec); + } else { + if (n==2) { + return -Li_series(2, 1-x, prec) - cln::log(x) * cln::log(1-x) + cln::zeta(2); + } else { + cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1); + for (int j=0; j check for vector boundaries and do missing calculations + + // check if precalculated values are sufficient + if (p > ynsize+1) { + for (int i=ynsize; i