multiple zeta values: make crandall_Y_loop helper function reentrant.

[ginac.git] / ginac / inifcns_nstdsums.cpp

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index dd8625a9713a076b89a6c3f665666d9e280506cb..d85044d04bd1f311b2d26ad5ab8ccf0572d1ebe7 100644 (file)
--- a/ginac/inifcns_nstdsums.cpp
+++ b/ginac/inifcns_nstdsums.cpp
@@ -1219,7 +1219,12 @@ static ex G2_evalf(const ex& x_, const ex& y)
                if (*it != _ex0) {
                        all_zero = false;
                }
-               s.append(+1);
+               if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
+                       s.append(-1);
+               }
+               else {
+                       s.append(+1);
+               }
        }
        if (all_zero) {
                return pow(log(y), x.nops()) / factorial(x.nops());
@@ -1255,7 +1260,12 @@ static ex G2_eval(const ex& x_, const ex& y)
                if (*it != _ex0) {
                        all_zero = false;
                }
-               s.append(+1);
+               if ( !ex_to<numeric>(*it).is_real() && ex_to<numeric>(*it).imag() < 0 ) {
+                       s.append(-1);
+               }
+               else {
+                       s.append(+1);
+               }
        }
        if (all_zero) {
                return pow(log(y), x.nops()) / factorial(x.nops());
@@ -1308,10 +1318,21 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
                if (*itx != _ex0) {
                        all_zero = false;
                }
-               if (*its >= 0) {
-                       sn.append(+1);
-               } else {
-                       sn.append(-1);
+               if ( ex_to<numeric>(*itx).is_real() ) {
+                       if ( *its >= 0 ) {
+                               sn.append(+1);
+                       }
+                       else {
+                               sn.append(-1);
+                       }
+               }
+               else {
+                       if ( ex_to<numeric>(*itx).imag() > 0 ) {
+                               sn.append(+1);
+                       }
+                       else {
+                               sn.append(-1);
+                       }
                }
        }
        if (all_zero) {
@@ -1355,10 +1376,21 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
                if (*itx != _ex0) {
                        all_zero = false;
                }
-               if (*its >= 0) {
-                       sn.append(+1);
-               } else {
-                       sn.append(-1);
+               if ( ex_to<numeric>(*itx).is_real() ) {
+                       if ( *its >= 0 ) {
+                               sn.append(+1);
+                       }
+                       else {
+                               sn.append(-1);
+                       }
+               }
+               else {
+                       if ( ex_to<numeric>(*itx).imag() > 0 ) {
+                               sn.append(+1);
+                       }
+                       else {
+                               sn.append(-1);
+                       }
                }
        }
        if (all_zero) {
@@ -3385,13 +3417,6 @@ namespace {
 
 // 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)
 {
@@ -3406,44 +3431,39 @@ void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln
 
 
 // [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);
@@ -3461,11 +3481,9 @@ cln::cl_N crandall_Y_loop(const cln::cl_N& 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();
@@ -3491,7 +3509,8 @@ void calc_f(int maxr)
 
 
 // [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();
 
@@ -3532,6 +3551,8 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
        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;
@@ -3539,6 +3560,7 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
                L1 = Digits * 3 + j*2;
        }
 
+       std::size_t L2;
        // decide on maximal size of crX for crandall_Y
        if (Digits < 38) {
                L2 = 63;
@@ -3565,7 +3587,8 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
                }
        }
 
-       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);
 
@@ -3579,12 +3602,13 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
                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()--;
                        
@@ -3598,9 +3622,11 @@ cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
        }
        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;
 }