* 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.
*
*/
/*
- * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
// lookup table for factors built from Bernoulli numbers
// see fill_Xn()
std::vector<std::vector<cln::cl_N> > Xn;
+// initial size of Xn that should suffice for 32bit machines (must be even)
+const int xninitsizestep = 26;
+int xninitsize = xninitsizestep;
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<cln::cl_N> buf(initsize);
+ std::vector<cln::cl_N> buf(xninitsize);
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++) {
+ for (int i=2; i<=xninitsize; i++) {
if (i&1) {
result = 0; // k == 0
} else {
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> buf(xninitsize);
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++) {
+ 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;
Xn.push_back(buf);
} else {
// calculate X_0
- std::vector<cln::cl_N> buf(initsize/2);
+ std::vector<cln::cl_N> buf(xninitsize/2);
std::vector<cln::cl_N>::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++;
}
}
+// 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<i/2; k++) {
+ result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
+ }
+ Xn[1].push_back(result);
+ }
+ }
+ // X_n
+ for (int n=2; n<Xn.size(); ++n) {
+ for (int i=xninitsize+1; i<=xend; ++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
+ Xn[n].push_back(result);
+ }
+ }
+ }
+ xninitsize += xninitsizestep;
+}
+
+
// calculates Li(2,x) without Xn
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<cln::cl_N>::const_iterator it = Xn[0].begin();
+ std::vector<cln::cl_N>::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;
+ cln::cl_N uu = cln::square(u);
+ cln::cl_N res = u - uu/4;
cln::cl_N resbuf;
unsigned i = 1;
do {
resbuf = res;
- factor = factor * u*u / (2*i * (2*i+1));
+ factor = factor * uu / (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;
}
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();
+ std::vector<cln::cl_N>::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;
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;
}
if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
}
+ if (is_a<numeric>(x1) && !is_a<lst>(x2)) {
+ // try to numerically evaluate second argument
+ ex x2_val = x2.evalf();
+ if (is_a<numeric>(x2_val)) {
+ return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2_val));
+ } else {
+ return Li(x1, x2).hold();
+ }
+ }
// multiple polylogs
else if (is_a<lst>(x1) && is_a<lst>(x2)) {
ex conv = 1;
else if (!x.imag().is_rational())
prec = cln::float_format(cln::the<cln::cl_F>(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);
static ex S_evalf(const ex& n, const ex& p, const ex& x)
{
- if (n.info(info_flags::posint) && p.info(info_flags::posint) && is_a<numeric>(x)) {
- return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
+ if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
+ if (is_a<numeric>(x)) {
+ return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
+ } else {
+ ex x_val = x.evalf();
+ if (is_a<numeric>(x_val)) {
+ return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x_val));
+ }
+ }
}
return S(n, p, x).hold();
}
return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
}
}
+ if (n.is_zero()) {
+ // [Kol] (5.3)
+ return pow(-log(1-x), p) / factorial(p);
+ }
return S(n, p, x).hold();
}
// anonymous namespace for helper functions
namespace {
+
+// regulates the pole (used by 1/x-transformation)
+symbol H_polesign("IMSIGN");
+
// convert parameters from H to Li representation
// parameters are expected to be in expanded form, i.e. only 0, 1 and -1
}
if (allthesame) {
map_trafo_H_mult unify;
- return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() - I*Pi, parameter.nops())
+ return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() + H_polesign, parameter.nops())
/ factorial(parameter.nops())).expand());
}
}
static ex H_evalf(const ex& x1, const ex& x2)
{
- if (is_a<lst>(x1) && is_a<numeric>(x2)) {
+ if (is_a<lst>(x1)) {
+
+ cln::cl_N x;
+ if (is_a<numeric>(x2)) {
+ x = ex_to<numeric>(x2).to_cl_N();
+ } else {
+ ex x2_val = x2.evalf();
+ if (is_a<numeric>(x2_val)) {
+ x = ex_to<numeric>(x2_val).to_cl_N();
+ }
+ }
+
for (int i=0; i<x1.nops(); i++) {
if (!x1.op(i).info(info_flags::integer)) {
- return H(x1,x2).hold();
+ return H(x1, x2).hold();
}
}
if (x1.nops() < 1) {
- return H(x1,x2).hold();
+ return H(x1, x2).hold();
}
- cln::cl_N x = ex_to<numeric>(x2).to_cl_N();
-
const lst& morg = ex_to<lst>(x1);
// remove trailing zeros ...
if (*(--morg.end()) == 0) {
// x -> 1/x
map_trafo_H_1overx trafo;
res *= trafo(H(m, xtemp));
+ if (cln::imagpart(x) <= 0) {
+ res = res.subs(H_polesign == -I*Pi);
+ } else {
+ res = res.subs(H_polesign == I*Pi);
+ }
}
// simplify result
}
-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).
}
-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).
git://www.ginac.de / ginac.git / blobdiff