1 /** @file inifcns_nstdsums.cpp
3 * Implementation of some special functions that have a representation as nested sums.
5 * classical polylogarithm Li(n,x)
6 * multiple polylogarithm Li(lst(n_1,...,n_k),lst(x_1,...,x_k)
7 * nielsen's generalized polylogarithm S(n,p,x)
8 * harmonic polylogarithm H(lst(m_1,...,m_k),x)
9 * multiple zeta value mZeta(lst(m_1,...,m_k))
12 * - All formulae used can be looked up in the following publication:
13 * Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
14 * This document will be referenced as [Kol] throughout this source code.
15 * - Classical polylogarithms (Li) and nielsen's generalized polylogarithms (S) can be numerically
16 * evaluated in the whole complex plane. And of course, there is still room for speed optimizations ;-).
17 * - The calculation of classical polylogarithms is speed up by using Euler-MacLaurin summation (EuMac).
18 * - The remaining functions can only be numerically evaluated with arguments lying in the unit sphere
19 * at the moment. Sorry. The evaluation especially for mZeta is very slow ... better not use it
21 * - The functions have no series expansion. To do it, you have to convert these functions
22 * into the appropriate objects from the nestedsums library, do the expansion and convert the
24 * - Numerical testing of this implementation has been performed by doing a comparison of results
25 * between this software and the commercial M.......... 4.1.
30 * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
32 * This program is free software; you can redistribute it and/or modify
33 * it under the terms of the GNU General Public License as published by
34 * the Free Software Foundation; either version 2 of the License, or
35 * (at your option) any later version.
37 * This program is distributed in the hope that it will be useful,
38 * but WITHOUT ANY WARRANTY; without even the implied warranty of
39 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
40 * GNU General Public License for more details.
42 * You should have received a copy of the GNU General Public License
43 * along with this program; if not, write to the Free Software
44 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
54 #include "operators.h"
55 #include "relational.h"
62 // lookup table for Euler-MacLaurin optimization
64 std::vector<std::vector<cln::cl_N> > Xn;
68 // lookup table for Euler-Zagier-Sums (used for S_n,p(x))
70 std::vector<std::vector<cln::cl_N> > Yn;
71 int ynsize = 0; // number of Yn[]
72 int ynlength = 100; // initial length of all Yn[i]
75 //////////////////////
76 // helper functions //
77 //////////////////////
80 // This function calculates the X_n. The X_n are needed for the Euler-MacLaurin summation (EMS) of
81 // classical polylogarithms.
82 // With EMS the polylogs can be calculated as follows:
83 // Li_p (x) = \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with u = -log(1-x)
84 // X_0(n) = B_n (Bernoulli numbers)
85 // X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
86 // The calculation of Xn depends on X0 and X{n-1}.
87 // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
88 // This results in a slightly more complicated algorithm for the X_n.
89 // The first index in Xn corresponds to the index of the polylog minus 2.
90 // The second index in Xn corresponds to the index from the EMS.
91 static void fill_Xn(int n)
93 // rule of thumb. needs to be improved. TODO
94 const int initsize = Digits * 3 / 2;
97 // calculate X_2 and higher (corresponding to Li_4 and higher)
98 std::vector<cln::cl_N> buf(initsize);
99 std::vector<cln::cl_N>::iterator it = buf.begin();
101 *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
103 for (int i=2; i<=initsize; i++) {
105 result = 0; // k == 0
107 result = Xn[0][i/2-1]; // k == 0
109 for (int k=1; k<i-1; k++) {
110 if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
111 result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
114 result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
115 result = result + Xn[n-1][i-1] / (i+1); // k == i
122 // special case to handle the X_0 correct
123 std::vector<cln::cl_N> buf(initsize);
124 std::vector<cln::cl_N>::iterator it = buf.begin();
126 *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
128 *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
130 for (int i=3; i<=initsize; i++) {
132 result = -Xn[0][(i-3)/2]/2;
133 *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
136 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
137 for (int k=1; k<i/2; k++) {
138 result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
147 std::vector<cln::cl_N> buf(initsize/2);
148 std::vector<cln::cl_N>::iterator it = buf.begin();
149 for (int i=1; i<=initsize/2; i++) {
150 *it = bernoulli(i*2).to_cl_N();
160 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
161 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
162 // representing S_{n,p}(x).
163 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
165 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
166 // representing S_{n,p}(x).
167 // The calculation of Y_n uses the values from Y_{n-1}.
168 static void fill_Yn(int n, const cln::float_format_t& prec)
170 // TODO -> get rid of the magic number
171 const int initsize = ynlength;
172 //const int initsize = initsize_Yn;
173 cln::cl_N one = cln::cl_float(1, prec);
176 std::vector<cln::cl_N> buf(initsize);
177 std::vector<cln::cl_N>::iterator it = buf.begin();
178 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
179 *it = (*itprev) / cln::cl_N(n+1) * one;
182 // sums with an index smaller than the depth are zero and need not to be calculated.
183 // calculation starts with depth, which is n+2)
184 for (int i=n+2; i<=initsize+n; i++) {
185 *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
191 std::vector<cln::cl_N> buf(initsize);
192 std::vector<cln::cl_N>::iterator it = buf.begin();
195 for (int i=2; i<=initsize; i++) {
196 *it = *(it-1) + 1 / cln::cl_N(i) * one;
205 // make Yn longer ...
206 static void make_Yn_longer(int newsize, const cln::float_format_t& prec)
209 cln::cl_N one = cln::cl_float(1, prec);
211 Yn[0].resize(newsize);
212 std::vector<cln::cl_N>::iterator it = Yn[0].begin();
214 for (int i=ynlength+1; i<=newsize; i++) {
215 *it = *(it-1) + 1 / cln::cl_N(i) * one;
219 for (int n=1; n<ynsize; n++) {
220 Yn[n].resize(newsize);
221 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
222 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
225 for (int i=ynlength+n+1; i<=newsize+n; i++) {
226 *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
236 // calculates Li(2,x) without EuMac
237 static cln::cl_N Li2_series(const cln::cl_N& x)
242 cln::cl_I den = 1; // n^2 = 1
247 den = den + i; // n^2 = 4, 9, 16, ...
249 res = res + num / den;
250 } while (res != resbuf);
255 // calculates Li(2,x) with EuMac
256 static cln::cl_N Li2_series_EuMac(const cln::cl_N& x)
258 std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
259 cln::cl_N u = -cln::log(1-x);
260 cln::cl_N factor = u;
261 cln::cl_N res = u - u*u/4;
266 factor = factor * u*u / (2*i * (2*i+1));
267 res = res + (*it) * factor;
268 it++; // should we check it? or rely on initsize? ...
270 } while (res != resbuf);
275 // calculates Li(n,x), n>2 without EuMac
276 static cln::cl_N Lin_series(int n, const cln::cl_N& x)
278 cln::cl_N factor = x;
285 res = res + factor / cln::expt(cln::cl_I(i),n);
287 } while (res != resbuf);
292 // calculates Li(n,x), n>2 with EuMac
293 static cln::cl_N Lin_series_EuMac(int n, const cln::cl_N& x)
295 std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
296 cln::cl_N u = -cln::log(1-x);
297 cln::cl_N factor = u;
303 factor = factor * u / i;
304 res = res + (*it) * factor;
305 it++; // should we check it? or rely on initsize? ...
307 } while (res != resbuf);
312 // forward declaration needed by function Li_projection and C below
313 static numeric S_num(int n, int p, const numeric& x);
316 // helper function for classical polylog Li
317 static cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
319 // treat n=2 as special case
321 // check if precalculated X0 exists
326 if (cln::realpart(x) < 0.5) {
327 // choose the faster algorithm
328 // the switching point was empirically determined. the optimal point
329 // depends on hardware, Digits, ... so an approx value is okay.
330 // it solves also the problem with precision due to the u=-log(1-x) transformation
331 if (cln::abs(cln::realpart(x)) < 0.25) {
333 return Li2_series(x);
335 return Li2_series_EuMac(x);
338 // choose the faster algorithm
339 if (cln::abs(cln::realpart(x)) > 0.75) {
340 return -Li2_series(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
342 return -Li2_series_EuMac(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
346 // check if precalculated Xn exist
348 for (int i=xnsize; i<n-1; i++) {
353 if (cln::realpart(x) < 0.5) {
354 // choose the faster algorithm
355 // with n>=12 the "normal" summation always wins against EuMac
356 if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
357 return Lin_series(n, x);
359 return Lin_series_EuMac(n, x);
362 cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
363 for (int j=0; j<n-1; j++) {
364 result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
365 * cln::expt(cln::log(x), j) / cln::factorial(j);
373 // helper function for classical polylog Li
374 static numeric Li_num(int n, const numeric& x)
378 return -cln::log(1-x.to_cl_N());
389 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
392 // what is the desired float format?
393 // first guess: default format
394 cln::float_format_t prec = cln::default_float_format;
395 const cln::cl_N value = x.to_cl_N();
396 // second guess: the argument's format
397 if (!x.real().is_rational())
398 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
399 else if (!x.imag().is_rational())
400 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
403 if (cln::abs(value) > 1) {
404 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
405 // check if argument is complex. if it is real, the new polylog has to be conjugated.
406 if (cln::zerop(cln::imagpart(value))) {
408 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
411 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
416 result = result + Li_projection(n, cln::recip(value), prec);
419 result = result - Li_projection(n, cln::recip(value), prec);
423 for (int j=0; j<n-1; j++) {
424 add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
425 * Li_num(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
427 result = result - add;
431 return Li_projection(n, value, prec);
436 // helper function for S(n,p,x)
437 static cln::cl_N numeric_nielsen(int n, int step)
441 for (int i=1; i<n; i++) {
442 res = res + numeric_nielsen(i, step-1) / cln::cl_I(i);
452 // helper function for S(n,p,x)
454 static cln::cl_N C(int n, int p)
458 for (int k=0; k<p; k++) {
459 for (int j=0; j<=(n+k-1)/2; j++) {
463 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
466 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
473 result = result + cln::factorial(n+k-1)
474 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
475 / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
478 result = result - cln::factorial(n+k-1)
479 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
480 / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
485 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()
486 / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
489 result = result + cln::factorial(n+k-1)
490 * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
491 / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
499 if (((np)/2+n) & 1) {
500 result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
503 result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
511 // helper function for S(n,p,x)
512 // [Kol] remark to (9.1)
513 static cln::cl_N a_k(int k)
522 for (int m=2; m<=k; m++) {
523 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
530 // helper function for S(n,p,x)
531 // [Kol] remark to (9.1)
532 static cln::cl_N b_k(int k)
541 for (int m=2; m<=k; m++) {
542 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
549 // helper function for S(n,p,x)
550 static cln::cl_N S_series(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
553 return Li_projection(n+1, x, prec);
556 // TODO -> check for vector boundaries and do missing calculations
558 // check if precalculated values are sufficient
560 for (int i=ynsize; i<p-1; i++) {
565 // should be done otherwise
566 cln::cl_N xf = x * cln::cl_float(1, prec);
569 cln::cl_N resultbuffer;
571 for (i=p; true; i++) {
572 resultbuffer = result;
573 if (i-p >= ynlength) {
575 make_Yn_longer(ynlength*2, prec);
577 result = result + cln::expt(xf,i) / cln::expt(cln::cl_I(i),n+1) * Yn[p-2][i-p]; // should we check it? or rely on magic number? ...
578 if (cln::zerop(result-resultbuffer)) {
587 // helper function for S(n,p,x)
588 static cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
591 if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
593 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
594 * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
596 for (int s=0; s<n; s++) {
598 for (int r=0; r<p; r++) {
599 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
600 * S_series(p-r,n-s,1-x,prec) / cln::factorial(r);
602 result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
608 return S_series(n, p, x, prec);
612 // helper function for S(n,p,x)
613 static numeric S_num(int n, int p, const numeric& x)
617 // [Kol] (2.22) with (2.21)
618 return cln::zeta(p+1);
623 return cln::zeta(n+1);
628 for (int nu=0; nu<n; nu++) {
629 for (int rho=0; rho<=p; rho++) {
630 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
631 * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
634 result = result * cln::expt(cln::cl_I(-1),n+p-1);
641 return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
643 // throw std::runtime_error("don't know how to evaluate this function!");
646 // what is the desired float format?
647 // first guess: default format
648 cln::float_format_t prec = cln::default_float_format;
649 const cln::cl_N value = x.to_cl_N();
650 // second guess: the argument's format
651 if (!x.real().is_rational())
652 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
653 else if (!x.imag().is_rational())
654 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
658 if (cln::realpart(value) < -0.5) {
660 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
661 * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
663 for (int s=0; s<n; s++) {
665 for (int r=0; r<p; r++) {
666 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
667 * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
669 result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
676 if (cln::abs(value) > 1) {
680 for (int s=0; s<p; s++) {
681 for (int r=0; r<=s; r++) {
682 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
683 / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
684 * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
687 result = result * cln::expt(cln::cl_I(-1),n);
690 for (int r=0; r<n; r++) {
691 res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
693 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
695 result = result + cln::expt(cln::cl_I(-1),p) * res2;
700 return S_projection(n, p, value, prec);
705 // helper function for multiple polylogarithm
706 static cln::cl_N numeric_zsum(int n, std::vector<cln::cl_N>& x, std::vector<cln::cl_N>& m)
712 for (int i=1; i<n; i++) {
713 std::vector<cln::cl_N>::iterator be;
714 std::vector<cln::cl_N>::iterator en;
718 std::vector<cln::cl_N> xbuf(be, en);
722 std::vector<cln::cl_N> mbuf(be, en);
723 res = res + cln::expt(x[0],i) / cln::expt(i,m[0]) * numeric_zsum(i, xbuf, mbuf);
729 // helper function for harmonic polylogarithm
730 static cln::cl_N numeric_harmonic(int n, std::vector<cln::cl_N>& m)
736 for (int i=1; i<n; i++) {
737 std::vector<cln::cl_N>::iterator be;
738 std::vector<cln::cl_N>::iterator en;
742 std::vector<cln::cl_N> mbuf(be, en);
743 res = res + cln::recip(cln::expt(i,m[0])) * numeric_harmonic(i, mbuf);
749 /////////////////////////////
750 // end of helper functions //
751 /////////////////////////////
754 // Polylogarithm and multiple polylogarithm
756 static ex Li_eval(const ex& x1, const ex& x2)
762 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational)))
763 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
764 return Li(x1,x2).hold();
768 static ex Li_evalf(const ex& x1, const ex& x2)
770 // classical polylogs
771 if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
772 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
775 else if (is_a<lst>(x1) && is_a<lst>(x2)) {
776 for (int i=0; i<x1.nops(); i++) {
777 if (!is_a<numeric>(x1.op(i)))
778 return Li(x1,x2).hold();
779 if (!is_a<numeric>(x2.op(i)))
780 return Li(x1,x2).hold();
782 return Li(x1,x2).hold();
785 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
786 cln::cl_N x_1 = ex_to<numeric>(x2.op(x2.nops()-1)).to_cl_N();
787 std::vector<cln::cl_N> x;
788 std::vector<cln::cl_N> m;
789 const int nops = ex_to<numeric>(x1.nops()).to_int();
790 for (int i=nops-2; i>=0; i--) {
791 m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
792 x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
797 for (int i=nops; true; i++) {
799 res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_zsum(i, x, m);
800 if (cln::zerop(res-resbuf))
808 return Li(x1,x2).hold();
811 static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
814 seq.push_back(expair(Li(x1,x2), 0));
815 return pseries(rel,seq);
818 REGISTER_FUNCTION(Li, eval_func(Li_eval).evalf_func(Li_evalf).do_not_evalf_params().series_func(Li_series));
821 // Nielsen's generalized polylogarithm
823 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
828 if (x3.info(info_flags::numeric) && (!x3.info(info_flags::crational)) &&
829 x1.info(info_flags::posint) && x2.info(info_flags::posint)) {
830 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
832 return S(x1,x2,x3).hold();
835 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
837 if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
838 if ((x3 == -1) && (x2 != 1)) {
839 // no formula to evaluate this ... sorry
840 // return S(x1,x2,x3).hold();
842 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
844 return S(x1,x2,x3).hold();
847 static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options)
850 seq.push_back(expair(S(x1,x2,x3), 0));
851 return pseries(rel,seq);
854 REGISTER_FUNCTION(S, eval_func(S_eval).evalf_func(S_evalf).do_not_evalf_params().series_func(S_series));
857 // Harmonic polylogarithm
859 static ex H_eval(const ex& x1, const ex& x2)
861 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational))) {
862 return H(x1,x2).evalf();
864 return H(x1,x2).hold();
867 static ex H_evalf(const ex& x1, const ex& x2)
869 if (is_a<lst>(x1) && is_a<numeric>(x2)) {
870 for (int i=0; i<x1.nops(); i++) {
871 if (!is_a<numeric>(x1.op(i)))
872 return H(x1,x2).hold();
875 return H(x1,x2).hold();
878 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
879 cln::cl_N x_1 = ex_to<numeric>(x2).to_cl_N();
880 std::vector<cln::cl_N> m;
881 const int nops = ex_to<numeric>(x1.nops()).to_int();
882 for (int i=nops-2; i>=0; i--) {
883 m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
888 for (int i=nops; true; i++) {
890 res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_harmonic(i, m);
891 if (cln::zerop(res-resbuf))
899 return H(x1,x2).hold();
902 static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
905 seq.push_back(expair(H(x1,x2), 0));
906 return pseries(rel,seq);
909 REGISTER_FUNCTION(H, eval_func(H_eval).evalf_func(H_evalf).do_not_evalf_params().series_func(H_series));
912 // Multiple zeta value
914 static ex mZeta_eval(const ex& x1)
916 return mZeta(x1).hold();
919 static ex mZeta_evalf(const ex& x1)
922 for (int i=0; i<x1.nops(); i++) {
923 if (!is_a<numeric>(x1.op(i)))
924 return mZeta(x1).hold();
927 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
929 // check for divergence
931 return mZeta(x1).hold();
934 std::vector<cln::cl_N> m;
935 const int nops = ex_to<numeric>(x1.nops()).to_int();
936 for (int i=nops-2; i>=0; i--) {
937 m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
940 cln::float_format_t prec = cln::default_float_format;
941 cln::cl_N res = cln::complex(cln::cl_float(0, prec), 0);
943 for (int i=nops; true; i++) {
944 // to infinity and beyond ... timewise
946 res = res + cln::recip(cln::expt(i,m_1)) * numeric_harmonic(i, m);
947 if (cln::zerop(res-resbuf))
955 return mZeta(x1).hold();
958 static ex mZeta_series(const ex& x1, const relational& rel, int order, unsigned options)
961 seq.push_back(expair(mZeta(x1), 0));
962 return pseries(rel,seq);
965 REGISTER_FUNCTION(mZeta, eval_func(mZeta_eval).evalf_func(mZeta_evalf).do_not_evalf_params().series_func(mZeta_series));