* - 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.
namespace GiNaC {
+// lookup table for Euler-MacLaurin optimization
+// see fill_Xn()
+std::vector<std::vector<cln::cl_N> > Xn;
+int xnsize = 0;
+
+
+// lookup table for Euler-Zagier-Sums (used for S_n,p(x))
+// see fill_Yn()
+std::vector<std::vector<cln::cl_N> > 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<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)
+{
+ // TODO -> get rid of the magic number
+ const int initsize = initsize_Yn;
+
+ 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);
+ 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<cln::cl_N> buf(initsize);
+ std::vector<cln::cl_N>::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<n-1; i++) {
+ fill_Xn(i);
+ }
+ }
+
+ // using Euler-MacLaurin summation
+ if (n==2) {
+ // Li_2. X_0 is special ...
+ std::vector<cln::cl_N>::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<cln::cl_N>::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<n-1; j++) {
+ result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
+ * cln::expt(cln::log(x), j) / cln::factorial(j) ;
+ }
+ return result;
+ }
+ }
}
}
-// forward declaration needed by function C below
-static numeric S_num(int n, int p, const numeric& x);
-
-
// helper function for S(n,p,x)
// [Kol] (7.2)
static cln::cl_N C(int n, int p)
// helper function for S(n,p,x)
static cln::cl_N S_series(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
{
- n++;
- int i = p;
- p--;
- cln::cl_N aug, acc;
- cln::cl_N num = cln::expt(x,p);
- cln::cl_N converter = cln::complex(cln::cl_float(1, prec), 0);
- cln::cl_N den = 0;
- do {
- num = num * x;
- den = cln::expt(cln::cl_I(i), n);
- aug = num / den * numeric_nielsen(i, p);
- i++;
- acc = acc + aug;
- } while (acc != acc+aug);
-
- return acc;
+ if (p==1) {
+ return Li_series(n+1, x, prec);
+ }
+
+ // TODO -> check for vector boundaries and do missing calculations
+
+ // check if precalculated values are sufficient
+ if (p > ynsize+1) {
+ for (int i=ynsize; i<p-1; i++) {
+ fill_Yn(i);
+ }
+ }
+
+ cln::cl_N result;
+ cln::cl_N resultbuffer;
+ for (int i=p; true; i++) {
+ resultbuffer = result;
+ result = result + cln::expt(x,i) / cln::expt(cln::cl_I(i),n+1) * Yn[p-2][i-p]; // should we check it? or rely on magic number? ...
+ if (cln::zerop(result-resultbuffer)) {
+ break;
+ }
+ }
+
+ return result;
}
git://www.ginac.de / ginac.git / commitdiff