X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=17db2c885dfd63a8ea880d61ff85ec4be5e8ea0f;hp=3b7fc556892398c64180d6ffc78f4cdf4d242593;hb=e5e2ac399c225a8c7248e1879e192eba00067538;hpb=0a2231e31d5d0eb780ebb7ffb8d0628c78b74248 diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index 3b7fc556..17db2c88 100644 --- a/ginac/inifcns_nstdsums.cpp +++ b/ginac/inifcns_nstdsums.cpp @@ -29,7 +29,7 @@ * 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. * - * - The calculation of classical polylogarithms is speed up by using Bernoulli numbers and + * - The calculation of classical polylogarithms is speeded up by using Bernoulli numbers and * look-up tables. S uses look-up tables as well. The zeta function applies the algorithms in * [Cra] and [BBB] for speed up. * @@ -102,6 +102,9 @@ namespace { // lookup table for factors built from Bernoulli numbers // see fill_Xn() std::vector > Xn; +// initial size of Xn that should suffice for 32bit machines (must be even) +const int xninitsizestep = 26; +int xninitsize = xninitsizestep; int xnsize = 0; @@ -117,17 +120,14 @@ int xnsize = 0; // 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; - if (n>1) { // calculate X_2 and higher (corresponding to Li_4 and higher) - std::vector buf(initsize); + std::vector buf(xninitsize); 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++) { + for (int i=2; i<=xninitsize; i++) { if (i&1) { result = 0; // k == 0 } else { @@ -147,14 +147,14 @@ void fill_Xn(int n) Xn.push_back(buf); } else if (n==1) { // special case to handle the X_0 correct - std::vector buf(initsize); + std::vector buf(xninitsize); 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++) { + for (int i=3; i<=xninitsize; 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; @@ -171,9 +171,9 @@ void fill_Xn(int n) Xn.push_back(buf); } else { // calculate X_0 - std::vector buf(initsize/2); + std::vector buf(xninitsize/2); std::vector::iterator it = buf.begin(); - for (int i=1; i<=initsize/2; i++) { + for (int i=1; i<=xninitsize/2; i++) { *it = bernoulli(i*2).to_cl_N(); it++; } @@ -184,6 +184,53 @@ void fill_Xn(int n) } +// doubles the number of entries in each Xn[] +void double_Xn() +{ + const int pos0 = xninitsize / 2; + // X_0 + for (int i=1; i<=xninitsizestep/2; ++i) { + Xn[0].push_back(bernoulli((i+pos0)*2).to_cl_N()); + } + if (Xn.size() > 0) { + int xend = xninitsize + xninitsizestep; + cln::cl_N result; + // X_1 + for (int i=xninitsize+1; i<=xend; ++i) { + if (i & 1) { + result = -Xn[0][(i-3)/2]/2; + Xn[1].push_back((cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result); + } else { + result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1); + 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 + Xn[n].push_back(result); + } + } + } + xninitsize += xninitsizestep; +} + + // calculates Li(2,x) without Xn cln::cl_N Li2_do_sum(const cln::cl_N& x) { @@ -207,6 +254,7 @@ cln::cl_N Li2_do_sum(const cln::cl_N& x) cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x) { std::vector::const_iterator it = Xn[0].begin(); + std::vector::const_iterator xend = Xn[0].end(); cln::cl_N u = -cln::log(1-x); cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits)); cln::cl_N res = u - u*u/4; @@ -216,8 +264,12 @@ cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x) resbuf = res; factor = factor * u*u / (2*i * (2*i+1)); res = res + (*it) * factor; - it++; // should we check it? or rely on initsize? ... i++; + if (++it == xend) { + double_Xn(); + it = Xn[0].begin() + (i-1); + xend = Xn[0].end(); + } } while (res != resbuf); return res; } @@ -244,6 +296,7 @@ cln::cl_N Lin_do_sum(int n, const cln::cl_N& x) cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x) { std::vector::const_iterator it = Xn[n-2].begin(); + std::vector::const_iterator xend = Xn[n-2].end(); cln::cl_N u = -cln::log(1-x); cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits)); cln::cl_N res = u; @@ -253,8 +306,12 @@ cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x) resbuf = res; factor = factor * u / i; res = res + (*it) * factor; - it++; // should we check it? or rely on initsize? ... i++; + if (++it == xend) { + double_Xn(); + it = Xn[n-2].begin() + (i-2); + xend = Xn[n-2].end(); + } } while (res != resbuf); return res; } @@ -939,9 +996,8 @@ numeric S_num(int n, int p, const numeric& x) else if (!x.imag().is_rational()) prec = cln::float_format(cln::the(cln::imagpart(value))); - // [Kol] (5.3) - if (cln::realpart(value) < -0.5) { + if ((cln::realpart(value) < -0.5) || (n == 0)) { 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); @@ -1029,6 +1085,10 @@ static ex S_eval(const ex& n, const ex& p, const ex& x) return S_num(ex_to(n).to_int(), ex_to(p).to_int(), ex_to(x)); } } + if (n.is_zero()) { + // [Kol] (5.3) + return pow(-log(1-x), p) / factorial(p); + } return S(n, p, x).hold(); } @@ -2632,7 +2692,7 @@ static void zeta1_print_latex(const ex& m_, const print_context& c) } -unsigned zeta1_SERIAL::serial = function::register_new(function_options("zeta"). +unsigned zeta1_SERIAL::serial = function::register_new(function_options("zeta", 1). evalf_func(zeta1_evalf). eval_func(zeta1_eval). derivative_func(zeta1_deriv). @@ -2769,7 +2829,7 @@ static void zeta2_print_latex(const ex& m_, const ex& s_, const print_context& c } -unsigned zeta2_SERIAL::serial = function::register_new(function_options("zeta"). +unsigned zeta2_SERIAL::serial = function::register_new(function_options("zeta", 2). evalf_func(zeta2_evalf). eval_func(zeta2_eval). derivative_func(zeta2_deriv).