X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=c20526bafd5a38fb301d9470e8c4138ba8540d27;hb=0eaae44cd9eb9fa987bb9cbd4341b0f4c8d2f495;hp=e21e311d2fd44c6b19a3ebbfb14ff61dcf4461dc;hpb=bb68d77b93ad9aab1be3705e052d7091cee94d56;p=ginac.git diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index e21e311d..c20526ba 100644 --- a/ginac/inifcns_nstdsums.cpp +++ b/ginac/inifcns_nstdsums.cpp @@ -47,7 +47,7 @@ */ /* - * GiNaC Copyright (C) 1999-2016 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2019 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 @@ -83,6 +83,7 @@ #include #include #include +#include namespace GiNaC { @@ -337,10 +338,10 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr // the switching point was empirically determined. the optimal point // depends on hardware, Digits, ... so an approx value is okay. // it solves also the problem with precision due to the u=-log(1-x) transformation - if (cln::abs(cln::realpart(x)) < 0.25) { - + if (cln::abs(x) < 0.25) { return Li2_do_sum(x); } else { + // Li2_do_sum practically doesn't converge near x == ±I return Li2_do_sum_Xn(x); } } else { @@ -366,9 +367,10 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr if (cln::realpart(x) < 0.5) { // choose the faster algorithm // with n>=12 the "normal" summation always wins against the method with Xn - if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) { + if ((cln::abs(x) < 0.3) || (n >= 12)) { return Lin_do_sum(n, x); } else { + // Li2_do_sum practically doesn't converge near x == ±I return Lin_do_sum_Xn(n, x); } } else { @@ -604,6 +606,27 @@ ex G_eval(const Gparameter& a, int scale, const exvector& gsyms) return pow(-1, x.nops()) * Li(m, x); } +// convert back to standard G-function, keep information on small imaginary parts +ex G_eval_to_G(const Gparameter& a, int scale, const exvector& gsyms) +{ + lst z; + lst s; + for (const auto & it : a) { + if (it != 0) { + z.append(gsyms[std::abs(it)]); + if ( it<0 ) { + s.append(-1); + } else { + s.append(1); + } + } else { + z.append(0); + s.append(1); + } + } + return G(z,s,gsyms[std::abs(scale)]); +} + // converts data for G: pending_integrals -> a Gparameter convert_pending_integrals_G(const Gparameter& pending_integrals) @@ -843,8 +866,12 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, return result / trailing_zeros; } - // convergence case or flag_trailing_zeros_only - if (convergent || flag_trailing_zeros_only) { + // flag_trailing_zeros_only: in this case we don't have pending integrals + if (flag_trailing_zeros_only) + return G_eval_to_G(a, scale, gsyms); + + // convergence case + if (convergent) { if (pendint.size() > 0) { return G_eval(convert_pending_integrals_G(pendint), pendint.front(), gsyms) * @@ -1005,7 +1032,7 @@ 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) { + if (imagpart(x[j-1])==0 && realpart(x[j-1]) >= 1) { qlsts.push_back(1); } else { qlsts.push_back(-s[j-1]); @@ -3196,7 +3223,6 @@ static ex H_evalf(const ex& x1, const ex& x2) return filter(H(x1, xtemp).hold()).subs(xtemp==x2).evalf(); } // ... and expand parameter notation - bool has_minus_one = false; lst m; for (const auto & it : morg) { if (it > 1) { @@ -3209,7 +3235,6 @@ static ex H_evalf(const ex& x1, const ex& x2) m.append(0); } m.append(-1); - has_minus_one = true; } else { m.append(it); } @@ -3270,7 +3295,14 @@ static ex H_evalf(const ex& x1, const ex& x2) } return res.subs(xtemp == numeric(x)).evalf(); } - + + // check for letters (-1) + bool has_minus_one = false; + for (const auto & it : m) { + if (it == -1) + has_minus_one = true; + } + // check transformations for 0.95 <= |x| < 2.0 // |(1-x)/(1+x)| < 0.9 -> circular area with center=9.53+0i and radius=9.47 @@ -3566,7 +3598,7 @@ static cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk, factor = factor * lambda; N++; res = res + crX[N] * factor / (N+Sqk); - } while ((res != resbuf) || cln::zerop(crX[N])); + } while (((res != resbuf) || cln::zerop(crX[N])) && (N+1 < crX.size())); return res; } @@ -3613,7 +3645,7 @@ static cln::cl_N crandall_Z(const std::vector& s, t0buf = t0; q++; t0 = t0 + f_kj[q+j-2][s[0]-1]; - } while (t0 != t0buf); + } while ((t0 != t0buf) && (q+j-1 < f_kj.size())); return t0 / cln::factorial(s[0]-1); } @@ -3630,7 +3662,7 @@ static cln::cl_N crandall_Z(const std::vector& s, t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]); } t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1]; - } while (t[0] != t0buf); + } while ((t[0] != t0buf) && (q+j-1 < f_kj.size())); return t[0] / cln::factorial(s[0]-1); } @@ -3663,8 +3695,11 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector& s) L2 = 511; } else if (Digits < 808) { L2 = 1023; - } else { + } else if (Digits < 1636) { L2 = 2047; + } else { + // [Cra] section 6, log10(lambda/2/Pi) approx -0.79 for lambda=319/320, add some extra digits + L2 = std::pow(2, ceil( std::log2((long(Digits))/0.79 + 40 )) ) - 1; } cln::cl_N res; @@ -3852,17 +3887,21 @@ static ex zeta1_evalf(const ex& x) const int count = x.nops(); const lst& xlst = ex_to(x); std::vector r(count); + std::vector si(count); // check parameters and convert them auto it1 = xlst.begin(); auto it2 = r.begin(); + auto it_swrite = si.begin(); do { if (!(*it1).info(info_flags::posint)) { return zeta(x).hold(); } *it2 = ex_to(*it1).to_int(); + *it_swrite = 1; it1++; it2++; + it_swrite++; } while (it2 != r.end()); // check for divergence @@ -3870,6 +3909,10 @@ static ex zeta1_evalf(const ex& x) return zeta(x).hold(); } + // use Hoelder convolution if Digits is large + if (Digits>50) + return numeric(zeta_do_Hoelder_convolution(r, si)); + // decide on summation algorithm // this is still a bit clumsy int limit = (Digits>17) ? 10 : 6;