From 40b0071bbe4055bbae533634f008af38503ab696 Mon Sep 17 00:00:00 2001 From: Stefan Weinzierl Date: Sun, 12 Oct 2014 18:24:07 +0000 Subject: [PATCH] Fix evaluation of some G-polylogs. This patch corrects and supersedes the patches 97ef604e "G_numeric: fix numeric evaluation with trailing zeros and y != 1." and 9e80b0d3 "G_numeric: fix evaluation with y == 1". The original motivation for 97ef604e was to make sure that Hoelder convolution is only used if there no trailing zeros. Patch 97ef604e delegated the case of trailing zeros to G_do_trafo, which correctly removes the trailing zeros, but in addition also performed the transformation described in sect. 5.3 of hep-ph/0410259 ("convergence transformation"). The inappropriate call of the convergence transformation is the cause of the new bug reported at http://www.cebix.net/pipermail/ginac-list/2014-September/002011.html The patch 9e80b0d3 cures the symptons mentioned in the above bug report, but fails for other cases, like evalf( G({-2,2,-2,0},2) ). What should be done is the following: If trailing zeros are detected in G_numeric, these should be removed and the result should be returned to G_numeric. The routine G_numeric decides then what to do next: either calling the convergence transformation, or Hoelder convolution or direct summation. What is needed is a subroutine, which just removes trailing zeros, but does not perform the convergence transformation. With the present code the minimal modification to achieve this goal is to add an additional boolean parameter flag_trailing_zeros_only to G_do_trafo (and its dependent sub-routines), so G_do_trafo can be called for the removal of trailing zeros only. This patch implements this and uses G_do_trafo to remove trailing zeros only for the case at hand. --- ginac/inifcns_nstdsums.cpp | 61 ++++++++++++++++++-------------------- 1 file changed, 29 insertions(+), 32 deletions(-) diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index 977642d9..866472f4 100644 --- a/ginac/inifcns_nstdsums.cpp +++ b/ginac/inifcns_nstdsums.cpp @@ -803,12 +803,12 @@ ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, i // forward declaration ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2, const Gparameter& pendint, const Gparameter& a_old, int scale, - const exvector& gsyms); + const exvector& gsyms, bool flag_trailing_zeros_only); // G transformation [VSW] ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, - const exvector& gsyms) + const exvector& gsyms, bool flag_trailing_zeros_only) { // main recursion routine // @@ -848,18 +848,18 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, if (trailing_zeros > 0) { ex result; Gparameter new_a(a.begin(), a.end()-1); - result += G_eval1(0, scale, gsyms) * G_transform(pendint, new_a, scale, gsyms); + result += G_eval1(0, scale, gsyms) * G_transform(pendint, new_a, scale, gsyms, flag_trailing_zeros_only); for (Gparameter::const_iterator it = a.begin(); it != firstzero; ++it) { Gparameter new_a(a.begin(), it); new_a.push_back(0); new_a.insert(new_a.end(), it, a.end()-1); - result -= G_transform(pendint, new_a, scale, gsyms); + result -= G_transform(pendint, new_a, scale, gsyms, flag_trailing_zeros_only); } return result / trailing_zeros; } - // convergence case - if (convergent) { + // convergence case or flag_trailing_zeros_only + if (convergent || flag_trailing_zeros_only) { if (pendint.size() > 0) { return G_eval(convert_pending_integrals_G(pendint), pendint.front(), gsyms)* @@ -886,10 +886,10 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, Gparameter a1(a.begin(),min_it+1); Gparameter a2(min_it+1,a.end()); - ex result = G_transform(pendint, a2, scale, gsyms)* - G_transform(empty, a1, scale, gsyms); + ex result = G_transform(pendint, a2, scale, gsyms, flag_trailing_zeros_only)* + G_transform(empty, a1, scale, gsyms, flag_trailing_zeros_only); - result -= shuffle_G(empty, a1, a2, pendint, a, scale, gsyms); + result -= shuffle_G(empty, a1, a2, pendint, a, scale, gsyms, flag_trailing_zeros_only); return result; } @@ -900,7 +900,7 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, Gparameter new_pendint = prepare_pending_integrals(pendint, a[min_it_pos]); Gparameter new_a = a; new_a[min_it_pos] = 0; - ex result = G_transform(empty, new_a, scale, gsyms); + ex result = G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); if (pendint.size() > 0) { result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front(), gsyms); @@ -914,31 +914,31 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, new_pendint.push_back(*changeit); result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front(), gsyms)* - G_transform(empty, new_a, scale, gsyms); + G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); int buffer = *changeit; *changeit = *min_it; - result += G_transform(new_pendint, new_a, scale, gsyms); + result += G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only); *changeit = buffer; new_pendint.pop_back(); --changeit; new_pendint.push_back(*changeit); result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front(), gsyms)* - G_transform(empty, new_a, scale, gsyms); + G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); *changeit = *min_it; - result -= G_transform(new_pendint, new_a, scale, gsyms); + result -= G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only); } else { // smallest at the front new_pendint.push_back(scale); result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front(), gsyms)* - G_transform(empty, new_a, scale, gsyms); + G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); new_pendint.back() = *changeit; result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front(), gsyms)* - G_transform(empty, new_a, scale, gsyms); + G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); *changeit = *min_it; - result += G_transform(new_pendint, new_a, scale, gsyms); + result += G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only); } return result; } @@ -948,27 +948,27 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, // for the one that is equal to a_old ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2, const Gparameter& pendint, const Gparameter& a_old, int scale, - const exvector& gsyms) + const exvector& gsyms, bool flag_trailing_zeros_only) { if (a1.size()==0 && a2.size()==0) { // veto the one configuration we don't want if ( a0 == a_old ) return 0; - return G_transform(pendint, a0, scale, gsyms); + return G_transform(pendint, a0, scale, gsyms, flag_trailing_zeros_only); } if (a2.size()==0) { Gparameter empty; Gparameter aa0 = a0; aa0.insert(aa0.end(),a1.begin(),a1.end()); - return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms); + return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms, flag_trailing_zeros_only); } if (a1.size()==0) { Gparameter empty; Gparameter aa0 = a0; aa0.insert(aa0.end(),a2.begin(),a2.end()); - return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms); + return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms, flag_trailing_zeros_only); } Gparameter a1_removed(a1.begin()+1,a1.end()); @@ -980,8 +980,8 @@ ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2 a01.push_back( a1[0] ); a02.push_back( a2[0] ); - return shuffle_G(a01, a1_removed, a2, pendint, a_old, scale, gsyms) - + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms); + return shuffle_G(a01, a1_removed, a2, pendint, a_old, scale, gsyms, flag_trailing_zeros_only) + + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms, flag_trailing_zeros_only); } // handles the transformations and the numerical evaluation of G @@ -1047,7 +1047,7 @@ G_do_hoelder(std::vector x, /* yes, it's passed by value */ // the parameter x, s and y must only contain numerics static cln::cl_N G_do_trafo(const std::vector& x, const std::vector& s, - const cln::cl_N& y) + const cln::cl_N& y, bool flag_trailing_zeros_only) { // sort (|x|<->position) to determine indices typedef std::multimap sortmap_t; @@ -1115,7 +1115,7 @@ G_do_trafo(const std::vector& x, const std::vector& s, // do transformation Gparameter pendint; - ex result = G_transform(pendint, a, scale, gsyms); + ex result = G_transform(pendint, a, scale, gsyms, flag_trailing_zeros_only); // replace dummy symbols with their values result = result.eval().expand(); result = result.subs(subslst).evalf(); @@ -1149,24 +1149,21 @@ G_numeric(const std::vector& x, const std::vector& s, need_hoelder = true; } } - have_trailing_zero = zerop(x.back()); - if (have_trailing_zero) { + if (zerop(x.back())) { + have_trailing_zero = true; need_trafo = true; - if (y != 1) { - need_hoelder = false; - } } if (depth == 1 && x.size() == 2 && !need_trafo) return - Li_projection(2, y/x[1], cln::float_format(Digits)); // do acceleration transformation (hoelder convolution [BBB]) - if (need_hoelder) + if (need_hoelder && !have_trailing_zero) return G_do_hoelder(x, s, y); // convergence transformation if (need_trafo) - return G_do_trafo(x, s, y); + return G_do_trafo(x, s, y, have_trailing_zero); // do summation std::vector newx; -- 2.44.0