/** @file inifcns_nstdsums.cpp * * Implementation of some special functions that have a representation as nested sums. * The functions are: * classical polylogarithm Li(n,x) * multiple polylogarithm Li(lst(n_1,...,n_k),lst(x_1,...,x_k) * nielsen's generalized polylogarithm S(n,p,x) * harmonic polylogarithm H(lst(m_1,...,m_k),x) * multiple zeta value mZeta(lst(m_1,...,m_k)) * * Some remarks: * - All formulae used can be looked up in the following publication: * Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258. * This document will be referenced as [Kol] throughout this source code. * - 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. * - The functions have no series expansion. To do it, you have to convert these functions * into the appropriate objects from the nestedsums library, do the expansion and convert the * result back. * - Numerical testing of this implementation has been performed by doing a comparison of results * between this software and the commercial M.......... 4.1. * */ /* * GiNaC Copyright (C) 1999-2003 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 * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include #include #include #include "inifcns.h" #include "lst.h" #include "numeric.h" #include "operators.h" #include "relational.h" #include "pseries.h" 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 // ////////////////////// // 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) { // 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) { 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(cln::realpart(value))); else if (!x.imag().is_rational()) prec = cln::float_format(cln::the(cln::imagpart(value))); // [Kol] (5.15) if (cln::abs(value) > 1) { cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n); // check if argument is complex. if it is real, the new polylog has to be conjugated. if (cln::zerop(cln::imagpart(value))) { if (n & 1) { result = result + conjugate(Li_projection(n, cln::recip(value), prec)); } else { result = result - conjugate(Li_projection(n, cln::recip(value), prec)); } } else { if (n & 1) { result = result + Li_projection(n, cln::recip(value), prec); } else { result = result - Li_projection(n, cln::recip(value), prec); } } cln::cl_N add; 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 cln::cl_F("0.5")) { cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n) * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p); for (int s=0; s(cln::realpart(value))); else if (!x.imag().is_rational()) prec = cln::float_format(cln::the(cln::imagpart(value))); // [Kol] (5.3) if (cln::realpart(value) < -0.5) { 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); for (int s=0; s 1) { cln::cl_N result; for (int s=0; s& x, std::vector& m) { cln::cl_N res; if (x.empty()) { return 1; } for (int i=1; i::iterator be; std::vector::iterator en; be = x.begin(); be++; en = x.end(); std::vector xbuf(be, en); be = m.begin(); be++; en = m.end(); std::vector mbuf(be, en); res = res + cln::expt(x[0],i) / cln::expt(i,m[0]) * numeric_zsum(i, xbuf, mbuf); } return res; } // helper function for harmonic polylogarithm static cln::cl_N numeric_harmonic(int n, std::vector& m) { cln::cl_N res; if (m.empty()) { return 1; } for (int i=1; i::iterator be; std::vector::iterator en; be = m.begin(); be++; en = m.end(); std::vector mbuf(be, en); res = res + cln::recip(cln::expt(i,m[0])) * numeric_harmonic(i, mbuf); } return res; } ///////////////////////////// // end of helper functions // ///////////////////////////// // Polylogarithm and multiple polylogarithm static ex Li_eval(const ex& x1, const ex& x2) { if (x2.is_zero()) { return 0; } else { return Li(x1,x2).hold(); } } static ex Li_evalf(const ex& x1, const ex& x2) { // classical polylogs if (is_a(x1) && is_a(x2)) { return Li_num(ex_to(x1).to_int(), ex_to(x2)); } // multiple polylogs else if (is_a(x1) && is_a(x2)) { for (int i=0; i(x1.op(i))) return Li(x1,x2).hold(); if (!is_a(x2.op(i))) return Li(x1,x2).hold(); if (x2.op(i) >= 1) return Li(x1,x2).hold(); } cln::cl_N m_1 = ex_to(x1.op(x1.nops()-1)).to_cl_N(); cln::cl_N x_1 = ex_to(x2.op(x2.nops()-1)).to_cl_N(); std::vector x; std::vector m; const int nops = ex_to(x1.nops()).to_int(); for (int i=nops-2; i>=0; i--) { m.push_back(ex_to(x1.op(i)).to_cl_N()); x.push_back(ex_to(x2.op(i)).to_cl_N()); } cln::cl_N res; cln::cl_N resbuf; for (int i=nops; true; i++) { resbuf = res; res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_zsum(i, x, m); if (cln::zerop(res-resbuf)) break; } return numeric(res); } return Li(x1,x2).hold(); } static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options) { epvector seq; seq.push_back(expair(Li(x1,x2), 0)); return pseries(rel,seq); } REGISTER_FUNCTION(Li, eval_func(Li_eval).evalf_func(Li_evalf).do_not_evalf_params().series_func(Li_series)); // Nielsen's generalized polylogarithm static ex S_eval(const ex& x1, const ex& x2, const ex& x3) { if (x2 == 1) { return Li(x1+1,x3); } return S(x1,x2,x3).hold(); } static ex S_evalf(const ex& x1, const ex& x2, const ex& x3) { if (is_a(x1) && is_a(x2) && is_a(x3)) { if ((x3 == -1) && (x2 != 1)) { // no formula to evaluate this ... sorry return S(x1,x2,x3).hold(); } return S_num(ex_to(x1).to_int(), ex_to(x2).to_int(), ex_to(x3)); } return S(x1,x2,x3).hold(); } static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options) { epvector seq; seq.push_back(expair(S(x1,x2,x3), 0)); return pseries(rel,seq); } REGISTER_FUNCTION(S, eval_func(S_eval).evalf_func(S_evalf).do_not_evalf_params().series_func(S_series)); // Harmonic polylogarithm static ex H_eval(const ex& x1, const ex& x2) { return H(x1,x2).hold(); } static ex H_evalf(const ex& x1, const ex& x2) { if (is_a(x1) && is_a(x2)) { for (int i=0; i(x1.op(i))) return H(x1,x2).hold(); } if (x2 >= 1) { return H(x1,x2).hold(); } cln::cl_N m_1 = ex_to(x1.op(x1.nops()-1)).to_cl_N(); cln::cl_N x_1 = ex_to(x2).to_cl_N(); std::vector m; const int nops = ex_to(x1.nops()).to_int(); for (int i=nops-2; i>=0; i--) { m.push_back(ex_to(x1.op(i)).to_cl_N()); } cln::cl_N res; cln::cl_N resbuf; for (int i=nops; true; i++) { resbuf = res; res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_harmonic(i, m); if (cln::zerop(res-resbuf)) break; } return numeric(res); } return H(x1,x2).hold(); } static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options) { epvector seq; seq.push_back(expair(H(x1,x2), 0)); return pseries(rel,seq); } REGISTER_FUNCTION(H, eval_func(H_eval).evalf_func(H_evalf).do_not_evalf_params().series_func(H_series)); // Multiple zeta value static ex mZeta_eval(const ex& x1) { return mZeta(x1).hold(); } static ex mZeta_evalf(const ex& x1) { if (is_a(x1)) { for (int i=0; i(x1.op(i))) return mZeta(x1).hold(); } cln::cl_N m_1 = ex_to(x1.op(x1.nops()-1)).to_cl_N(); std::vector m; const int nops = ex_to(x1.nops()).to_int(); for (int i=nops-2; i>=0; i--) { m.push_back(ex_to(x1.op(i)).to_cl_N()); } cln::float_format_t prec = cln::default_float_format; cln::cl_N res = cln::complex(cln::cl_float(0, prec), 0); cln::cl_N resbuf; for (int i=nops; true; i++) { // to infinity and beyond ... timewise resbuf = res; res = res + cln::recip(cln::expt(i,m_1)) * numeric_harmonic(i, m); if (cln::zerop(res-resbuf)) break; } return numeric(res); } return mZeta(x1).hold(); } static ex mZeta_series(const ex& x1, const relational& rel, int order, unsigned options) { epvector seq; seq.push_back(expair(mZeta(x1), 0)); return pseries(rel,seq); } REGISTER_FUNCTION(mZeta, eval_func(mZeta_eval).evalf_func(mZeta_evalf).do_not_evalf_params().series_func(mZeta_series)); } // namespace GiNaC