From 8803eb15c4067dc7378b58f6100a5d11bcc95e80 Mon Sep 17 00:00:00 2001 From: Jens Vollinga Date: Thu, 14 Aug 2003 18:47:45 +0000 Subject: [PATCH] Added new functions (polylogarithm[classical/multiple/harmonic/nielsens] and multiple zeta value) --- ginac/Makefile.am | 3 +- ginac/inifcns.h | 12 + ginac/inifcns_nstdsums.cpp | 631 +++++++++++++++++++++++++++++++++++++ 3 files changed, 645 insertions(+), 1 deletion(-) create mode 100644 ginac/inifcns_nstdsums.cpp diff --git a/ginac/Makefile.am b/ginac/Makefile.am index f0896bd4..df67f197 100644 --- a/ginac/Makefile.am +++ b/ginac/Makefile.am @@ -3,7 +3,8 @@ lib_LTLIBRARIES = libginac.la libginac_la_SOURCES = add.cpp archive.cpp basic.cpp constant.cpp ex.cpp \ expair.cpp expairseq.cpp exprseq.cpp fail.cpp function.cpp inifcns.cpp \ - inifcns_trans.cpp inifcns_zeta.cpp inifcns_gamma.cpp matrix.cpp mul.cpp \ + inifcns_trans.cpp inifcns_zeta.cpp inifcns_gamma.cpp inifcns_nstdsums.cpp \ + matrix.cpp mul.cpp \ normal.cpp numeric.cpp operators.cpp power.cpp registrar.cpp relational.cpp \ symbol.cpp pseries.cpp utils.cpp ncmul.cpp structure.cpp exprseq_suppl.cpp \ lst.cpp lst_suppl.cpp input_parser.yy input_lexer.ll input_lexer.h \ diff --git a/ginac/inifcns.h b/ginac/inifcns.h index f4c00dd0..258ee2e2 100644 --- a/ginac/inifcns.h +++ b/ginac/inifcns.h @@ -142,6 +142,18 @@ DECLARE_FUNCTION_2P(binomial) /** Order term function (for truncated power series). */ DECLARE_FUNCTION_1P(Order) +/** Polylogarithm and multiple polylogarithm. */ +DECLARE_FUNCTION_2P(Li) + +/** Nielsen's generalized polylogarithm. */ +DECLARE_FUNCTION_3P(S) + +/** Harmonic polylogarithm. */ +DECLARE_FUNCTION_2P(H) + +/** Multiple zeta value. */ +DECLARE_FUNCTION_1P(mZeta) + ex lsolve(const ex &eqns, const ex &symbols, unsigned options = solve_algo::automatic); /** Check whether a function is the Order (O(n)) function. */ diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp new file mode 100644 index 00000000..879bb639 --- /dev/null +++ b/ginac/inifcns_nstdsums.cpp @@ -0,0 +1,631 @@ +/** @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 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" + + +namespace GiNaC { + + +////////////////////// +// helper functions // +////////////////////// + + +// helper function for classical polylog Li +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; +} + + +// 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); +} + + +// helper function for classical polylog Li +static numeric Li_num(int n, const numeric& x) +{ + if (n == 1) { + // just a log + return -cln::log(1-x.to_cl_N()); + } + if (x.is_zero()) { + return 0; + } + if (x == 1) { + // [Kol] (2.22) + return cln::zeta(n); + } + else if (x == -1) { + // [Kol] (2.22) + return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n); + } + + // what is the desired float format? + // first guess: default format + cln::float_format_t prec = cln::default_float_format; + const cln::cl_N value = x.to_cl_N(); + // second guess: the argument's format + if (!x.real().is_rational()) + prec = cln::float_format(cln::the(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 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 >= 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(); +} + +REGISTER_FUNCTION(Li, eval_func(Li_eval).evalf_func(Li_evalf).do_not_evalf_params()); + + +// 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(); +} + +REGISTER_FUNCTION(S, eval_func(S_eval).evalf_func(S_evalf).do_not_evalf_params()); + + +// 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(); + } + + 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(); +} + +REGISTER_FUNCTION(H, eval_func(H_eval).evalf_func(H_evalf).do_not_evalf_params()); + + +// 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(); +} + +REGISTER_FUNCTION(mZeta, eval_func(mZeta_eval).evalf_func(mZeta_evalf).do_not_evalf_params()); + + +} // namespace GiNaC + -- 2.49.0