From ff604e21b9108236522fb59ba5f72ba1b88b374e Mon Sep 17 00:00:00 2001 From: Stefan Weinzierl Date: Sun, 12 Jan 2014 22:31:22 +0000 Subject: [PATCH] Fix S_num for arguments close to the sixth root of unity or its conjugate. The method S_num within the Nielsen polylogs used to map the region abs(x)<=1 && abs(x)>0.95 && abs(1-x)<=1 && abs(1-x)>0.95 infinitely many times onto itself. This infinite recursion is now avoided. This however reveals the next problem: The numerical convergence in this region is very slow. Within the Nielsen polylogs there is no transformation available to improve the convergence. However we can use the (1-x)/(1+x) transformation within the harmonic polylogs. In order to avoid another infinite recursion I have inserted a few hold()'s in the method H_evalf, otherwise we would fall back immediately again to Nielsen polylogs. The hold()s should have been there anyway. --- ginac/inifcns_nstdsums.cpp | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index c33a4cf8..9876075c 100644 --- a/ginac/inifcns_nstdsums.cpp +++ b/ginac/inifcns_nstdsums.cpp @@ -2044,7 +2044,9 @@ const cln::cl_N S_num(int n, int p, const cln::cl_N& x) prec = cln::float_format(cln::the(cln::imagpart(value))); // [Kol] (5.3) - if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95))) { + // the condition abs(1-value)>1 avoids an infinite recursion in the region abs(value)<=1 && abs(value)>0.95 && abs(1-value)<=1 && abs(1-value)>0.95 + // we don't care here about abs(value)<1 && real(value)>0.5, this will be taken care of in S_projection + if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95) && (cln::abs(1-value) > 1) )) { 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); @@ -2085,6 +2087,16 @@ const cln::cl_N S_num(int n, int p, const cln::cl_N& x) return result; } + + if ((cln::abs(value) > 0.95) && (cln::abs(value-9.53) < 9.47)) { + lst m; + m.append(n+1); + for (int s=0; s(res).to_cl_N(); + } else { return S_projection(n, p, value, prec); } @@ -3230,7 +3242,7 @@ static ex H_evalf(const ex& x1, const ex& x2) // x -> 1/x if (cln::abs(x) >= 2.0) { map_trafo_H_1overx trafo; - res *= trafo(H(m, xtemp)); + res *= trafo(H(m, xtemp).hold()); if (cln::imagpart(x) <= 0) { res = res.subs(H_polesign == -I*Pi); } else { @@ -3245,7 +3257,7 @@ static ex H_evalf(const ex& x1, const ex& x2) if (cln::abs(x-9.53) <= 9.47) { // x -> (1-x)/(1+x) map_trafo_H_1mxt1px trafo; - res *= trafo(H(m, xtemp)); + res *= trafo(H(m, xtemp).hold()); } else { // x -> 1-x if (has_minus_one) { @@ -3253,7 +3265,7 @@ static ex H_evalf(const ex& x1, const ex& x2) return filter(H(m, numeric(x)).hold()).evalf(); } map_trafo_H_1mx trafo; - res *= trafo(H(m, xtemp)); + res *= trafo(H(m, xtemp).hold()); } return res.subs(xtemp == numeric(x)).evalf(); -- 2.44.0