X-Git-Url: https://ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=94d0fb0c1d3af9442977375b4f90120c972a0b1f;hb=f4f7f966d313edf65d46216c7eb48d0791c6537c;hp=4fe91ad391eeb7f41e40362c7fa67223339da3a6;hpb=9401a14b25939a415575aaf77aa824e90c92d0a9;p=ginac.git diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index 4fe91ad3..94d0fb0c 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. * @@ -46,7 +46,7 @@ */ /* - * 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 @@ -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,17 +254,23 @@ 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; + 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; } @@ -244,6 +297,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 +307,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; } @@ -450,6 +508,15 @@ static ex Li_evalf(const ex& x1, const ex& x2) if (is_a(x1) && is_a(x2)) { return Li_num(ex_to(x1).to_int(), ex_to(x2)); } + if (is_a(x1) && !is_a(x2)) { + // try to numerically evaluate second argument + ex x2_val = x2.evalf(); + if (is_a(x2_val)) { + return Li_num(ex_to(x1).to_int(), ex_to(x2_val)); + } else { + return Li(x1, x2).hold(); + } + } // multiple polylogs else if (is_a(x1) && is_a(x2)) { ex conv = 1; @@ -939,9 +1006,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); @@ -1002,8 +1068,15 @@ numeric S_num(int n, int p, const numeric& x) 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(x)) { - return S_num(ex_to(n).to_int(), ex_to(p).to_int(), ex_to(x)); + if (n.info(info_flags::posint) && p.info(info_flags::posint)) { + if (is_a(x)) { + return S_num(ex_to(n).to_int(), ex_to(p).to_int(), ex_to(x)); + } else { + ex x_val = x.evalf(); + if (is_a(x_val)) { + return S_num(ex_to(n).to_int(), ex_to(p).to_int(), ex_to(x_val)); + } + } } return S(n, p, x).hold(); } @@ -1029,6 +1102,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(); } @@ -1088,6 +1165,10 @@ REGISTER_FUNCTION(S, // 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 @@ -1601,7 +1682,7 @@ struct map_trafo_H_1overx : public map_function } 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()); } } @@ -1836,18 +1917,27 @@ cln::cl_N H_do_sum(const std::vector& m, const cln::cl_N& x) static ex H_evalf(const ex& x1, const ex& x2) { - if (is_a(x1) && is_a(x2)) { + if (is_a(x1)) { + + cln::cl_N x; + if (is_a(x2)) { + x = ex_to(x2).to_cl_N(); + } else { + ex x2_val = x2.evalf(); + if (is_a(x2_val)) { + x = ex_to(x2_val).to_cl_N(); + } + } + for (int i=0; i(x2).to_cl_N(); - const lst& morg = ex_to(x1); // remove trailing zeros ... if (*(--morg.end()) == 0) { @@ -1928,6 +2018,11 @@ static ex H_evalf(const ex& x1, const ex& x2) // 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 @@ -2033,25 +2128,25 @@ static ex H_eval(const ex& m_, const ex& x) if ((x == _ex1) && (*(--m.end()) != _ex0)) { return convert_H_to_zeta(m); } -// if (step == 0) { -// if (pos1 == _ex0) { -// // all zero -// if (x == _ex0) { -// return H(m_, x).hold(); -// } -// return pow(log(x), m.nops()) / factorial(m.nops()); -// } else { -// // all (minus) one -// return pow(-pos1*log(1-pos1*x), m.nops()) / factorial(m.nops()); -// } -// } else if ((step == 1) && (pos1 == _ex0)){ -// // convertible to S -// if (pos2 == _ex1) { -// return S(n, p, x); -// } else { -// return pow(-1, p) * S(n, p, -x); -// } -// } + if (step == 0) { + if (pos1 == _ex0) { + // all zero + if (x == _ex0) { + return H(m_, x).hold(); + } + return pow(log(x), m.nops()) / factorial(m.nops()); + } else { + // all (minus) one + return pow(-pos1*log(1-pos1*x), m.nops()) / factorial(m.nops()); + } + } else if ((step == 1) && (pos1 == _ex0)){ + // convertible to S + if (pos2 == _ex1) { + return S(n, p, x); + } else { + return pow(-1, p) * S(n, p, -x); + } + } if (x == _ex0) { return _ex0; } @@ -2623,7 +2718,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). @@ -2760,7 +2855,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).