inifcns_nstdsums.cpp: Lin_numeric takes cl_N as an argument instead of numeric.

[ginac.git] / ginac / inifcns_nstdsums.cpp

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index e92e8eb22f9b95e4ec85c855bf9b46efd5693358..01a57569125cf1675a7dfbd522e21740c6cd3bc0 100644 (file)
--- a/ginac/inifcns_nstdsums.cpp
+++ b/ginac/inifcns_nstdsums.cpp
@@ -320,7 +320,7 @@ cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
 
 
 // forward declaration needed by function Li_projection and C below
-numeric S_num(int n, int p, const numeric& x);
+const cln::cl_N S_num(int n, int p, const cln::cl_N& x);
 
 
 // helper function for classical polylog Li
@@ -371,7 +371,7 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
                } else {
                        cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
                        for (int j=0; j<n-1; j++) {
-                               result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
+                               result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
                                                  * cln::expt(cln::log(x), j) / cln::factorial(j);
                        }
                        return result;
@@ -379,15 +379,14 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
        }
 }
 
-
 // helper function for classical polylog Li
-numeric Lin_numeric(int n, const numeric& x)
+const cln::cl_N Lin_numeric(const int n, const cln::cl_N& x)
 {
        if (n == 1) {
                // just a log
-               return -cln::log(1-x.to_cl_N());
+               return -cln::log(1-x);
        }
-       if (x.is_zero()) {
+       if (zerop(x)) {
                return 0;
        }
        if (x == 1) {
@@ -398,12 +397,11 @@ numeric Lin_numeric(int n, const numeric& x)
                // [Kol] (2.22)
                return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
        }
-       if (abs(x.real()) < 0.4 && abs(abs(x)-1) < 0.01) {
-               cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
-               cln::cl_N result = -cln::expt(cln::log(x_), n-1) * cln::log(1-x_) / cln::factorial(n-1);
+       if (cln::abs(realpart(x)) < 0.4 && cln::abs(cln::abs(x)-1) < 0.01) {
+               cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
                for (int j=0; j<n-1; j++) {
-                       result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x_).to_cl_N())
-                               * cln::expt(cln::log(x_), j) / cln::factorial(j);
+                       result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
+                               * cln::expt(cln::log(x), j) / cln::factorial(j);
                }
                return result;
        }
@@ -411,11 +409,11 @@ numeric Lin_numeric(int n, const numeric& x)
        // what is the desired float format?
        // first guess: default format
        cln::float_format_t prec = cln::default_float_format;
-       const cln::cl_N value = x.to_cl_N();
+       const cln::cl_N value = x;
        // second guess: the argument's format
-       if (!x.real().is_rational())
+       if (!instanceof(realpart(x), cln::cl_RA_ring))
                prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
-       else if (!x.imag().is_rational())
+       else if (!instanceof(imagpart(x), cln::cl_RA_ring))
                prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
        
        // [Kol] (5.15)
@@ -441,7 +439,7 @@ numeric Lin_numeric(int n, const numeric& x)
                cln::cl_N add;
                for (int j=0; j<n-1; j++) {
                        add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
-                                   * Lin_numeric(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
+                                   * Lin_numeric(n-j,1) * cln::expt(cln::log(-value),j) / cln::factorial(j);
                }
                result = result - add;
                return result;
@@ -488,11 +486,16 @@ cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl
                t0buf = t[0];
                q++;
                t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
+               for (int k=j-2; k>=0; k--) {
+                       t[k] = t[k] + t[k+1] * cln::expt(x[k], q+j-1-k) / cln::expt(cln::cl_I(q+j-1-k), s[k]);
+               }
+               q++;
+               t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
                for (int k=j-2; k>=0; k--) {
                        flag_accidental_zero = cln::zerop(t[k+1]);
                        t[k] = t[k] + t[k+1] * cln::expt(x[k], q+j-1-k) / cln::expt(cln::cl_I(q+j-1-k), s[k]);
                }
-       } while ( (t[0] != t0buf) || flag_accidental_zero );
+       } while ( (t[0] != t0buf) || cln::zerop(t[0]) || flag_accidental_zero );
 
        return t[0];
 }
@@ -991,7 +994,7 @@ ex G_numeric(const lst& x, const lst& s, const ex& y)
        if (x.op(x.nops()-1).is_zero()) {
                need_trafo = true;
        }
-       if (depth == 1 && !need_trafo) {
+       if (depth == 1 && x.nops() == 2 && !need_trafo) {
                return -Li(x.nops(), y / x.op(x.nops()-1)).evalf();
        }
        
@@ -1375,12 +1378,18 @@ static ex Li_evalf(const ex& m_, const ex& x_)
        // classical polylogs
        if (m_.info(info_flags::posint)) {
                if (x_.info(info_flags::numeric)) {
-                       return Lin_numeric(ex_to<numeric>(m_).to_int(), ex_to<numeric>(x_));
+                       int m__ = ex_to<numeric>(m_).to_int();
+                       const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
+                       const cln::cl_N result = Lin_numeric(m__, x__);
+                       return numeric(result);
                } else {
                        // try to numerically evaluate second argument
                        ex x_val = x_.evalf();
                        if (x_val.info(info_flags::numeric)) {
-                               return Lin_numeric(ex_to<numeric>(m_).to_int(), ex_to<numeric>(x_val));
+                               int m__ = ex_to<numeric>(m_).to_int();
+                               const cln::cl_N x__ = ex_to<numeric>(x_val).to_cl_N();
+                               const cln::cl_N result = Lin_numeric(m__, x__);
+                               return numeric(result);
                        }
                }
        }
@@ -1494,7 +1503,10 @@ static ex Li_eval(const ex& m_, const ex& x_)
                }
        }
        if (m_.info(info_flags::posint) && x_.info(info_flags::numeric) && !x_.info(info_flags::crational)) {
-               return Lin_numeric(ex_to<numeric>(m_).to_int(), ex_to<numeric>(x_));
+               int m__ = ex_to<numeric>(m_).to_int();
+               const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
+               const cln::cl_N result = Lin_numeric(m__, x__);
+               return numeric(result);
        }
 
        return Li(m_, x_).hold();
