*/
/*
- * GiNaC Copyright (C) 1999-2018 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
#include <sstream>
#include <stdexcept>
#include <vector>
+#include <cmath>
namespace GiNaC {
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)
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) *
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) {
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]);
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) {
m.append(0);
}
m.append(-1);
- has_minus_one = true;
} else {
m.append(it);
}
}
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
// 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());
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;
}
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);
}
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);
}
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;
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
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;
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();
}
git://www.ginac.de / ginac.git / blobdiff