Bug fix in the routine H_evalf (a minus sign should not be forgotten).

[ginac.git] / ginac / inifcns_nstdsums.cpp

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index 0a60d780ded9f427703ee4942d8c8a2fb8510bb2..51ba0641447ef7e33682818cb22beafdeeb9e6ce 100644 (file)
--- a/ginac/inifcns_nstdsums.cpp
+++ b/ginac/inifcns_nstdsums.cpp
@@ -47,7 +47,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2019 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2021 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
@@ -1013,21 +1013,25 @@ G_do_hoelder(std::vector<cln::cl_N> x, /* yes, it's passed by value */
        for (std::size_t i = 0; i < size; ++i)
                x[i] = x[i]/y;
 
+        // 24.03.2021: this block can be outside the loop over r 
+       cln::cl_RA p(2);
+       bool adjustp;
+       do {
+               adjustp = false;
+               for (std::size_t i = 0; i < size; ++i) {
+                        // 24.03.2021: replaced (x[i] == cln::cl_RA(1)/p) by (cln::zerop(x[i] - cln::cl_RA(1)/p)
+                        //             in the case where we compare a float with a rational, CLN behaves differently in the two versions
+                       if (cln::zerop(x[i] - cln::cl_RA(1)/p) ) {
+                               p = p/2 + cln::cl_RA(3)/2;
+                               adjustp = true;
+                               continue;
+                       }
+               }
+       } while (adjustp);
+       cln::cl_RA q = p/(p-1);
+
        for (std::size_t r = 0; r <= size; ++r) {
                cln::cl_N buffer(1 & r ? -1 : 1);
-               cln::cl_RA p(2);
-               bool adjustp;
-               do {
-                       adjustp = false;
-                       for (std::size_t i = 0; i < size; ++i) {
-                               if (x[i] == cln::cl_RA(1)/p) {
-                                       p = p/2 + cln::cl_RA(3)/2;
-                                       adjustp = true;
-                                       continue;
-                               }
-                       }
-               } while (adjustp);
-               cln::cl_RA q = p/(p-1);
                std::vector<cln::cl_N> qlstx;
                std::vector<int> qlsts;
                for (std::size_t j = r; j >= 1; --j) {
@@ -3223,7 +3227,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) {
@@ -3236,7 +3239,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);
                        }
@@ -3297,7 +3299,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
@@ -3309,7 +3318,9 @@ static ex H_evalf(const ex& x1, const ex& x2)
                        // x -> 1-x
                        if (has_minus_one) {
                                map_trafo_H_convert_to_Li filter;
-                               return filter(H(m, numeric(x)).hold()).evalf();
+                                // 09.06.2021: bug fix: don't forget a possible minus sign from the case realpart(x) < 0
+                               res *= filter(H(m, numeric(x)).hold()).evalf();
+                               return res;
                        }
                        map_trafo_H_1mx trafo;
                        res *= trafo(H(m, xtemp).hold());
@@ -3882,17 +3893,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
@@ -3900,6 +3915,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;
@@ -4053,7 +4072,8 @@ static ex zeta2_evalf(const ex& x, const ex& s)
                return numeric(zeta_do_Hoelder_convolution(xi, si));
        }
 
-       return zeta(x, s).hold();
+       // x and s are not lists: convert to lists
+       return zeta(lst{x}, lst{s}).evalf();
 }