@@ -1715,10 +1727,10 @@ cln::cl_N C(int n, int p)
                        if (k == 0) {
                                if (n & 1) {
                                        if (j & 1) {
-                                               result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
+                                               result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
                                        }
                                        else {
-                                               result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
+                                               result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
                                        }
                                }
                        }
@@ -1726,23 +1738,23 @@ cln::cl_N C(int n, int p)
                                if (k & 1) {
                                        if (j & 1) {
                                                result = result + cln::factorial(n+k-1)
-                                                                 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+                                                                 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
                                                                  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
                                        }
                                        else {
                                                result = result - cln::factorial(n+k-1)
-                                                                 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+                                                                 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
                                                                  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
                                        }
                                }
                                else {
                                        if (j & 1) {
-                                               result = result - cln::factorial(n+k-1) * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+                                               result = result - cln::factorial(n+k-1) * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
                                                                  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
                                        }
                                        else {
                                                result = result + cln::factorial(n+k-1)
-                                                                 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
+                                                                 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
                                                                  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
                                        }
                                }
@@ -1804,10 +1816,20 @@ cln::cl_N b_k(int k)
 // helper function for S(n,p,x)
 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 {
+       static cln::float_format_t oldprec = cln::default_float_format;
+
        if (p==1) {
                return Li_projection(n+1, x, prec);
        }
-       
+
+       // precision has changed, we need to clear lookup table Yn
+       if ( oldprec != prec ) {
+               Yn.clear();
+               ynsize = 0;
+               ynlength = 100;
+               oldprec = prec;
+       }
+               
        // check if precalculated values are sufficient
        if (p > ynsize+1) {
                for (int i=ynsize; i<p-1; i++) {
@@ -1855,7 +1877,7 @@ cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format
                                res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
                                              * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
                        }
-                       result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
+                       result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
                }
 
                return result;
@@ -1866,7 +1888,7 @@ cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format
 
 
 // helper function for S(n,p,x)
-numeric S_num(int n, int p, const numeric& x)
+const cln::cl_N S_num(int n, int p, const cln::cl_N& x)
 {
        if (x == 1) {
                if (n == 1) {
@@ -1902,11 +1924,11 @@ numeric S_num(int n, int p, const numeric& x)
        // what is the desired float format?
        // first guess: default format
        cln::float_format_t prec = cln::default_float_format;
-       const cln::cl_N value = x.to_cl_N();
+       const cln::cl_N value = x;
        // second guess: the argument's format
-       if (!x.real().is_rational())
+       if (!instanceof(realpart(value), cln::cl_RA_ring))
                prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
-       else if (!x.imag().is_rational())
+       else if (!instanceof(imagpart(value), cln::cl_RA_ring))
                prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 
        // [Kol] (5.3)
@@ -1919,9 +1941,9 @@ numeric S_num(int n, int p, const numeric& x)
                        cln::cl_N res2;
                        for (int r=0; r<p; r++) {
                                res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
-                                             * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
+                                             * S_num(p-r,n-s,1-value) / cln::factorial(r);
                        }
-                       result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
+                       result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
                }
 
                return result;
@@ -1936,7 +1958,7 @@ numeric S_num(int n, int p, const numeric& x)
                        for (int r=0; r<=s; r++) {
                                result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
                                                  / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
-                                                 * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
+                                                 * S_num(n+s-r,p-s,cln::recip(value));
                        }
                }
                result = result * cln::expt(cln::cl_I(-1),n);
@@ -1972,12 +1994,18 @@ numeric S_num(int n, int p, const numeric& x)
 static ex S_evalf(const ex& n, const ex& p, const ex& x)
 {
        if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
+               const int n_ = ex_to<numeric>(n).to_int();
+               const int p_ = ex_to<numeric>(p).to_int();
                if (is_a<numeric>(x)) {
-                       return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
+                       const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
+                       const cln::cl_N result = S_num(n_, p_, x_);
+                       return numeric(result);
                } else {
                        ex x_val = x.evalf();
                        if (is_a<numeric>(x_val)) {
-                               return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x_val));
+                               const cln::cl_N x_val_ = ex_to<numeric>(x_val).to_cl_N();
+                               const cln::cl_N result = S_num(n_, p_, x_val_);
+                               return numeric(result);
                        }
                }
        }
@@ -2002,7 +2030,11 @@ static ex S_eval(const ex& n, const ex& p, const ex& x)
                        return Li(n+1, x);
                }
                if (x.info(info_flags::numeric) && (!x.info(info_flags::crational))) {
-                       return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
+                       int n_ = ex_to<numeric>(n).to_int();
+                       int p_ = ex_to<numeric>(p).to_int();
+                       const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
+                       const cln::cl_N result = S_num(n_, p_, x_);
+                       return numeric(result);
                }
        }
        if (n.is_zero()) {