G_numeric: fix numeric evaluation with trailing zeros and y != 1.

[ginac.git] / ginac / inifcns_nstdsums.cpp

diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp
index 4312377f10cf2d868972db41a56dac127973855d..62c6c62f2737cf6b87a34db73bc51060103c1ff5 100644 (file)
--- a/ginac/inifcns_nstdsums.cpp
+++ b/ginac/inifcns_nstdsums.cpp
@@ -47,7 +47,7 @@
  */
 
 /*
- *  GiNaC Copyright (C) 1999-2010 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2011 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
@@ -346,7 +346,11 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
                } else {
                        // choose the faster algorithm
                        if (cln::abs(cln::realpart(x)) > 0.75) {
-                               return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
+                               if ( x == 1 ) {
+                                       return cln::zeta(2);
+                               } else {
+                                       return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
+                               }
                        } else {
                                return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
                        }
@@ -368,7 +372,8 @@ cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& pr
                                return Lin_do_sum_Xn(n, x);
                        }
                } else {
-                       cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
+                       cln::cl_N result = 0;
+                       if ( x != 1 ) 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) - S_num(1, n-j-1, 1-x))
                                                  * cln::expt(cln::log(x), j) / cln::factorial(j);
@@ -1131,6 +1136,7 @@ G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
        // check for convergence and necessary accelerations
        bool need_trafo = false;
        bool need_hoelder = false;
+       bool have_trailing_zero = false;
        std::size_t depth = 0;
        for (std::size_t i = 0; i < x.size(); ++i) {
                if (!zerop(x[i])) {
@@ -1144,14 +1150,16 @@ G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
                                need_hoelder = true;
                }
        }
-       if (zerop(x[x.size() - 1]))
+       if (zerop(x[x.size() - 1])) {
+               have_trailing_zero = true;
                need_trafo = true;
+       }
 
        if (depth == 1 && x.size() == 2 && !need_trafo)
                return - Li_projection(2, y/x[1], cln::float_format(Digits));
        
        // do acceleration transformation (hoelder convolution [BBB])
-       if (need_hoelder)
+       if (need_hoelder && !have_trailing_zero)
                return G_do_hoelder(x, s, y);
        
        // convergence transformation
@@ -1196,9 +1204,14 @@ ex mLi_numeric(const lst& m, const lst& x)
                        s.push_back(1);
                }
                const cln::cl_N xi = ex_to<numeric>(*itx).to_cl_N();
-               newx.push_back(factor/xi);
                factor = factor/xi;
-               s.push_back(1);
+               newx.push_back(factor);
+               if ( !instanceof(factor, cln::cl_R_ring) && imagpart(factor) < 0 ) {
+                       s.push_back(-1);
+               }
+               else {
+                       s.push_back(1);
+               }
        }
        return numeric(cln::cl_N(1 & m.nops() ? - 1 : 1)*G_numeric(newx, s, cln::cl_N(1)));
 }
@@ -1350,12 +1363,16 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
                        all_zero = false;
                }
                if ( ex_to<numeric>(*itx).is_real() ) {
-                       if ( *its >= 0 ) {
+                       if ( ex_to<numeric>(*itx).is_positive() ) {
+                               if ( *its >= 0 ) {
+                                       sn.push_back(1);
+                               }
+                               else {
+                                       sn.push_back(-1);
+                               }
+                       } else {
                                sn.push_back(1);
                        }
-                       else {
-                               sn.push_back(-1);
-                       }
                }
                else {
                        if ( ex_to<numeric>(*itx).imag() > 0 ) {
@@ -1414,12 +1431,16 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
                        all_zero = false;
                }
                if ( ex_to<numeric>(*itx).is_real() ) {
-                       if ( *its >= 0 ) {
+                       if ( ex_to<numeric>(*itx).is_positive() ) {
+                               if ( *its >= 0 ) {
+                                       sn.push_back(1);
+                               }
+                               else {
+                                       sn.push_back(-1);
+                               }
+                       } else {
                                sn.push_back(1);
                        }
-                       else {
-                               sn.push_back(-1);
-                       }
                }
                else {
                        if ( ex_to<numeric>(*itx).imag() > 0 ) {
@@ -1559,7 +1580,17 @@ static ex Li_eval(const ex& m_, const ex& x_)
                                }
                        }
                        if (is_zeta) {
-                               return zeta(m_,x_);
+                               lst newx;
+                               for (lst::const_iterator itx = x.begin(); itx != x.end(); ++itx) {
+                                       GINAC_ASSERT((*itx == _ex1) || (*itx == _ex_1));
+                                       // XXX: 1 + 0.0*I is considered equal to 1. However
+                                       // the former is a not automatically converted
+                                       // to a real number. Do the conversion explicitly
+                                       // to avoid the "numeric::operator>(): complex inequality"
+                                       // exception (and similar problems).
+                                       newx.append(*itx != _ex_1 ? _ex1 : _ex_1);
+                               }
+                               return zeta(m_, newx);
                        }
                        if (is_H) {
                                ex prefactor;
@@ -2026,7 +2057,9 @@ const cln::cl_N S_num(int n, int p, const cln::cl_N& x)
                prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 
        // [Kol] (5.3)
-       if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95))) {
+       // the condition abs(1-value)>1 avoids an infinite recursion in the region abs(value)<=1 && abs(value)>0.95 && abs(1-value)<=1 && abs(1-value)>0.95
+       // we don't care here about abs(value)<1 && real(value)>0.5, this will be taken care of in S_projection
+       if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95) && (cln::abs(1-value) > 1) )) {
 
                cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
                                   * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
@@ -2067,6 +2100,16 @@ const cln::cl_N S_num(int n, int p, const cln::cl_N& x)
 
                return result;
        }
+
+       if ((cln::abs(value) > 0.95) && (cln::abs(value-9.53) < 9.47)) {
+               lst m;
+               m.append(n+1);
+               for (int s=0; s<p-1; s++)
+                       m.append(1);
+
+               ex res = H(m,numeric(value)).evalf();
+               return ex_to<numeric>(res).to_cl_N();
+       }
        else {
                return S_projection(n, p, value, prec);
        }
@@ -2464,7 +2507,12 @@ lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf)
        res.append(*itm);
        itm++;
        while (itx != x.end()) {
-               signum *= (*itx > 0) ? 1 : -1;
+               GINAC_ASSERT((*itx == _ex1) || (*itx == _ex_1));
+               // XXX: 1 + 0.0*I is considered equal to 1. However the former
+               // is not automatically converted to a real number.
+               // Do the conversion explicitly to avoid the
+               // "numeric::operator>(): complex inequality" exception.
+               signum *= (*itx != _ex_1) ? 1 : -1;
                pf *= signum;
                res.append((*itm) * signum);
                itm++;
@@ -2601,7 +2649,7 @@ ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
                ex addzeta = convert_H_to_zeta(newparameter);
                return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
        } else {
-               return e * (-H(lst(0),1/arg).hold());
+               return e * (-H(lst(ex(0)),1/arg).hold());
        }
 }
 
@@ -2632,7 +2680,7 @@ ex trafo_H_prepend_one(const ex& e, const ex& arg)
                newparameter.prepend(1);
                return e.subs(h == H(newparameter, h.op(1)).hold());
        } else {
-               return e * H(lst(1),1-arg).hold();
+               return e * H(lst(ex(1)),1-arg).hold();
        }
 }
 
@@ -2664,8 +2712,8 @@ ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
                ex addzeta = convert_H_to_zeta(newparameter);
                return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
        } else {
-               ex addzeta = convert_H_to_zeta(lst(-1));
-               return (e * (addzeta - H(lst(-1),1/arg).hold())).expand();
+               ex addzeta = convert_H_to_zeta(lst(ex(-1)));
+               return (e * (addzeta - H(lst(ex(-1)),1/arg).hold())).expand();
        }
 }
 
@@ -2696,7 +2744,7 @@ ex trafo_H_1mxt1px_prepend_minusone(const ex& e, const ex& arg)
                newparameter.prepend(-1);
                return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
        } else {
-               return (e * H(lst(-1),(1-arg)/(1+arg)).hold()).expand();
+               return (e * H(lst(ex(-1)),(1-arg)/(1+arg)).hold()).expand();
        }
 }
 
@@ -2727,7 +2775,7 @@ ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
                newparameter.prepend(1);
                return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
        } else {
-               return (e * H(lst(1),(1-arg)/(1+arg)).hold()).expand();
+               return (e * H(lst(ex(1)),(1-arg)/(1+arg)).hold()).expand();
        }
 }
 
@@ -2806,7 +2854,7 @@ struct map_trafo_H_1mx : public map_function
                                        // leading one
                                        map_trafo_H_1mx recursion;
                                        map_trafo_H_mult unify;
-                                       ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
+                                       ex res = H(lst(ex(1)), arg).hold() * H(newparameter, arg).hold();
                                        std::size_t firstzero = 0;
                                        while (parameter.op(firstzero) == 1) {
                                                firstzero++;
@@ -2870,7 +2918,7 @@ struct map_trafo_H_1overx : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(H(lst(-1),1/arg).hold() - H(lst(0),1/arg).hold(), parameter.nops())
+                                               return unify((pow(H(lst(ex(-1)),1/arg).hold() - H(lst(ex(0)),1/arg).hold(), parameter.nops())
                                                       / factorial(parameter.nops())).expand());
                                        }
                                } else {
@@ -2882,7 +2930,7 @@ struct map_trafo_H_1overx : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() + H_polesign, parameter.nops())
+                                               return unify((pow(H(lst(ex(1)),1/arg).hold() + H(lst(ex(0)),1/arg).hold() + H_polesign, parameter.nops())
                                                       / factorial(parameter.nops())).expand());
                                        }
                                }
@@ -2925,7 +2973,7 @@ struct map_trafo_H_1overx : public map_function
                                        // leading one
                                        map_trafo_H_1overx recursion;
                                        map_trafo_H_mult unify;
-                                       ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
+                                       ex res = H(lst(ex(1)), arg).hold() * H(newparameter, arg).hold();
                                        std::size_t firstzero = 0;
                                        while (parameter.op(firstzero) == 1) {
                                                firstzero++;
@@ -2981,7 +3029,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(-H(lst(1),(1-arg)/(1+arg)).hold() - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
+                                               return unify((pow(-H(lst(ex(1)),(1-arg)/(1+arg)).hold() - H(lst(ex(-1)),(1-arg)/(1+arg)).hold(), parameter.nops())
                                                       / factorial(parameter.nops())).expand());
                                        }
                                } else if (parameter.op(0) == -1) {
@@ -2993,7 +3041,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(log(2) - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
+                                               return unify((pow(log(2) - H(lst(ex(-1)),(1-arg)/(1+arg)).hold(), parameter.nops())
                                                       / factorial(parameter.nops())).expand());
                                        }
                                } else {
@@ -3005,7 +3053,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                        }
                                        if (allthesame) {
                                                map_trafo_H_mult unify;
-                                               return unify((pow(-log(2) - H(lst(0),(1-arg)/(1+arg)).hold() + H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
+                                               return unify((pow(-log(2) - H(lst(ex(0)),(1-arg)/(1+arg)).hold() + H(lst(ex(-1)),(1-arg)/(1+arg)).hold(), parameter.nops())
                                                       / factorial(parameter.nops())).expand());
                                        }
                                }
@@ -3048,7 +3096,7 @@ struct map_trafo_H_1mxt1px : public map_function
                                        // leading one
                                        map_trafo_H_1mxt1px recursion;
                                        map_trafo_H_mult unify;
-                                       ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
+                                       ex res = H(lst(ex(1)), arg).hold() * H(newparameter, arg).hold();
                                        std::size_t firstzero = 0;
                                        while (parameter.op(firstzero) == 1) {
                                                firstzero++;
@@ -3212,7 +3260,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                // x -> 1/x
                if (cln::abs(x) >= 2.0) {
                        map_trafo_H_1overx trafo;
-                       res *= trafo(H(m, xtemp));
+                       res *= trafo(H(m, xtemp).hold());
                        if (cln::imagpart(x) <= 0) {
                                res = res.subs(H_polesign == -I*Pi);
                        } else {
@@ -3227,7 +3275,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                if (cln::abs(x-9.53) <= 9.47) {
                        // x -> (1-x)/(1+x)
                        map_trafo_H_1mxt1px trafo;
-                       res *= trafo(H(m, xtemp));
+                       res *= trafo(H(m, xtemp).hold());
                } else {
                        // x -> 1-x
                        if (has_minus_one) {
@@ -3235,7 +3283,7 @@ static ex H_evalf(const ex& x1, const ex& x2)
                                return filter(H(m, numeric(x)).hold()).evalf();
                        }
                        map_trafo_H_1mx trafo;
-                       res *= trafo(H(m, xtemp));
+                       res *= trafo(H(m, xtemp).hold());
                }
 
                return res.subs(xtemp == numeric(x)).evalf();