1 /** @file inifcns_nstdsums.cpp
3 * Implementation of some special functions that have a representation as nested sums.
6 * classical polylogarithm Li(n,x)
7 * multiple polylogarithm Li(lst(n_1,...,n_k),lst(x_1,...,x_k))
8 * nielsen's generalized polylogarithm S(n,p,x)
9 * harmonic polylogarithm H(n,x) or H(lst(n_1,...,n_k),x)
10 * multiple zeta value zeta(n) or zeta(lst(n_1,...,n_k))
14 * - All formulae used can be looked up in the following publications:
15 * [Kol] Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
16 * [Cra] Fast Evaluation of Multiple Zeta Sums, R.E.Crandall, Math.Comp. 67 (1998), pp. 1163-1172.
17 * [ReV] Harmonic Polylogarithms, E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
19 * - The order of parameters and arguments of H, Li and zeta is defined according to their order in the
20 * nested sums representation.
22 * - Except for the multiple polylogarithm all functions can be nummerically evaluated with arguments in
23 * the whole complex plane. Multiple polylogarithms evaluate only if each argument x_i is smaller than
24 * one. The parameters for every function (n, p or n_i) must be positive integers.
26 * - The calculation of classical polylogarithms is speed up by using Bernoulli numbers and
27 * look-up tables. S uses look-up tables as well. The zeta function applies the algorithm in
30 * - The functions have no series expansion as nested sums. To do it, you have to convert these functions
31 * into the appropriate objects from the nestedsums library, do the expansion and convert the
34 * - Numerical testing of this implementation has been performed by doing a comparison of results
35 * between this software and the commercial M.......... 4.1. Multiple zeta values have been checked
36 * by means of evaluations into simple zeta values. Harmonic polylogarithms have been checked by
37 * comparison to S(n,p,x) for corresponding parameter combinations and by continuity checks
38 * around |x|=1 along with comparisons to corresponding zeta functions.
43 * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
45 * This program is free software; you can redistribute it and/or modify
46 * it under the terms of the GNU General Public License as published by
47 * the Free Software Foundation; either version 2 of the License, or
48 * (at your option) any later version.
50 * This program is distributed in the hope that it will be useful,
51 * but WITHOUT ANY WARRANTY; without even the implied warranty of
52 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
53 * GNU General Public License for more details.
55 * You should have received a copy of the GNU General Public License
56 * along with this program; if not, write to the Free Software
57 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
71 #include "operators.h"
74 #include "relational.h"
83 //////////////////////////////////////////////////////////////////////
85 // Classical polylogarithm Li
89 //////////////////////////////////////////////////////////////////////
92 // anonymous namespace for helper functions
96 // lookup table for factors built from Bernoulli numbers
98 std::vector<std::vector<cln::cl_N> > Xn;
102 // This function calculates the X_n. The X_n are needed for speed up of classical polylogarithms.
103 // With these numbers the polylogs can be calculated as follows:
104 // Li_p (x) = \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with u = -log(1-x)
105 // X_0(n) = B_n (Bernoulli numbers)
106 // X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
107 // The calculation of Xn depends on X0 and X{n-1}.
108 // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
109 // This results in a slightly more complicated algorithm for the X_n.
110 // The first index in Xn corresponds to the index of the polylog minus 2.
111 // The second index in Xn corresponds to the index from the actual sum.
114 // rule of thumb. needs to be improved. TODO
115 const int initsize = Digits * 3 / 2;
118 // calculate X_2 and higher (corresponding to Li_4 and higher)
119 std::vector<cln::cl_N> buf(initsize);
120 std::vector<cln::cl_N>::iterator it = buf.begin();
122 *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
124 for (int i=2; i<=initsize; i++) {
126 result = 0; // k == 0
128 result = Xn[0][i/2-1]; // k == 0
130 for (int k=1; k<i-1; k++) {
131 if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
132 result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
135 result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
136 result = result + Xn[n-1][i-1] / (i+1); // k == i
143 // special case to handle the X_0 correct
144 std::vector<cln::cl_N> buf(initsize);
145 std::vector<cln::cl_N>::iterator it = buf.begin();
147 *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
149 *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
151 for (int i=3; i<=initsize; i++) {
153 result = -Xn[0][(i-3)/2]/2;
154 *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
157 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
158 for (int k=1; k<i/2; k++) {
159 result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
168 std::vector<cln::cl_N> buf(initsize/2);
169 std::vector<cln::cl_N>::iterator it = buf.begin();
170 for (int i=1; i<=initsize/2; i++) {
171 *it = bernoulli(i*2).to_cl_N();
181 // calculates Li(2,x) without Xn
182 cln::cl_N Li2_do_sum(const cln::cl_N& x)
186 cln::cl_N num = x * cln::cl_float(1, cln::float_format(Digits));
187 cln::cl_I den = 1; // n^2 = 1
192 den = den + i; // n^2 = 4, 9, 16, ...
194 res = res + num / den;
195 } while (res != resbuf);
200 // calculates Li(2,x) with Xn
201 cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x)
203 std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
204 cln::cl_N u = -cln::log(1-x);
205 cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
206 cln::cl_N res = u - u*u/4;
211 factor = factor * u*u / (2*i * (2*i+1));
212 res = res + (*it) * factor;
213 it++; // should we check it? or rely on initsize? ...
215 } while (res != resbuf);
220 // calculates Li(n,x), n>2 without Xn
221 cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
223 cln::cl_N factor = x * cln::cl_float(1, cln::float_format(Digits));
230 res = res + factor / cln::expt(cln::cl_I(i),n);
232 } while (res != resbuf);
237 // calculates Li(n,x), n>2 with Xn
238 cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
240 std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
241 cln::cl_N u = -cln::log(1-x);
242 cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
248 factor = factor * u / i;
249 res = res + (*it) * factor;
250 it++; // should we check it? or rely on initsize? ...
252 } while (res != resbuf);
257 // forward declaration needed by function Li_projection and C below
258 numeric S_num(int n, int p, const numeric& x);
261 // helper function for classical polylog Li
262 cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
264 // treat n=2 as special case
266 // check if precalculated X0 exists
271 if (cln::realpart(x) < 0.5) {
272 // choose the faster algorithm
273 // the switching point was empirically determined. the optimal point
274 // depends on hardware, Digits, ... so an approx value is okay.
275 // it solves also the problem with precision due to the u=-log(1-x) transformation
276 if (cln::abs(cln::realpart(x)) < 0.25) {
278 return Li2_do_sum(x);
280 return Li2_do_sum_Xn(x);
283 // choose the faster algorithm
284 if (cln::abs(cln::realpart(x)) > 0.75) {
285 return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
287 return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
291 // check if precalculated Xn exist
293 for (int i=xnsize; i<n-1; i++) {
298 if (cln::realpart(x) < 0.5) {
299 // choose the faster algorithm
300 // with n>=12 the "normal" summation always wins against the method with Xn
301 if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
302 return Lin_do_sum(n, x);
304 return Lin_do_sum_Xn(n, x);
307 cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
308 for (int j=0; j<n-1; j++) {
309 result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
310 * cln::expt(cln::log(x), j) / cln::factorial(j);
318 // helper function for classical polylog Li
319 numeric Li_num(int n, const numeric& x)
323 return -cln::log(1-x.to_cl_N());
334 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
337 // what is the desired float format?
338 // first guess: default format
339 cln::float_format_t prec = cln::default_float_format;
340 const cln::cl_N value = x.to_cl_N();
341 // second guess: the argument's format
342 if (!x.real().is_rational())
343 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
344 else if (!x.imag().is_rational())
345 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
348 if (cln::abs(value) > 1) {
349 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
350 // check if argument is complex. if it is real, the new polylog has to be conjugated.
351 if (cln::zerop(cln::imagpart(value))) {
353 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
356 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
361 result = result + Li_projection(n, cln::recip(value), prec);
364 result = result - Li_projection(n, cln::recip(value), prec);
368 for (int j=0; j<n-1; j++) {
369 add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
370 * Li_num(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
372 result = result - add;
376 return Li_projection(n, value, prec);
381 } // end of anonymous namespace
384 //////////////////////////////////////////////////////////////////////
386 // Multiple polylogarithm Li
390 //////////////////////////////////////////////////////////////////////
393 // anonymous namespace for helper function
397 cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
399 const int j = s.size();
401 std::vector<cln::cl_N> t(j);
402 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
409 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
410 for (int k=j-2; k>=0; k--) {
411 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]);
413 } while (t[0] != t0buf);
419 } // end of anonymous namespace
422 //////////////////////////////////////////////////////////////////////
424 // Classical polylogarithm and multiple polylogarithm Li
428 //////////////////////////////////////////////////////////////////////
431 static ex Li_eval(const ex& x1, const ex& x2)
437 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational)))
438 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
440 for (int i=0; i<x2.nops(); i++) {
441 if (!is_a<numeric>(x2.op(i))) {
442 return Li(x1,x2).hold();
445 return Li(x1,x2).evalf();
447 return Li(x1,x2).hold();
452 static ex Li_evalf(const ex& x1, const ex& x2)
454 // classical polylogs
455 if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
456 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
459 else if (is_a<lst>(x1) && is_a<lst>(x2)) {
460 for (int i=0; i<x1.nops(); i++) {
461 if (!x1.op(i).info(info_flags::posint)) {
462 return Li(x1,x2).hold();
464 if (!is_a<numeric>(x2.op(i))) {
465 return Li(x1,x2).hold();
467 if (abs(x2.op(i)) >= 1) {
468 return Li(x1,x2).hold();
473 std::vector<cln::cl_N> x;
474 for (int i=0; i<ex_to<numeric>(x1.nops()).to_int(); i++) {
475 m.push_back(ex_to<numeric>(x1.op(i)).to_int());
476 x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
479 return numeric(multipleLi_do_sum(m, x));
482 return Li(x1,x2).hold();
486 static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
489 seq.push_back(expair(Li(x1,x2), 0));
490 return pseries(rel,seq);
494 static ex Li_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
496 GINAC_ASSERT(deriv_param < 2);
497 if (deriv_param == 0) {
501 return Li(x1-1, x2) / x2;
508 REGISTER_FUNCTION(Li,
510 evalf_func(Li_evalf).
511 do_not_evalf_params().
512 series_func(Li_series).
513 derivative_func(Li_deriv));
516 //////////////////////////////////////////////////////////////////////
518 // Nielsen's generalized polylogarithm S
522 //////////////////////////////////////////////////////////////////////
525 // anonymous namespace for helper functions
529 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
531 std::vector<std::vector<cln::cl_N> > Yn;
532 int ynsize = 0; // number of Yn[]
533 int ynlength = 100; // initial length of all Yn[i]
536 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
537 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
538 // representing S_{n,p}(x).
539 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
541 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
542 // representing S_{n,p}(x).
543 // The calculation of Y_n uses the values from Y_{n-1}.
544 void fill_Yn(int n, const cln::float_format_t& prec)
546 const int initsize = ynlength;
547 //const int initsize = initsize_Yn;
548 cln::cl_N one = cln::cl_float(1, prec);
551 std::vector<cln::cl_N> buf(initsize);
552 std::vector<cln::cl_N>::iterator it = buf.begin();
553 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
554 *it = (*itprev) / cln::cl_N(n+1) * one;
557 // sums with an index smaller than the depth are zero and need not to be calculated.
558 // calculation starts with depth, which is n+2)
559 for (int i=n+2; i<=initsize+n; i++) {
560 *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
566 std::vector<cln::cl_N> buf(initsize);
567 std::vector<cln::cl_N>::iterator it = buf.begin();
570 for (int i=2; i<=initsize; i++) {
571 *it = *(it-1) + 1 / cln::cl_N(i) * one;
580 // make Yn longer ...
581 void make_Yn_longer(int newsize, const cln::float_format_t& prec)
584 cln::cl_N one = cln::cl_float(1, prec);
586 Yn[0].resize(newsize);
587 std::vector<cln::cl_N>::iterator it = Yn[0].begin();
589 for (int i=ynlength+1; i<=newsize; i++) {
590 *it = *(it-1) + 1 / cln::cl_N(i) * one;
594 for (int n=1; n<ynsize; n++) {
595 Yn[n].resize(newsize);
596 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
597 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
600 for (int i=ynlength+n+1; i<=newsize+n; i++) {
601 *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
611 // helper function for S(n,p,x)
613 cln::cl_N C(int n, int p)
617 for (int k=0; k<p; k++) {
618 for (int j=0; j<=(n+k-1)/2; j++) {
622 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
625 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
632 result = result + cln::factorial(n+k-1)
633 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
634 / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
637 result = result - cln::factorial(n+k-1)
638 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
639 / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
644 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()
645 / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
648 result = result + cln::factorial(n+k-1)
649 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
650 / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
658 if (((np)/2+n) & 1) {
659 result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
662 result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
670 // helper function for S(n,p,x)
671 // [Kol] remark to (9.1)
681 for (int m=2; m<=k; m++) {
682 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
689 // helper function for S(n,p,x)
690 // [Kol] remark to (9.1)
700 for (int m=2; m<=k; m++) {
701 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
708 // helper function for S(n,p,x)
709 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
712 return Li_projection(n+1, x, prec);
715 // check if precalculated values are sufficient
717 for (int i=ynsize; i<p-1; i++) {
722 // should be done otherwise
723 cln::cl_N xf = x * cln::cl_float(1, prec);
727 cln::cl_N factor = cln::expt(xf, p);
731 if (i-p >= ynlength) {
733 make_Yn_longer(ynlength*2, prec);
735 res = res + factor / cln::expt(cln::cl_I(i),n+1) * Yn[p-2][i-p]; // should we check it? or rely on magic number? ...
736 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
737 factor = factor * xf;
739 } while (res != resbuf);
745 // helper function for S(n,p,x)
746 cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
749 if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
751 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
752 * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
754 for (int s=0; s<n; s++) {
756 for (int r=0; r<p; r++) {
757 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
758 * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
760 result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
766 return S_do_sum(n, p, x, prec);
770 // helper function for S(n,p,x)
771 numeric S_num(int n, int p, const numeric& x)
775 // [Kol] (2.22) with (2.21)
776 return cln::zeta(p+1);
781 return cln::zeta(n+1);
786 for (int nu=0; nu<n; nu++) {
787 for (int rho=0; rho<=p; rho++) {
788 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
789 * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
792 result = result * cln::expt(cln::cl_I(-1),n+p-1);
799 return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
801 // throw std::runtime_error("don't know how to evaluate this function!");
804 // what is the desired float format?
805 // first guess: default format
806 cln::float_format_t prec = cln::default_float_format;
807 const cln::cl_N value = x.to_cl_N();
808 // second guess: the argument's format
809 if (!x.real().is_rational())
810 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
811 else if (!x.imag().is_rational())
812 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
816 if (cln::realpart(value) < -0.5) {
818 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
819 * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
821 for (int s=0; s<n; s++) {
823 for (int r=0; r<p; r++) {
824 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
825 * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
827 result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
834 if (cln::abs(value) > 1) {
838 for (int s=0; s<p; s++) {
839 for (int r=0; r<=s; r++) {
840 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
841 / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
842 * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
845 result = result * cln::expt(cln::cl_I(-1),n);
848 for (int r=0; r<n; r++) {
849 res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
851 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
853 result = result + cln::expt(cln::cl_I(-1),p) * res2;
858 return S_projection(n, p, value, prec);
863 } // end of anonymous namespace
866 //////////////////////////////////////////////////////////////////////
868 // Nielsen's generalized polylogarithm S
872 //////////////////////////////////////////////////////////////////////
875 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
880 if (x3.info(info_flags::numeric) && (!x3.info(info_flags::crational)) &&
881 x1.info(info_flags::posint) && x2.info(info_flags::posint)) {
882 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
884 return S(x1,x2,x3).hold();
888 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
890 if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
891 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
893 return S(x1,x2,x3).hold();
897 static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options)
900 seq.push_back(expair(S(x1,x2,x3), 0));
901 return pseries(rel,seq);
905 static ex S_deriv(const ex& x1, const ex& x2, const ex& x3, unsigned deriv_param)
907 GINAC_ASSERT(deriv_param < 3);
908 if (deriv_param < 2) {
912 return S(x1-1, x2, x3) / x3;
914 return S(x1, x2-1, x3) / (1-x3);
922 do_not_evalf_params().
923 series_func(S_series).
924 derivative_func(S_deriv));
927 //////////////////////////////////////////////////////////////////////
929 // Harmonic polylogarithm H
933 //////////////////////////////////////////////////////////////////////
936 // anonymous namespace for helper functions
940 // forward declaration
941 ex convert_from_RV(const lst& parameterlst, const ex& arg);
944 // multiplies an one-dimensional H with another H
946 ex trafo_H_mult(const ex& h1, const ex& h2)
951 ex h1nops = h1.op(0).nops();
952 ex h2nops = h2.op(0).nops();
954 hshort = h2.op(0).op(0);
955 hlong = ex_to<lst>(h1.op(0));
957 hshort = h1.op(0).op(0);
959 hlong = ex_to<lst>(h2.op(0));
961 hlong = h2.op(0).op(0);
964 for (int i=0; i<=hlong.nops(); i++) {
968 newparameter.append(hlong[j]);
970 newparameter.append(hshort);
971 for (; j<hlong.nops(); j++) {
972 newparameter.append(hlong[j]);
974 res += H(newparameter, h1.op(1)).hold();
980 // applies trafo_H_mult recursively on expressions
981 struct map_trafo_H_mult : public map_function
983 ex operator()(const ex& e)
994 for (int pos=0; pos<e.nops(); pos++) {
995 if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
996 std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
998 for (ex i=0; i<e.op(pos).op(1); i++) {
999 Hlst.append(e.op(pos).op(0));
1003 } else if (is_a<function>(e.op(pos))) {
1004 std::string name = ex_to<function>(e.op(pos)).get_name();
1006 if (e.op(pos).op(0).nops() > 1) {
1009 Hlst.append(e.op(pos));
1014 result *= e.op(pos);
1017 if (Hlst.nops() > 0) {
1018 firstH = Hlst[Hlst.nops()-1];
1025 if (Hlst.nops() > 0) {
1026 ex buffer = trafo_H_mult(firstH, Hlst.op(0));
1028 for (int i=1; i<Hlst.nops(); i++) {
1029 result *= Hlst.op(i);
1031 result = result.expand();
1032 map_trafo_H_mult recursion;
1033 return recursion(result);
1044 // do integration [ReV] (49)
1045 // put parameter 1 in front of existing parameters
1046 ex trafo_H_prepend_one(const ex& e, const ex& arg)
1050 if (is_a<function>(e)) {
1051 name = ex_to<function>(e).get_name();
1056 for (int i=0; i<e.nops(); i++) {
1057 if (is_a<function>(e.op(i))) {
1058 std::string name = ex_to<function>(e.op(i)).get_name();
1066 lst newparameter = ex_to<lst>(h.op(0));
1067 newparameter.prepend(1);
1068 return e.subs(h == H(newparameter, h.op(1)).hold());
1070 return e * H(lst(1),1-arg).hold();
1075 // do integration [ReV] (55)
1076 // put parameter 0 in front of existing parameters
1077 ex trafo_H_prepend_zero(const ex& e, const ex& arg)
1081 if (is_a<function>(e)) {
1082 name = ex_to<function>(e).get_name();
1087 for (int i=0; i<e.nops(); i++) {
1088 if (is_a<function>(e.op(i))) {
1089 std::string name = ex_to<function>(e.op(i)).get_name();
1097 lst newparameter = ex_to<lst>(h.op(0));
1098 newparameter.prepend(0);
1099 ex addzeta = convert_from_RV(newparameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1100 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
1102 return e * (-H(lst(0),1/arg).hold());
1107 // do x -> 1-x transformation
1108 struct map_trafo_H_1mx : public map_function
1110 ex operator()(const ex& e)
1112 if (is_a<add>(e) || is_a<mul>(e)) {
1113 return e.map(*this);
1116 if (is_a<function>(e)) {
1117 std::string name = ex_to<function>(e).get_name();
1120 lst parameter = ex_to<lst>(e.op(0));
1123 // if all parameters are either zero or one return the transformed function
1124 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1126 for (int i=parameter.nops(); i>0; i--) {
1127 newparameter.append(0);
1129 return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1130 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1132 for (int i=parameter.nops(); i>0; i--) {
1133 newparameter.append(1);
1135 return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1138 lst newparameter = parameter;
1139 newparameter.remove_first();
1141 if (parameter.op(0) == 0) {
1144 ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1145 map_trafo_H_1mx recursion;
1146 ex buffer = recursion(H(newparameter, arg).hold());
1147 if (is_a<add>(buffer)) {
1148 for (int i=0; i<buffer.nops(); i++) {
1149 res -= trafo_H_prepend_one(buffer.op(i), arg);
1152 res -= trafo_H_prepend_one(buffer, arg);
1159 map_trafo_H_1mx recursion;
1160 map_trafo_H_mult unify;
1163 while (parameter.op(firstzero) == 1) {
1166 for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1170 newparameter.append(parameter[j+1]);
1172 newparameter.append(1);
1173 for (; j<parameter.nops()-1; j++) {
1174 newparameter.append(parameter[j+1]);
1176 res -= H(newparameter, arg).hold();
1178 return (unify((-H(lst(0), 1-arg).hold() * recursion(H(newparameter, arg).hold())).expand()) +
1179 recursion(res)) / firstzero;
1190 // do x -> 1/x transformation
1191 struct map_trafo_H_1overx : public map_function
1193 ex operator()(const ex& e)
1195 if (is_a<add>(e) || is_a<mul>(e)) {
1196 return e.map(*this);
1199 if (is_a<function>(e)) {
1200 std::string name = ex_to<function>(e).get_name();
1203 lst parameter = ex_to<lst>(e.op(0));
1206 // if all parameters are either zero or one return the transformed function
1207 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1208 map_trafo_H_mult unify;
1209 return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() - I*Pi, parameter.nops()) /
1210 factorial(parameter.nops())).expand());
1211 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1212 return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
1215 lst newparameter = parameter;
1216 newparameter.remove_first();
1218 if (parameter.op(0) == 0) {
1221 ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1222 map_trafo_H_1overx recursion;
1223 ex buffer = recursion(H(newparameter, arg).hold());
1224 if (is_a<add>(buffer)) {
1225 for (int i=0; i<buffer.nops(); i++) {
1226 res += trafo_H_prepend_zero(buffer.op(i), arg);
1229 res += trafo_H_prepend_zero(buffer, arg);
1236 map_trafo_H_1overx recursion;
1237 map_trafo_H_mult unify;
1238 ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
1240 while (parameter.op(firstzero) == 1) {
1243 for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1247 newparameter.append(parameter[j+1]);
1249 newparameter.append(1);
1250 for (; j<parameter.nops()-1; j++) {
1251 newparameter.append(parameter[j+1]);
1253 res -= H(newparameter, arg).hold();
1255 res = recursion(res).expand() / firstzero;
1267 // remove trailing zeros from H-parameters
1268 struct map_trafo_H_reduce_trailing_zeros : public map_function
1270 ex operator()(const ex& e)
1272 if (is_a<add>(e) || is_a<mul>(e)) {
1273 return e.map(*this);
1275 if (is_a<function>(e)) {
1276 std::string name = ex_to<function>(e).get_name();
1279 if (is_a<lst>(e.op(0))) {
1280 parameter = ex_to<lst>(e.op(0));
1282 parameter = lst(e.op(0));
1285 if (parameter.op(parameter.nops()-1) == 0) {
1288 if (parameter.nops() == 1) {
1293 lst::const_iterator it = parameter.begin();
1294 while ((it != parameter.end()) && (*it == 0)) {
1297 if (it == parameter.end()) {
1298 return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
1302 parameter.remove_last();
1303 int lastentry = parameter.nops();
1304 while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
1309 ex result = log(arg) * H(parameter,arg).hold();
1310 for (ex i=0; i<lastentry; i++) {
1312 result -= (parameter[i]-1) * H(parameter, arg).hold();
1316 if (lastentry < parameter.nops()) {
1317 result = result / (parameter.nops()-lastentry+1);
1318 return result.map(*this);
1330 // recursively call convert_from_RV on expression
1331 struct map_trafo_H_convert : public map_function
1333 ex operator()(const ex& e)
1335 if (is_a<add>(e) || is_a<mul>(e) || is_a<power>(e)) {
1336 return e.map(*this);
1338 if (is_a<function>(e)) {
1339 std::string name = ex_to<function>(e).get_name();
1341 lst parameter = ex_to<lst>(e.op(0));
1343 return convert_from_RV(parameter, arg);
1351 // translate notation from nested sums to Remiddi/Vermaseren
1352 lst convert_to_RV(const lst& o)
1355 for (lst::const_iterator it = o.begin(); it != o.end(); it++) {
1356 for (ex i=0; i<(*it)-1; i++) {
1365 // translate notation from Remiddi/Vermaseren to nested sums
1366 ex convert_from_RV(const lst& parameterlst, const ex& arg)
1368 lst newparameterlst;
1370 lst::const_iterator it = parameterlst.begin();
1372 while (it != parameterlst.end()) {
1376 newparameterlst.append((*it>0) ? count : -count);
1381 for (int i=1; i<count; i++) {
1382 newparameterlst.append(0);
1385 map_trafo_H_reduce_trailing_zeros filter;
1386 return filter(H(newparameterlst, arg).hold());
1390 // do the actual summation.
1391 cln::cl_N H_do_sum(const std::vector<int>& s, const cln::cl_N& x)
1393 const int j = s.size();
1395 std::vector<cln::cl_N> t(j);
1397 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1398 cln::cl_N factor = cln::expt(x, j) * one;
1404 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1405 for (int k=j-2; k>=1; k--) {
1406 t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1408 t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), s[0]);
1409 factor = factor * x;
1410 } while (t[0] != t0buf);
1416 } // end of anonymous namespace
1419 //////////////////////////////////////////////////////////////////////
1421 // Harmonic polylogarithm H
1425 //////////////////////////////////////////////////////////////////////
1428 static ex H_eval(const ex& x1, const ex& x2)
1436 if (x1.nops() == 1) {
1437 return Li(x1.op(0), x2);
1439 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational))) {
1440 return H(x1,x2).evalf();
1442 return H(x1,x2).hold();
1446 static ex H_evalf(const ex& x1, const ex& x2)
1448 if (is_a<lst>(x1) && is_a<numeric>(x2)) {
1449 for (int i=0; i<x1.nops(); i++) {
1450 if (!x1.op(i).info(info_flags::posint)) {
1451 return H(x1,x2).hold();
1454 if (x1.nops() < 1) {
1457 if (x1.nops() == 1) {
1458 return Li(x1.op(0), x2).evalf();
1460 cln::cl_N x = ex_to<numeric>(x2).to_cl_N();
1462 return zeta(x1).evalf();
1466 if (cln::abs(x) > 1) {
1467 symbol xtemp("xtemp");
1468 map_trafo_H_1overx trafo;
1469 ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1470 map_trafo_H_convert converter;
1471 res = converter(res);
1472 return res.subs(xtemp==x2).evalf();
1475 // since the x->1-x transformation produces a lot of terms, it is only
1476 // efficient for argument near one.
1477 if (cln::realpart(x) > 0.95) {
1478 symbol xtemp("xtemp");
1479 map_trafo_H_1mx trafo;
1480 ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1481 map_trafo_H_convert converter;
1482 res = converter(res);
1483 return res.subs(xtemp==x2).evalf();
1486 // no trafo -> do summation
1487 int count = x1.nops();
1488 std::vector<int> r(count);
1489 for (int i=0; i<count; i++) {
1490 r[i] = ex_to<numeric>(x1.op(i)).to_int();
1493 return numeric(H_do_sum(r,x));
1496 return H(x1,x2).hold();
1500 static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
1503 seq.push_back(expair(H(x1,x2), 0));
1504 return pseries(rel,seq);
1508 static ex H_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
1510 GINAC_ASSERT(deriv_param < 2);
1511 if (deriv_param == 0) {
1514 if (is_a<lst>(x1)) {
1515 lst newparameter = ex_to<lst>(x1);
1516 if (x1.op(0) == 1) {
1517 newparameter.remove_first();
1518 return 1/(1-x2) * H(newparameter, x2);
1521 return H(newparameter, x2).hold() / x2;
1527 return H(x1-1, x2).hold() / x2;
1533 REGISTER_FUNCTION(H,
1535 evalf_func(H_evalf).
1536 do_not_evalf_params().
1537 series_func(H_series).
1538 derivative_func(H_deriv));
1541 //////////////////////////////////////////////////////////////////////
1543 // Multiple zeta values zeta
1547 //////////////////////////////////////////////////////////////////////
1550 // anonymous namespace for helper functions
1554 // parameters and data for [Cra] algorithm
1555 const cln::cl_N lambda = cln::cl_N("319/320");
1558 std::vector<std::vector<cln::cl_N> > f_kj;
1559 std::vector<cln::cl_N> crB;
1560 std::vector<std::vector<cln::cl_N> > crG;
1561 std::vector<cln::cl_N> crX;
1564 void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
1566 const int size = a.size();
1567 for (int n=0; n<size; n++) {
1569 for (int m=0; m<=n; m++) {
1570 c[n] = c[n] + a[m]*b[n-m];
1577 void initcX(const std::vector<int>& s)
1579 const int k = s.size();
1585 for (int i=0; i<=L2; i++) {
1586 crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
1591 for (int m=0; m<k-1; m++) {
1592 std::vector<cln::cl_N> crGbuf;
1595 for (int i=0; i<=L2; i++) {
1596 crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
1598 crG.push_back(crGbuf);
1603 for (int m=0; m<k-1; m++) {
1604 std::vector<cln::cl_N> Xbuf;
1605 for (int i=0; i<=L2; i++) {
1606 Xbuf.push_back(crX[i] * crG[m][i]);
1608 halfcyclic_convolute(Xbuf, crB, crX);
1614 cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
1616 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1617 cln::cl_N factor = cln::expt(lambda, Sqk);
1618 cln::cl_N res = factor / Sqk * crX[0] * one;
1623 factor = factor * lambda;
1625 res = res + crX[N] * factor / (N+Sqk);
1626 } while ((res != resbuf) || cln::zerop(crX[N]));
1632 void calc_f(int maxr)
1637 cln::cl_N t0, t1, t2, t3, t4;
1639 std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
1640 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1642 t0 = cln::exp(-lambda);
1644 for (k=1; k<=L1; k++) {
1647 for (j=1; j<=maxr; j++) {
1650 for (i=2; i<=j; i++) {
1654 (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one);
1662 cln::cl_N crandall_Z(const std::vector<int>& s)
1664 const int j = s.size();
1673 t0 = t0 + f_kj[q+j-2][s[0]-1];
1674 } while (t0 != t0buf);
1676 return t0 / cln::factorial(s[0]-1);
1679 std::vector<cln::cl_N> t(j);
1686 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1687 for (int k=j-2; k>=1; k--) {
1688 t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1690 t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
1691 } while (t[0] != t0buf);
1693 return t[0] / cln::factorial(s[0]-1);
1698 cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
1700 std::vector<int> r = s;
1701 const int j = r.size();
1703 // decide on maximal size of f_kj for crandall_Z
1707 L1 = Digits * 3 + j*2;
1710 // decide on maximal size of crX for crandall_Y
1713 } else if (Digits < 86) {
1715 } else if (Digits < 192) {
1717 } else if (Digits < 394) {
1719 } else if (Digits < 808) {
1729 for (int i=0; i<j; i++) {
1738 const cln::cl_N r0factorial = cln::factorial(r[0]-1);
1740 std::vector<int> rz;
1743 for (int k=r.size()-1; k>0; k--) {
1745 rz.insert(rz.begin(), r.back());
1746 skp1buf = rz.front();
1752 for (int q=0; q<skp1buf; q++) {
1754 cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
1755 cln::cl_N pp2 = crandall_Z(rz);
1760 res = res - pp1 * pp2 / cln::factorial(q);
1762 res = res + pp1 * pp2 / cln::factorial(q);
1765 rz.front() = skp1buf;
1767 rz.insert(rz.begin(), r.back());
1771 res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
1777 cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
1779 const int j = r.size();
1781 // buffer for subsums
1782 std::vector<cln::cl_N> t(j);
1783 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1790 t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
1791 for (int k=j-2; k>=0; k--) {
1792 t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
1794 } while (t[0] != t0buf);
1800 } // end of anonymous namespace
1803 //////////////////////////////////////////////////////////////////////
1805 // Multiple zeta values zeta
1809 //////////////////////////////////////////////////////////////////////
1812 static ex zeta1_evalf(const ex& x)
1814 if (is_exactly_a<lst>(x) && (x.nops()>1)) {
1816 // multiple zeta value
1817 const int count = x.nops();
1818 const lst& xlst = ex_to<lst>(x);
1819 std::vector<int> r(count);
1821 // check parameters and convert them
1822 lst::const_iterator it1 = xlst.begin();
1823 std::vector<int>::iterator it2 = r.begin();
1825 if (!(*it1).info(info_flags::posint)) {
1826 return zeta(x).hold();
1828 *it2 = ex_to<numeric>(*it1).to_int();
1831 } while (it2 != r.end());
1833 // check for divergence
1835 return zeta(x).hold();
1838 // decide on summation algorithm
1839 // this is still a bit clumsy
1840 int limit = (Digits>17) ? 10 : 6;
1841 if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) {
1842 return numeric(zeta_do_sum_Crandall(r));
1844 return numeric(zeta_do_sum_simple(r));
1848 // single zeta value
1849 if (is_exactly_a<numeric>(x) && (x != 1)) {
1851 return zeta(ex_to<numeric>(x));
1852 } catch (const dunno &e) { }
1855 return zeta(x).hold();
1859 static ex zeta1_eval(const ex& x)
1861 if (is_exactly_a<lst>(x)) {
1862 if (x.nops() == 1) {
1863 return zeta(x.op(0));
1865 return zeta(x).hold();
1868 if (x.info(info_flags::numeric)) {
1869 const numeric& y = ex_to<numeric>(x);
1870 // trap integer arguments:
1871 if (y.is_integer()) {
1875 if (y.is_equal(_num1)) {
1876 return zeta(x).hold();
1878 if (y.info(info_flags::posint)) {
1879 if (y.info(info_flags::odd)) {
1880 return zeta(x).hold();
1882 return abs(bernoulli(y)) * pow(Pi, y) * pow(_num2, y-_num1) / factorial(y);
1885 if (y.info(info_flags::odd)) {
1886 return -bernoulli(_num1-y) / (_num1-y);
1893 if (y.info(info_flags::numeric) && !y.info(info_flags::crational))
1894 return zeta1_evalf(x);
1896 return zeta(x).hold();
1900 static ex zeta1_deriv(const ex& x, unsigned deriv_param)
1902 GINAC_ASSERT(deriv_param==0);
1904 if (is_exactly_a<lst>(x)) {
1907 return zeta(_ex1, x);
1912 unsigned zeta1_SERIAL::serial =
1913 function::register_new(function_options("zeta").
1914 eval_func(zeta1_eval).
1915 evalf_func(zeta1_evalf).
1916 do_not_evalf_params().
1917 derivative_func(zeta1_deriv).
1918 latex_name("\\zeta").
1922 //////////////////////////////////////////////////////////////////////
1924 // Multiple zeta values mZeta
1926 // The use of mZeta is deprecated! This function will be removed
1927 // from GiNaC source soon. Use zeta instead!!
1931 //////////////////////////////////////////////////////////////////////
1934 static ex mZeta_eval(const ex& x1)
1936 return mZeta(x1).hold();
1940 static ex mZeta_evalf(const ex& x1)
1942 if (is_a<lst>(x1)) {
1943 for (int i=0; i<x1.nops(); i++) {
1944 if (!x1.op(i).info(info_flags::posint))
1945 return mZeta(x1).hold();
1948 const int j = x1.nops();
1950 std::vector<int> r(j);
1951 for (int i=0; i<j; i++) {
1952 r[j-1-i] = ex_to<numeric>(x1.op(i)).to_int();
1955 // check for divergence
1957 return mZeta(x1).hold();
1960 // if only one argument, use cln::zeta
1962 return numeric(cln::zeta(r[0]));
1965 // decide on summation algorithm
1966 // this is still a bit clumsy
1967 int limit = (Digits>17) ? 10 : 6;
1968 if (r[0]<limit || (j>3 && r[1]<limit/2)) {
1969 return numeric(zeta_do_sum_Crandall(r));
1971 return numeric(zeta_do_sum_simple(r));
1973 } else if (x1.info(info_flags::posint) && (x1 != 1)) {
1974 return numeric(cln::zeta(ex_to<numeric>(x1).to_int()));
1977 return mZeta(x1).hold();
1981 static ex mZeta_deriv(const ex& x, unsigned deriv_param)
1987 REGISTER_FUNCTION(mZeta,
1988 eval_func(mZeta_eval).
1989 evalf_func(mZeta_evalf).
1990 do_not_evalf_params().
1991 derivative_func(mZeta_deriv));
1994 } // namespace GiNaC