+ shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms);
}
+// handles the transformations and the numerical evaluation of G
+// the parameter x, s and y must only contain numerics
+ex G_numeric(const lst& x, const lst& s, const ex& y);
+
+// do acceleration transformation (hoelder convolution [BBB])
+// the parameter x, s and y must only contain numerics
+ex G_do_hoelder(const lst& x, const lst& s, const ex& y)
+{
+ ex result;
+ const int size = x.nops();
+ lst newx;
+ for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
+ newx.append(*it / y);
+ }
+
+ for (int r=0; r<=size; ++r) {
+ ex buffer = pow(-1, r);
+ ex p = 2;
+ bool adjustp;
+ do {
+ adjustp = false;
+ for (lst::const_iterator it = newx.begin(); it != newx.end(); ++it) {
+ if (*it == 1/p) {
+ p += (3-p)/2;
+ adjustp = true;
+ continue;
+ }
+ }
+ } while (adjustp);
+ ex q = p / (p-1);
+ lst qlstx;
+ lst qlsts;
+ for (int j=r; j>=1; --j) {
+ qlstx.append(1-newx.op(j-1));
+ if (newx.op(j-1).info(info_flags::real) && newx.op(j-1) > 1 && newx.op(j-1) <= 2) {
+ qlsts.append( s.op(j-1));
+ } else {
+ qlsts.append( -s.op(j-1));
+ }
+ }
+ if (qlstx.nops() > 0) {
+ buffer *= G_numeric(qlstx, qlsts, 1/q);
+ }
+ lst plstx;
+ lst plsts;
+ for (int j=r+1; j<=size; ++j) {
+ plstx.append(newx.op(j-1));
+ plsts.append(s.op(j-1));
+ }
+ if (plstx.nops() > 0) {
+ buffer *= G_numeric(plstx, plsts, 1/p);
+ }
+ result += buffer;
+ }
+ return result;
+}
+
+// convergence transformation, used for numerical evaluation of G function.
+// the parameter x, s and y must only contain numerics
+static ex G_do_trafo(const lst& x, const lst& s, const ex& y)
+{
+ // sort (|x|<->position) to determine indices
+ std::multimap<ex,int> sortmap;
+ int size = 0;
+ for (int i=0; i<x.nops(); ++i) {
+ if (!x[i].is_zero()) {
+ sortmap.insert(std::pair<ex,int>(abs(x[i]), i));
+ ++size;
+ }
+ }
+ // include upper limit (scale)
+ sortmap.insert(std::pair<ex,int>(abs(y), x.nops()));
+
+ // generate missing dummy-symbols
+ int i = 1;
+ // holding dummy-symbols for the G/Li transformations
+ exvector gsyms;
+ gsyms.push_back(symbol("GSYMS_ERROR"));
+ ex lastentry;
+ for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
+ if (it != sortmap.begin()) {
+ if (it->second < x.nops()) {
+ if (x[it->second] == lastentry) {
+ gsyms.push_back(gsyms.back());
+ continue;
+ }
+ } else {
+ if (y == lastentry) {
+ gsyms.push_back(gsyms.back());
+ continue;
+ }
+ }
+ }
+ std::ostringstream os;
+ os << "a" << i;
+ gsyms.push_back(symbol(os.str()));
+ ++i;
+ if (it->second < x.nops()) {
+ lastentry = x[it->second];
+ } else {
+ lastentry = y;
+ }
+ }
+
+ // fill position data according to sorted indices and prepare substitution list
+ Gparameter a(x.nops());
+ lst subslst;
+ int pos = 1;
+ int scale;
+ for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
+ if (it->second < x.nops()) {
+ if (s[it->second] > 0) {
+ a[it->second] = pos;
+ } else {
+ a[it->second] = -pos;
+ }
+ subslst.append(gsyms[pos] == x[it->second]);
+ } else {
+ scale = pos;
+ subslst.append(gsyms[pos] == y);
+ }
+ ++pos;
+ }
+
+ // do transformation
+ Gparameter pendint;
+ ex result = G_transform(pendint, a, scale, gsyms);
+ // replace dummy symbols with their values
+ result = result.eval().expand();
+ result = result.subs(subslst).evalf();
+
+ return result;
+}
// handles the transformations and the numerical evaluation of G
// the parameter x, s and y must only contain numerics
}
// do acceleration transformation (hoelder convolution [BBB])
- if (need_hoelder) {
-
- ex result;
- const int size = x.nops();
- lst newx;
- for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
- newx.append(*it / y);
- }
-
- for (int r=0; r<=size; ++r) {
- ex buffer = pow(-1, r);
- ex p = 2;
- bool adjustp;
- do {
- adjustp = false;
- for (lst::const_iterator it = newx.begin(); it != newx.end(); ++it) {
- if (*it == 1/p) {
- p += (3-p)/2;
- adjustp = true;
- continue;
- }
- }
- } while (adjustp);
- ex q = p / (p-1);
- lst qlstx;
- lst qlsts;
- for (int j=r; j>=1; --j) {
- qlstx.append(1-newx.op(j-1));
- if (newx.op(j-1).info(info_flags::real) && newx.op(j-1) > 1 && newx.op(j-1) <= 2) {
- qlsts.append( s.op(j-1));
- } else {
- qlsts.append( -s.op(j-1));
- }
- }
- if (qlstx.nops() > 0) {
- buffer *= G_numeric(qlstx, qlsts, 1/q);
- }
- lst plstx;
- lst plsts;
- for (int j=r+1; j<=size; ++j) {
- plstx.append(newx.op(j-1));
- plsts.append(s.op(j-1));
- }
- if (plstx.nops() > 0) {
- buffer *= G_numeric(plstx, plsts, 1/p);
- }
- result += buffer;
- }
- return result;
- }
+ if (need_hoelder)
+ return G_do_hoelder(x, s, y);
// convergence transformation
- if (need_trafo) {
-
- // sort (|x|<->position) to determine indices
- std::multimap<ex,int> sortmap;
- int size = 0;
- for (int i=0; i<x.nops(); ++i) {
- if (!x[i].is_zero()) {
- sortmap.insert(std::pair<ex,int>(abs(x[i]), i));
- ++size;
- }
- }
- // include upper limit (scale)
- sortmap.insert(std::pair<ex,int>(abs(y), x.nops()));
-
- // generate missing dummy-symbols
- int i = 1;
- // holding dummy-symbols for the G/Li transformations
- exvector gsyms;
- gsyms.push_back(symbol("GSYMS_ERROR"));
- ex lastentry;
- for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
- if (it != sortmap.begin()) {
- if (it->second < x.nops()) {
- if (x[it->second] == lastentry) {
- gsyms.push_back(gsyms.back());
- continue;
- }
- } else {
- if (y == lastentry) {
- gsyms.push_back(gsyms.back());
- continue;
- }
- }
- }
- std::ostringstream os;
- os << "a" << i;
- gsyms.push_back(symbol(os.str()));
- ++i;
- if (it->second < x.nops()) {
- lastentry = x[it->second];
- } else {
- lastentry = y;
- }
- }
-
- // fill position data according to sorted indices and prepare substitution list
- Gparameter a(x.nops());
- lst subslst;
- int pos = 1;
- int scale;
- for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
- if (it->second < x.nops()) {
- if (s[it->second] > 0) {
- a[it->second] = pos;
- } else {
- a[it->second] = -pos;
- }
- subslst.append(gsyms[pos] == x[it->second]);
- } else {
- scale = pos;
- subslst.append(gsyms[pos] == y);
- }
- ++pos;
- }
-
- // do transformation
- Gparameter pendint;
- ex result = G_transform(pendint, a, scale, gsyms);
- // replace dummy symbols with their values
- result = result.eval().expand();
- result = result.subs(subslst).evalf();
-
- return result;
- }
+ if (need_trafo)
+ return G_do_trafo(x, s, y);
// do summation
lst newx;
// parameters and data for [Cra] algorithm
const cln::cl_N lambda = cln::cl_N("319/320");
-int L1;
-int L2;
-std::vector<std::vector<cln::cl_N> > f_kj;
-std::vector<cln::cl_N> crB;
-std::vector<std::vector<cln::cl_N> > crG;
-std::vector<cln::cl_N> crX;
-
void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
{
// [Cra] section 4
-void initcX(const std::vector<int>& s)
+static void initcX(std::vector<cln::cl_N>& crX,
+ const std::vector<int>& s,
+ const int L2)
{
- const int k = s.size();
-
- crX.clear();
- crG.clear();
- crB.clear();
-
- for (int i=0; i<=L2; i++) {
- crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
- }
+ std::vector<cln::cl_N> crB(L2 + 1);
+ for (int i=0; i<=L2; i++)
+ crB[i] = bernoulli(i).to_cl_N() / cln::factorial(i);
int Sm = 0;
int Smp1 = 0;
- for (int m=0; m<k-1; m++) {
- std::vector<cln::cl_N> crGbuf;
- Sm = Sm + s[m];
+ std::vector<std::vector<cln::cl_N> > crG(s.size() - 1, std::vector<cln::cl_N>(L2 + 1));
+ for (int m=0; m < s.size() - 1; m++) {
+ Sm += s[m];
Smp1 = Sm + s[m+1];
- for (int i=0; i<=L2; i++) {
- crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
- }
- crG.push_back(crGbuf);
+ for (int i = 0; i <= L2; i++)
+ crG[m][i] = cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2);
}
crX = crB;
- for (int m=0; m<k-1; m++) {
- std::vector<cln::cl_N> Xbuf;
- for (int i=0; i<=L2; i++) {
- Xbuf.push_back(crX[i] * crG[m][i]);
- }
+ for (std::size_t m = 0; m < s.size() - 1; m++) {
+ std::vector<cln::cl_N> Xbuf(L2 + 1);
+ for (int i = 0; i <= L2; i++)
+ Xbuf[i] = crX[i] * crG[m][i];
+
halfcyclic_convolute(Xbuf, crB, crX);
}
}
// [Cra] section 4
-cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
+static cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk,
+ const std::vector<cln::cl_N>& crX)
{
cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
cln::cl_N factor = cln::expt(lambda, Sqk);
// [Cra] section 4
-void calc_f(int maxr)
+static void calc_f(std::vector<std::vector<cln::cl_N> >& f_kj,
+ const int maxr, const int L1)
{
- f_kj.clear();
- f_kj.resize(L1);
-
cln::cl_N t0, t1, t2, t3, t4;
int i, j, k;
std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
// [Cra] (3.1)
-cln::cl_N crandall_Z(const std::vector<int>& s)
+static cln::cl_N crandall_Z(const std::vector<int>& s,
+ const std::vector<std::vector<cln::cl_N> >& f_kj)
{
const int j = s.size();
std::vector<int> r = s;
const int j = r.size();
+ std::size_t L1;
+
// decide on maximal size of f_kj for crandall_Z
if (Digits < 50) {
L1 = 150;
L1 = Digits * 3 + j*2;
}
+ std::size_t L2;
// decide on maximal size of crX for crandall_Y
if (Digits < 38) {
L2 = 63;
}
}
- calc_f(maxr);
+ std::vector<std::vector<cln::cl_N> > f_kj(L1);
+ calc_f(f_kj, maxr, L1);
const cln::cl_N r0factorial = cln::factorial(r[0]-1);
Srun -= skp1buf;
r.pop_back();
- initcX(r);
+ std::vector<cln::cl_N> crX;
+ initcX(crX, r, L2);
for (int q=0; q<skp1buf; q++) {
- cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
- cln::cl_N pp2 = crandall_Z(rz);
+ cln::cl_N pp1 = crandall_Y_loop(Srun+q-k, crX);
+ cln::cl_N pp2 = crandall_Z(rz, f_kj);
rz.front()--;
}
rz.insert(rz.begin(), r.back());
- initcX(rz);
+ std::vector<cln::cl_N> crX;
+ initcX(crX, rz, L2);
- res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
+ res = (res + crandall_Y_loop(S-j, crX)) / r0factorial
+ + crandall_Z(rz, f_kj);
return res;
}
git://www.ginac.de / ginac.git / blobdiff