Fix bug in H_evalf: Flag has_minus_one is now computed where it is needed.

[ginac.git] / ginac / inifcns_nstdsums.cpp

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index e21e311d2fd44c6b19a3ebbfb14ff61dcf4461dc..c20526bafd5a38fb301d9470e8c4138ba8540d27 100644 (file)
--- 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 <sstream>
 #include <stdexcept>
 #include <vector>
+#include <cmath>
 
 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<cln::cl_N> x, /* yes, it's passed by value */
                std::vector<int> 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<int>& 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<int>& 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<int>& 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<lst>(x);
                std::vector<int> r(count);
+               std::vector<int> 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<numeric>(*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;