X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Finifcns_nstdsums.cpp;h=a42c1bd3245895a697ca46108fe206afef175542;hp=4312377f10cf2d868972db41a56dac127973855d;hb=cda7fd319d75f95fbdac21213a10a04de724d228;hpb=c28e61da33905ddc69abf893aaffec98aa30d053 diff --git a/ginac/inifcns_nstdsums.cpp b/ginac/inifcns_nstdsums.cpp index 4312377f..a42c1bd3 100644 --- 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 0) { ex result; Gparameter new_a(a.begin(), a.end()-1); - result += G_eval1(0, scale, gsyms) * G_transform(pendint, new_a, scale, gsyms); + result += G_eval1(0, scale, gsyms) * G_transform(pendint, new_a, scale, gsyms, flag_trailing_zeros_only); for (Gparameter::const_iterator it = a.begin(); it != firstzero; ++it) { Gparameter new_a(a.begin(), it); new_a.push_back(0); new_a.insert(new_a.end(), it, a.end()-1); - result -= G_transform(pendint, new_a, scale, gsyms); + result -= G_transform(pendint, new_a, scale, gsyms, flag_trailing_zeros_only); } return result / trailing_zeros; } - // convergence case - if (convergent) { + // convergence case or flag_trailing_zeros_only + if (convergent || flag_trailing_zeros_only) { if (pendint.size() > 0) { return G_eval(convert_pending_integrals_G(pendint), pendint.front(), gsyms)* @@ -881,10 +886,10 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, Gparameter a1(a.begin(),min_it+1); Gparameter a2(min_it+1,a.end()); - ex result = G_transform(pendint, a2, scale, gsyms)* - G_transform(empty, a1, scale, gsyms); + ex result = G_transform(pendint, a2, scale, gsyms, flag_trailing_zeros_only)* + G_transform(empty, a1, scale, gsyms, flag_trailing_zeros_only); - result -= shuffle_G(empty, a1, a2, pendint, a, scale, gsyms); + result -= shuffle_G(empty, a1, a2, pendint, a, scale, gsyms, flag_trailing_zeros_only); return result; } @@ -895,7 +900,7 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, Gparameter new_pendint = prepare_pending_integrals(pendint, a[min_it_pos]); Gparameter new_a = a; new_a[min_it_pos] = 0; - ex result = G_transform(empty, new_a, scale, gsyms); + ex result = G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); if (pendint.size() > 0) { result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front(), gsyms); @@ -909,31 +914,31 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, new_pendint.push_back(*changeit); result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front(), gsyms)* - G_transform(empty, new_a, scale, gsyms); + G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); int buffer = *changeit; *changeit = *min_it; - result += G_transform(new_pendint, new_a, scale, gsyms); + result += G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only); *changeit = buffer; new_pendint.pop_back(); --changeit; new_pendint.push_back(*changeit); result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front(), gsyms)* - G_transform(empty, new_a, scale, gsyms); + G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); *changeit = *min_it; - result -= G_transform(new_pendint, new_a, scale, gsyms); + result -= G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only); } else { // smallest at the front new_pendint.push_back(scale); result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front(), gsyms)* - G_transform(empty, new_a, scale, gsyms); + G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); new_pendint.back() = *changeit; result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front(), gsyms)* - G_transform(empty, new_a, scale, gsyms); + G_transform(empty, new_a, scale, gsyms, flag_trailing_zeros_only); *changeit = *min_it; - result += G_transform(new_pendint, new_a, scale, gsyms); + result += G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only); } return result; } @@ -943,27 +948,27 @@ ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale, // for the one that is equal to a_old ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2, const Gparameter& pendint, const Gparameter& a_old, int scale, - const exvector& gsyms) + const exvector& gsyms, bool flag_trailing_zeros_only) { if (a1.size()==0 && a2.size()==0) { // veto the one configuration we don't want if ( a0 == a_old ) return 0; - return G_transform(pendint, a0, scale, gsyms); + return G_transform(pendint, a0, scale, gsyms, flag_trailing_zeros_only); } if (a2.size()==0) { Gparameter empty; Gparameter aa0 = a0; aa0.insert(aa0.end(),a1.begin(),a1.end()); - return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms); + return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms, flag_trailing_zeros_only); } if (a1.size()==0) { Gparameter empty; Gparameter aa0 = a0; aa0.insert(aa0.end(),a2.begin(),a2.end()); - return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms); + return shuffle_G(aa0, empty, empty, pendint, a_old, scale, gsyms, flag_trailing_zeros_only); } Gparameter a1_removed(a1.begin()+1,a1.end()); @@ -975,8 +980,8 @@ ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2 a01.push_back( a1[0] ); a02.push_back( a2[0] ); - return shuffle_G(a01, a1_removed, a2, pendint, a_old, scale, gsyms) - + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms); + return shuffle_G(a01, a1_removed, a2, pendint, a_old, scale, gsyms, flag_trailing_zeros_only) + + shuffle_G(a02, a1, a2_removed, pendint, a_old, scale, gsyms, flag_trailing_zeros_only); } // handles the transformations and the numerical evaluation of G @@ -1015,9 +1020,8 @@ G_do_hoelder(std::vector x, /* yes, it's passed by value */ std::vector 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 && realpart(x[j-1]) <= 2) { - qlsts.push_back(s[j-1]); + if (instanceof(x[j-1], cln::cl_R_ring) && realpart(x[j-1]) > 1) { + qlsts.push_back(1); } else { qlsts.push_back(-s[j-1]); } @@ -1039,24 +1043,43 @@ G_do_hoelder(std::vector x, /* yes, it's passed by value */ return result; } +class less_object_for_cl_N +{ +public: + bool operator() (const cln::cl_N & a, const cln::cl_N & b) const + { + // absolute value? + if (abs(a) != abs(b)) + return (abs(a) < abs(b)) ? true : false; + + // complex phase? + if (phase(a) != phase(b)) + return (phase(a) < phase(b)) ? true : false; + + // equal, therefore "less" is not true + return false; + } +}; + + // convergence transformation, used for numerical evaluation of G function. // the parameter x, s and y must only contain numerics static cln::cl_N G_do_trafo(const std::vector& x, const std::vector& s, - const cln::cl_N& y) + const cln::cl_N& y, bool flag_trailing_zeros_only) { // sort (|x|<->position) to determine indices - typedef std::multimap sortmap_t; + typedef std::multimap sortmap_t; sortmap_t sortmap; std::size_t size = 0; for (std::size_t i = 0; i < x.size(); ++i) { if (!zerop(x[i])) { - sortmap.insert(std::make_pair(abs(x[i]), i)); + sortmap.insert(std::make_pair(x[i], i)); ++size; } } // include upper limit (scale) - sortmap.insert(std::make_pair(abs(y), x.size())); + sortmap.insert(std::make_pair(y, x.size())); // generate missing dummy-symbols int i = 1; @@ -1111,7 +1134,7 @@ G_do_trafo(const std::vector& x, const std::vector& s, // do transformation Gparameter pendint; - ex result = G_transform(pendint, a, scale, gsyms); + ex result = G_transform(pendint, a, scale, gsyms, flag_trailing_zeros_only); // replace dummy symbols with their values result = result.eval().expand(); result = result.subs(subslst).evalf(); @@ -1131,6 +1154,7 @@ G_numeric(const std::vector& x, const std::vector& 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,19 +1168,21 @@ G_numeric(const std::vector& x, const std::vector& s, need_hoelder = true; } } - if (zerop(x[x.size() - 1])) + if (zerop(x.back())) { + 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 if (need_trafo) - return G_do_trafo(x, s, y); + return G_do_trafo(x, s, y, have_trailing_zero); // do summation std::vector newx; @@ -1196,9 +1222,14 @@ ex mLi_numeric(const lst& m, const lst& x) s.push_back(1); } const cln::cl_N xi = ex_to(*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 +1381,16 @@ static ex G3_evalf(const ex& x_, const ex& s_, const ex& y) all_zero = false; } if ( ex_to(*itx).is_real() ) { - if ( *its >= 0 ) { + if ( ex_to(*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(*itx).imag() > 0 ) { @@ -1414,12 +1449,16 @@ static ex G3_eval(const ex& x_, const ex& s_, const ex& y) all_zero = false; } if ( ex_to(*itx).is_real() ) { - if ( *its >= 0 ) { + if ( ex_to(*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(*itx).imag() > 0 ) { @@ -1559,7 +1598,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 +2075,9 @@ const cln::cl_N S_num(int n, int p, const cln::cl_N& x) prec = cln::float_format(cln::the(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 +2118,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(res).to_cl_N(); + } else { return S_projection(n, p, value, prec); } @@ -2464,7 +2525,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 +2667,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 +2698,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 +2730,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 +2762,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 +2793,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 +2872,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 +2936,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 +2948,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 +2991,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 +3047,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 +3059,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 +3071,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 +3114,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 +3278,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 +3293,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 +3301,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();