X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=f040e8ad64df2aab86a5d8dda2d3b9c51c3e8b3d;hp=e2e48bc1eeeba031f67490d1ef6383110bf8e295;hb=HEAD;hpb=8cffcdf13d817a47f217f1a1043317d95969e070

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index e2e48bc1..e69cdb40 100644
--- 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-2024 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 {
 
@@ -605,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)
@@ -844,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) *
@@ -987,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) {
@@ -3197,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) {
@@ -3210,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);
 			}
@@ -3271,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
@@ -3283,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());
@@ -3567,7 +3604,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;
 }
 
@@ -3614,7 +3651,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);
 	}
@@ -3631,7 +3668,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);
 }
@@ -3664,8 +3701,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;
@@ -3853,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
@@ -3871,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;
@@ -4024,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();
 }