X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=54d12c9b511060c2b53b6d28ca332200dd56c0ca;hb=dfaba64cff32f7dfdd96c1a96d1744ce6c1d80ad;hp=975a81e685ddb18863f4ce6f0625455aef0d36b4;hpb=5acfced77964d811d66f83e83b8fa06d77f3edd0;p=ginac.git diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index 975a81e6..54d12c9b 100644 --- a/ginac/inifcns_nstdsums.cpp +++ b/ginac/inifcns_nstdsums.cpp @@ -25,7 +25,7 @@ * 0, 1 and -1 --- or in compactified --- a string with zeros in front of 1 or -1 is written as a single * number --- notation. * - * - All functions can be nummerically evaluated with arguments in the whole complex plane. The parameters + * - All functions can be numerically evaluated with arguments in the whole complex plane. The parameters * for Li, zeta and S must be positive integers. 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. * @@ -47,7 +47,7 @@ */ /* - * GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2015 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,7 +102,7 @@ namespace { // lookup table for factors built from Bernoulli numbers // see fill_Xn() -std::vector > Xn; +std::vector> Xn; // initial size of Xn that should suffice for 32bit machines (must be even) const int xninitsizestep = 26; int xninitsize = xninitsizestep; @@ -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 @@ -1020,9 +1020,8 @@ G_do_hoelder(std::vector x, /* yes, it's passed by value */ std::vector qlsts; for (std::size_t j = r; j >= 1; --j) { qlstx.push_back(cln::cl_N(1) - x[j-1]); - if (instanceof(x[j-1], cln::cl_R_ring) && - realpart(x[j-1]) > 1 && realpart(x[j-1]) <= 2) { - qlsts.push_back(s[j-1]); + if (instanceof(x[j-1], cln::cl_R_ring) && realpart(x[j-1]) > 1) { + qlsts.push_back(1); } else { qlsts.push_back(-s[j-1]); } @@ -1044,24 +1043,43 @@ G_do_hoelder(std::vector x, /* yes, it's passed by value */ return result; } +class less_object_for_cl_N +{ +public: + bool operator() (const cln::cl_N & a, const cln::cl_N & b) const + { + // absolute value? + if (abs(a) != abs(b)) + return (abs(a) < abs(b)) ? true : false; + + // complex phase? + if (phase(a) != phase(b)) + return (phase(a) < phase(b)) ? true : false; + + // equal, therefore "less" is not true + return false; + } +}; + + // convergence transformation, used for numerical evaluation of G function. // 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; + typedef std::multimap sortmap_t; sortmap_t sortmap; std::size_t size = 0; for (std::size_t i = 0; i < x.size(); ++i) { if (!zerop(x[i])) { - sortmap.insert(std::make_pair(abs(x[i]), i)); + sortmap.insert(std::make_pair(x[i], i)); ++size; } } // include upper limit (scale) - sortmap.insert(std::make_pair(abs(y), x.size())); + sortmap.insert(std::make_pair(y, x.size())); // generate missing dummy-symbols int i = 1; @@ -1116,7 +1134,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(); @@ -1136,6 +1154,7 @@ G_numeric(const std::vector& x, const std::vector& s, // check for convergence and necessary accelerations bool need_trafo = false; bool need_hoelder = false; + bool have_trailing_zero = false; std::size_t depth = 0; for (std::size_t i = 0; i < x.size(); ++i) { if (!zerop(x[i])) { @@ -1149,19 +1168,21 @@ G_numeric(const std::vector& x, const std::vector& s, need_hoelder = true; } } - if (zerop(x[x.size() - 1])) + if (zerop(x.back())) { + have_trailing_zero = true; need_trafo = true; + } 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; @@ -1228,7 +1249,7 @@ ex mLi_numeric(const lst& m, const lst& x) static ex G2_evalf(const ex& x_, const ex& y) { - if (!y.info(info_flags::positive)) { + if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) { return G(x_, y).hold(); } lst x = is_a(x_) ? ex_to(x_) : lst(x_); @@ -1271,7 +1292,7 @@ static ex G2_eval(const ex& x_, const ex& y) { //TODO eval to MZV or H or S or Lin - if (!y.info(info_flags::positive)) { + if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) { return G(x_, y).hold(); } lst x = is_a(x_) ? ex_to(x_) : lst(x_); @@ -1320,10 +1341,10 @@ static ex G2_eval(const ex& x_, const ex& y) } +// option do_not_evalf_params() removed. unsigned G2_SERIAL::serial = function::register_new(function_options("G", 2). evalf_func(G2_evalf). eval_func(G2_eval). - do_not_evalf_params(). overloaded(2)); //TODO // derivative_func(G2_deriv). @@ -1332,7 +1353,7 @@ unsigned G2_SERIAL::serial = function::register_new(function_options("G", 2). static ex G3_evalf(const ex& x_, const ex& s_, const ex& y) { - if (!y.info(info_flags::positive)) { + if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) { return G(x_, s_, y).hold(); } lst x = is_a(x_) ? ex_to(x_) : lst(x_); @@ -1396,7 +1417,7 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y) { //TODO eval to MZV or H or S or Lin - if (!y.info(info_flags::positive)) { + if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) { return G(x_, s_, y).hold(); } lst x = is_a(x_) ? ex_to(x_) : lst(x_); @@ -1466,10 +1487,12 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y) } +// option do_not_evalf_params() removed. +// This is safe: in the code above it only matters if s_ > 0 or s_ < 0, +// s_ is allowed to be of floating type. unsigned G3_SERIAL::serial = function::register_new(function_options("G", 3). evalf_func(G3_evalf). eval_func(G3_eval). - do_not_evalf_params(). overloaded(2)); //TODO // derivative_func(G3_deriv). @@ -1639,9 +1662,8 @@ static ex Li_series(const ex& m, const ex& x, const relational& rel, int order, { if (is_a(m) || is_a(x)) { // multiple polylog - epvector seq; - seq.push_back(expair(Li(m, x), 0)); - return pseries(rel, seq); + epvector seq { expair(Li(m, x), 0) }; + return pseries(rel, std::move(seq)); } // classical polylog @@ -1657,9 +1679,8 @@ static ex Li_series(const ex& m, const ex& x, const relational& rel, int order, // substitute the argument's series expansion ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern); // maybe that was terminating, so add a proper order term - epvector nseq; - nseq.push_back(expair(Order(_ex1), order)); - ser += pseries(rel, nseq); + epvector nseq { expair(Order(_ex1), order) }; + ser += pseries(rel, std::move(nseq)); // reexpanding it will collapse the series again return ser.series(rel, order); } @@ -1758,7 +1779,7 @@ namespace { // lookup table for special Euler-Zagier-Sums (used for S_n,p(x)) // see fill_Yn() -std::vector > Yn; +std::vector> Yn; int ynsize = 0; // number of Yn[] int ynlength = 100; // initial length of all Yn[i] @@ -2212,9 +2233,8 @@ static ex S_series(const ex& n, const ex& p, const ex& x, const relational& rel, // substitute the argument's series expansion ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern); // maybe that was terminating, so add a proper order term - epvector nseq; - nseq.push_back(expair(Order(_ex1), order)); - ser += pseries(rel, nseq); + epvector nseq { expair(Order(_ex1), order) }; + ser += pseries(rel, std::move(nseq)); // reexpanding it will collapse the series again return ser.series(rel, order); } @@ -3407,9 +3427,8 @@ static ex H_eval(const ex& m_, const ex& x) static ex H_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options) { - epvector seq; - seq.push_back(expair(H(m, x), 0)); - return pseries(rel, seq); + epvector seq { expair(H(m, x), 0) }; + return pseries(rel, std::move(seq)); } @@ -3528,7 +3547,7 @@ static void initcX(std::vector& crX, int Sm = 0; int Smp1 = 0; - std::vector > crG(s.size() - 1, std::vector(L2 + 1)); + std::vector> crG(s.size() - 1, std::vector(L2 + 1)); for (int m=0; m < (int)s.size() - 1; m++) { Sm += s[m]; Smp1 = Sm + s[m+1]; @@ -3568,12 +3587,12 @@ static cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk, // [Cra] section 4 -static void calc_f(std::vector >& f_kj, +static void calc_f(std::vector>& f_kj, const int maxr, const int L1) { cln::cl_N t0, t1, t2, t3, t4; int i, j, k; - std::vector >::iterator it = f_kj.begin(); + std::vector>::iterator it = f_kj.begin(); cln::cl_F one = cln::cl_float(1, cln::float_format(Digits)); t0 = cln::exp(-lambda); @@ -3597,7 +3616,7 @@ static void calc_f(std::vector >& f_kj, // [Cra] (3.1) static cln::cl_N crandall_Z(const std::vector& s, - const std::vector >& f_kj) + const std::vector>& f_kj) { const int j = s.size(); @@ -3674,7 +3693,7 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector& s) } } - std::vector > f_kj(L1); + std::vector> f_kj(L1); calc_f(f_kj, maxr, L1); const cln::cl_N r0factorial = cln::factorial(r[0]-1);