fbc4a3faba81e4f24cb2a53ab14a1dc6d03d8bda
[ginac.git] / ginac / inifcns_nstdsums.cpp
1 /** @file inifcns_nstdsums.cpp
2  *
3  *  Implementation of some special functions that have a representation as nested sums.
4  *  
5  *  The functions are: 
6  *    classical polylogarithm              Li(n,x)
7  *    multiple polylogarithm               Li(lst(m_1,...,m_k),lst(x_1,...,x_k))
8  *    nielsen's generalized polylogarithm  S(n,p,x)
9  *    harmonic polylogarithm               H(m,x) or H(lst(m_1,...,m_k),x)
10  *    multiple zeta value                  zeta(m) or zeta(lst(m_1,...,m_k))
11  *    alternating Euler sum                zeta(m,s) or zeta(lst(m_1,...,m_k),lst(s_1,...,s_k))
12  *
13  *  Some remarks:
14  *    
15  *    - All formulae used can be looked up in the following publications:
16  *      [Kol] Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
17  *      [Cra] Fast Evaluation of Multiple Zeta Sums, R.E.Crandall, Math.Comp. 67 (1998), pp. 1163-1172.
18  *      [ReV] Harmonic Polylogarithms, E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
19  *      [BBB] Special Values of Multiple Polylogarithms, J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001), pp. 907-941
20  *
21  *    - The order of parameters and arguments of H, Li and zeta is defined according to their order in the
22  *      nested sums representation.
23  *      
24  *    - Except for the multiple polylogarithm all functions can be nummerically evaluated with arguments in
25  *      the whole complex plane. Multiple polylogarithms evaluate only if each argument x_i is smaller than
26  *      one. The parameters for every function (n, p or n_i) must be positive integers.
27  *      
28  *    - The calculation of classical polylogarithms is speed up by using Bernoulli numbers and 
29  *      look-up tables. S uses look-up tables as well. The zeta function applies the algorithms in
30  *      [Cra] and [BBB] for speed up.
31  *      
32  *    - The functions have no series expansion as nested sums. To do it, you have to convert these functions
33  *      into the appropriate objects from the nestedsums library, do the expansion and convert the
34  *      result back. 
35  *      
36  *    - Numerical testing of this implementation has been performed by doing a comparison of results
37  *      between this software and the commercial M.......... 4.1. Multiple zeta values have been checked
38  *      by means of evaluations into simple zeta values. Harmonic polylogarithms have been checked by
39  *      comparison to S(n,p,x) for corresponding parameter combinations and by continuity checks
40  *      around |x|=1 along with comparisons to corresponding zeta functions.
41  *
42  */
43
44 /*
45  *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
46  *
47  *  This program is free software; you can redistribute it and/or modify
48  *  it under the terms of the GNU General Public License as published by
49  *  the Free Software Foundation; either version 2 of the License, or
50  *  (at your option) any later version.
51  *
52  *  This program is distributed in the hope that it will be useful,
53  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
54  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
55  *  GNU General Public License for more details.
56  *
57  *  You should have received a copy of the GNU General Public License
58  *  along with this program; if not, write to the Free Software
59  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
60  */
61
62 #include <stdexcept>
63 #include <vector>
64 #include <cln/cln.h>
65
66 #include "inifcns.h"
67
68 #include "add.h"
69 #include "constant.h"
70 #include "lst.h"
71 #include "mul.h"
72 #include "numeric.h"
73 #include "operators.h"
74 #include "power.h"
75 #include "pseries.h"
76 #include "relational.h"
77 #include "symbol.h"
78 #include "utils.h"
79 #include "wildcard.h"
80
81
82 namespace GiNaC {
83
84
85 //////////////////////////////////////////////////////////////////////
86 //
87 // Classical polylogarithm  Li(n,x)
88 //
89 // helper functions
90 //
91 //////////////////////////////////////////////////////////////////////
92
93
94 // anonymous namespace for helper functions
95 namespace {
96
97
98 // lookup table for factors built from Bernoulli numbers
99 // see fill_Xn()
100 std::vector<std::vector<cln::cl_N> > Xn;
101 int xnsize = 0;
102
103
104 // This function calculates the X_n. The X_n are needed for speed up of classical polylogarithms.
105 // With these numbers the polylogs can be calculated as follows:
106 //   Li_p (x)  =  \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with  u = -log(1-x)
107 //   X_0(n) = B_n (Bernoulli numbers)
108 //   X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
109 // The calculation of Xn depends on X0 and X{n-1}.
110 // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
111 // This results in a slightly more complicated algorithm for the X_n.
112 // The first index in Xn corresponds to the index of the polylog minus 2.
113 // The second index in Xn corresponds to the index from the actual sum.
114 void fill_Xn(int n)
115 {
116         // rule of thumb. needs to be improved. TODO
117         const int initsize = Digits * 3 / 2;
118
119         if (n>1) {
120                 // calculate X_2 and higher (corresponding to Li_4 and higher)
121                 std::vector<cln::cl_N> buf(initsize);
122                 std::vector<cln::cl_N>::iterator it = buf.begin();
123                 cln::cl_N result;
124                 *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
125                 it++;
126                 for (int i=2; i<=initsize; i++) {
127                         if (i&1) {
128                                 result = 0; // k == 0
129                         } else {
130                                 result = Xn[0][i/2-1]; // k == 0
131                         }
132                         for (int k=1; k<i-1; k++) {
133                                 if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
134                                         result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
135                                 }
136                         }
137                         result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
138                         result = result + Xn[n-1][i-1] / (i+1); // k == i
139                         
140                         *it = result;
141                         it++;
142                 }
143                 Xn.push_back(buf);
144         } else if (n==1) {
145                 // special case to handle the X_0 correct
146                 std::vector<cln::cl_N> buf(initsize);
147                 std::vector<cln::cl_N>::iterator it = buf.begin();
148                 cln::cl_N result;
149                 *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
150                 it++;
151                 *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
152                 it++;
153                 for (int i=3; i<=initsize; i++) {
154                         if (i & 1) {
155                                 result = -Xn[0][(i-3)/2]/2;
156                                 *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
157                                 it++;
158                         } else {
159                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
160                                 for (int k=1; k<i/2; k++) {
161                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
162                                 }
163                                 *it = result;
164                                 it++;
165                         }
166                 }
167                 Xn.push_back(buf);
168         } else {
169                 // calculate X_0
170                 std::vector<cln::cl_N> buf(initsize/2);
171                 std::vector<cln::cl_N>::iterator it = buf.begin();
172                 for (int i=1; i<=initsize/2; i++) {
173                         *it = bernoulli(i*2).to_cl_N();
174                         it++;
175                 }
176                 Xn.push_back(buf);
177         }
178
179         xnsize++;
180 }
181
182
183 // calculates Li(2,x) without Xn
184 cln::cl_N Li2_do_sum(const cln::cl_N& x)
185 {
186         cln::cl_N res = x;
187         cln::cl_N resbuf;
188         cln::cl_N num = x;
189         cln::cl_I den = 1; // n^2 = 1
190         unsigned i = 3;
191         do {
192                 resbuf = res;
193                 num = num * x;
194                 den = den + i;  // n^2 = 4, 9, 16, ...
195                 i += 2;
196                 res = res + num / den;
197         } while (res != resbuf);
198         return res;
199 }
200
201
202 // calculates Li(2,x) with Xn
203 cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x)
204 {
205         std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
206         cln::cl_N u = -cln::log(1-x);
207         cln::cl_N factor = u;
208         cln::cl_N res = u - u*u/4;
209         cln::cl_N resbuf;
210         unsigned i = 1;
211         do {
212                 resbuf = res;
213                 factor = factor * u*u / (2*i * (2*i+1));
214                 res = res + (*it) * factor;
215                 it++; // should we check it? or rely on initsize? ...
216                 i++;
217         } while (res != resbuf);
218         return res;
219 }
220
221
222 // calculates Li(n,x), n>2 without Xn
223 cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
224 {
225         cln::cl_N factor = x;
226         cln::cl_N res = x;
227         cln::cl_N resbuf;
228         int i=2;
229         do {
230                 resbuf = res;
231                 factor = factor * x;
232                 res = res + factor / cln::expt(cln::cl_I(i),n);
233                 i++;
234         } while (res != resbuf);
235         return res;
236 }
237
238
239 // calculates Li(n,x), n>2 with Xn
240 cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
241 {
242         std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
243         cln::cl_N u = -cln::log(1-x);
244         cln::cl_N factor = u;
245         cln::cl_N res = u;
246         cln::cl_N resbuf;
247         unsigned i=2;
248         do {
249                 resbuf = res;
250                 factor = factor * u / i;
251                 res = res + (*it) * factor;
252                 it++; // should we check it? or rely on initsize? ...
253                 i++;
254         } while (res != resbuf);
255         return res;
256 }
257
258
259 // forward declaration needed by function Li_projection and C below
260 numeric S_num(int n, int p, const numeric& x);
261
262
263 // helper function for classical polylog Li
264 cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
265 {
266         // treat n=2 as special case
267         if (n == 2) {
268                 // check if precalculated X0 exists
269                 if (xnsize == 0) {
270                         fill_Xn(0);
271                 }
272
273                 if (cln::realpart(x) < 0.5) {
274                         // choose the faster algorithm
275                         // the switching point was empirically determined. the optimal point
276                         // depends on hardware, Digits, ... so an approx value is okay.
277                         // it solves also the problem with precision due to the u=-log(1-x) transformation
278                         if (cln::abs(cln::realpart(x)) < 0.25) {
279                                 
280                                 return Li2_do_sum(x);
281                         } else {
282                                 return Li2_do_sum_Xn(x);
283                         }
284                 } else {
285                         // choose the faster algorithm
286                         if (cln::abs(cln::realpart(x)) > 0.75) {
287                                 return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
288                         } else {
289                                 return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
290                         }
291                 }
292         } else {
293                 // check if precalculated Xn exist
294                 if (n > xnsize+1) {
295                         for (int i=xnsize; i<n-1; i++) {
296                                 fill_Xn(i);
297                         }
298                 }
299
300                 if (cln::realpart(x) < 0.5) {
301                         // choose the faster algorithm
302                         // with n>=12 the "normal" summation always wins against the method with Xn
303                         if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
304                                 return Lin_do_sum(n, x);
305                         } else {
306                                 return Lin_do_sum_Xn(n, x);
307                         }
308                 } else {
309                         cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
310                         for (int j=0; j<n-1; j++) {
311                                 result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
312                                         * cln::expt(cln::log(x), j) / cln::factorial(j);
313                         }
314                         return result;
315                 }
316         }
317 }
318
319
320 // helper function for classical polylog Li
321 numeric Li_num(int n, const numeric& x)
322 {
323         if (n == 1) {
324                 // just a log
325                 return -cln::log(1-x.to_cl_N());
326         }
327         if (x.is_zero()) {
328                 return 0;
329         }
330         if (x == 1) {
331                 // [Kol] (2.22)
332                 return cln::zeta(n);
333         }
334         else if (x == -1) {
335                 // [Kol] (2.22)
336                 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
337         }
338         
339         // what is the desired float format?
340         // first guess: default format
341         cln::float_format_t prec = cln::default_float_format;
342         const cln::cl_N value = x.to_cl_N();
343         // second guess: the argument's format
344         if (!x.real().is_rational())
345                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
346         else if (!x.imag().is_rational())
347                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
348         
349         // [Kol] (5.15)
350         if (cln::abs(value) > 1) {
351                 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
352                 // check if argument is complex. if it is real, the new polylog has to be conjugated.
353                 if (cln::zerop(cln::imagpart(value))) {
354                         if (n & 1) {
355                                 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
356                         }
357                         else {
358                                 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
359                         }
360                 }
361                 else {
362                         if (n & 1) {
363                                 result = result + Li_projection(n, cln::recip(value), prec);
364                         }
365                         else {
366                                 result = result - Li_projection(n, cln::recip(value), prec);
367                         }
368                 }
369                 cln::cl_N add;
370                 for (int j=0; j<n-1; j++) {
371                         add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
372                                         * Li_num(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
373                 }
374                 result = result - add;
375                 return result;
376         }
377         else {
378                 return Li_projection(n, value, prec);
379         }
380 }
381
382
383 } // end of anonymous namespace
384
385
386 //////////////////////////////////////////////////////////////////////
387 //
388 // Multiple polylogarithm  Li(n,x)
389 //
390 // helper function
391 //
392 //////////////////////////////////////////////////////////////////////
393
394
395 // anonymous namespace for helper function
396 namespace {
397
398
399 cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
400 {
401         const int j = s.size();
402
403         std::vector<cln::cl_N> t(j);
404         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
405
406         cln::cl_N t0buf;
407         int q = 0;
408         do {
409                 t0buf = t[0];
410                 q++;
411                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
412                 for (int k=j-2; k>=0; k--) {
413                         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]);
414                 }
415         } while ((t[0] != t0buf) || (q<10));
416         
417         return t[0];
418 }
419
420
421 } // end of anonymous namespace
422
423
424 //////////////////////////////////////////////////////////////////////
425 //
426 // Classical polylogarithm and multiple polylogarithm  Li(n,x)
427 //
428 // GiNaC function
429 //
430 //////////////////////////////////////////////////////////////////////
431
432
433 static ex Li_eval(const ex& x1, const ex& x2)
434 {
435         if (x2.is_zero()) {
436                 return _ex0;
437         }
438         else {
439                 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational)))
440                         return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
441                 if (is_a<lst>(x2)) {
442                         for (int i=0; i<x2.nops(); i++) {
443                                 if (!is_a<numeric>(x2.op(i))) {
444                                         return Li(x1,x2).hold();
445                                 }
446                         }
447                         return Li(x1,x2).evalf();
448                 }
449                 return Li(x1,x2).hold();
450         }
451 }
452
453
454 static ex Li_evalf(const ex& x1, const ex& x2)
455 {
456         // classical polylogs
457         if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
458                 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
459         }
460         // multiple polylogs
461         else if (is_a<lst>(x1) && is_a<lst>(x2)) {
462                 ex conv = 1;
463                 for (int i=0; i<x1.nops(); i++) {
464                         if (!x1.op(i).info(info_flags::posint)) {
465                                 return Li(x1,x2).hold();
466                         }
467                         if (!is_a<numeric>(x2.op(i))) {
468                                 return Li(x1,x2).hold();
469                         }
470                         conv *= x2.op(i);
471                         if ((conv > 1) || ((conv == 1) && (x1.op(0) == 1))) {
472                                 return Li(x1,x2).hold();
473                         }
474                 }
475
476                 std::vector<int> m;
477                 std::vector<cln::cl_N> x;
478                 for (int i=0; i<ex_to<numeric>(x1.nops()).to_int(); i++) {
479                         m.push_back(ex_to<numeric>(x1.op(i)).to_int());
480                         x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
481                 }
482
483                 return numeric(multipleLi_do_sum(m, x));
484         }
485
486         return Li(x1,x2).hold();
487 }
488
489
490 static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
491 {
492         epvector seq;
493         seq.push_back(expair(Li(x1,x2), 0));
494         return pseries(rel,seq);
495 }
496
497
498 static ex Li_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
499 {
500         GINAC_ASSERT(deriv_param < 2);
501         if (deriv_param == 0) {
502                 return _ex0;
503         }
504         if (x1 > 0) {
505                 return Li(x1-1, x2) / x2;
506         } else {
507                 return 1/(1-x2);
508         }
509 }
510
511
512 REGISTER_FUNCTION(Li,
513                 eval_func(Li_eval).
514                 evalf_func(Li_evalf).
515                 do_not_evalf_params().
516                 series_func(Li_series).
517                 derivative_func(Li_deriv));
518
519
520 //////////////////////////////////////////////////////////////////////
521 //
522 // Nielsen's generalized polylogarithm  S(n,p,x)
523 //
524 // helper functions
525 //
526 //////////////////////////////////////////////////////////////////////
527
528
529 // anonymous namespace for helper functions
530 namespace {
531
532
533 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
534 // see fill_Yn()
535 std::vector<std::vector<cln::cl_N> > Yn;
536 int ynsize = 0; // number of Yn[]
537 int ynlength = 100; // initial length of all Yn[i]
538
539
540 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
541 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
542 // representing S_{n,p}(x).
543 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
544 // equivalent Z-sum.
545 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
546 // representing S_{n,p}(x).
547 // The calculation of Y_n uses the values from Y_{n-1}.
548 void fill_Yn(int n, const cln::float_format_t& prec)
549 {
550         const int initsize = ynlength;
551         //const int initsize = initsize_Yn;
552         cln::cl_N one = cln::cl_float(1, prec);
553
554         if (n) {
555                 std::vector<cln::cl_N> buf(initsize);
556                 std::vector<cln::cl_N>::iterator it = buf.begin();
557                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
558                 *it = (*itprev) / cln::cl_N(n+1) * one;
559                 it++;
560                 itprev++;
561                 // sums with an index smaller than the depth are zero and need not to be calculated.
562                 // calculation starts with depth, which is n+2)
563                 for (int i=n+2; i<=initsize+n; i++) {
564                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
565                         it++;
566                         itprev++;
567                 }
568                 Yn.push_back(buf);
569         } else {
570                 std::vector<cln::cl_N> buf(initsize);
571                 std::vector<cln::cl_N>::iterator it = buf.begin();
572                 *it = 1 * one;
573                 it++;
574                 for (int i=2; i<=initsize; i++) {
575                         *it = *(it-1) + 1 / cln::cl_N(i) * one;
576                         it++;
577                 }
578                 Yn.push_back(buf);
579         }
580         ynsize++;
581 }
582
583
584 // make Yn longer ... 
585 void make_Yn_longer(int newsize, const cln::float_format_t& prec)
586 {
587
588         cln::cl_N one = cln::cl_float(1, prec);
589
590         Yn[0].resize(newsize);
591         std::vector<cln::cl_N>::iterator it = Yn[0].begin();
592         it += ynlength;
593         for (int i=ynlength+1; i<=newsize; i++) {
594                 *it = *(it-1) + 1 / cln::cl_N(i) * one;
595                 it++;
596         }
597
598         for (int n=1; n<ynsize; n++) {
599                 Yn[n].resize(newsize);
600                 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
601                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
602                 it += ynlength;
603                 itprev += ynlength;
604                 for (int i=ynlength+n+1; i<=newsize+n; i++) {
605                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
606                         it++;
607                         itprev++;
608                 }
609         }
610         
611         ynlength = newsize;
612 }
613
614
615 // helper function for S(n,p,x)
616 // [Kol] (7.2)
617 cln::cl_N C(int n, int p)
618 {
619         cln::cl_N result;
620
621         for (int k=0; k<p; k++) {
622                 for (int j=0; j<=(n+k-1)/2; j++) {
623                         if (k == 0) {
624                                 if (n & 1) {
625                                         if (j & 1) {
626                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
627                                         }
628                                         else {
629                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
630                                         }
631                                 }
632                         }
633                         else {
634                                 if (k & 1) {
635                                         if (j & 1) {
636                                                 result = result + cln::factorial(n+k-1)
637                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
638                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
639                                         }
640                                         else {
641                                                 result = result - cln::factorial(n+k-1)
642                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
643                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
644                                         }
645                                 }
646                                 else {
647                                         if (j & 1) {
648                                                 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()
649                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
650                                         }
651                                         else {
652                                                 result = result + cln::factorial(n+k-1)
653                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
654                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
655                                         }
656                                 }
657                         }
658                 }
659         }
660         int np = n+p;
661         if ((np-1) & 1) {
662                 if (((np)/2+n) & 1) {
663                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
664                 }
665                 else {
666                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
667                 }
668         }
669
670         return result;
671 }
672
673
674 // helper function for S(n,p,x)
675 // [Kol] remark to (9.1)
676 cln::cl_N a_k(int k)
677 {
678         cln::cl_N result;
679
680         if (k == 0) {
681                 return 1;
682         }
683
684         result = result;
685         for (int m=2; m<=k; m++) {
686                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
687         }
688
689         return -result / k;
690 }
691
692
693 // helper function for S(n,p,x)
694 // [Kol] remark to (9.1)
695 cln::cl_N b_k(int k)
696 {
697         cln::cl_N result;
698
699         if (k == 0) {
700                 return 1;
701         }
702
703         result = result;
704         for (int m=2; m<=k; m++) {
705                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
706         }
707
708         return result / k;
709 }
710
711
712 // helper function for S(n,p,x)
713 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
714 {
715         if (p==1) {
716                 return Li_projection(n+1, x, prec);
717         }
718         
719         // check if precalculated values are sufficient
720         if (p > ynsize+1) {
721                 for (int i=ynsize; i<p-1; i++) {
722                         fill_Yn(i, prec);
723                 }
724         }
725
726         // should be done otherwise
727         cln::cl_N xf = x * cln::cl_float(1, prec);
728
729         cln::cl_N res;
730         cln::cl_N resbuf;
731         cln::cl_N factor = cln::expt(xf, p);
732         int i = p;
733         do {
734                 resbuf = res;
735                 if (i-p >= ynlength) {
736                         // make Yn longer
737                         make_Yn_longer(ynlength*2, prec);
738                 }
739                 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? ...
740                 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
741                 factor = factor * xf;
742                 i++;
743         } while (res != resbuf);
744         
745         return res;
746 }
747
748
749 // helper function for S(n,p,x)
750 cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
751 {
752         // [Kol] (5.3)
753         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
754
755                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
756                         * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
757
758                 for (int s=0; s<n; s++) {
759                         cln::cl_N res2;
760                         for (int r=0; r<p; r++) {
761                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
762                                         * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
763                         }
764                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
765                 }
766
767                 return result;
768         }
769         
770         return S_do_sum(n, p, x, prec);
771 }
772
773
774 // helper function for S(n,p,x)
775 numeric S_num(int n, int p, const numeric& x)
776 {
777         if (x == 1) {
778                 if (n == 1) {
779                     // [Kol] (2.22) with (2.21)
780                         return cln::zeta(p+1);
781                 }
782
783                 if (p == 1) {
784                     // [Kol] (2.22)
785                         return cln::zeta(n+1);
786                 }
787
788                 // [Kol] (9.1)
789                 cln::cl_N result;
790                 for (int nu=0; nu<n; nu++) {
791                         for (int rho=0; rho<=p; rho++) {
792                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
793                                         * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
794                         }
795                 }
796                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
797
798                 return result;
799         }
800         else if (x == -1) {
801                 // [Kol] (2.22)
802                 if (p == 1) {
803                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
804                 }
805 //              throw std::runtime_error("don't know how to evaluate this function!");
806         }
807
808         // what is the desired float format?
809         // first guess: default format
810         cln::float_format_t prec = cln::default_float_format;
811         const cln::cl_N value = x.to_cl_N();
812         // second guess: the argument's format
813         if (!x.real().is_rational())
814                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
815         else if (!x.imag().is_rational())
816                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
817
818
819         // [Kol] (5.3)
820         if (cln::realpart(value) < -0.5) {
821
822                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
823                         * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
824
825                 for (int s=0; s<n; s++) {
826                         cln::cl_N res2;
827                         for (int r=0; r<p; r++) {
828                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
829                                         * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
830                         }
831                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
832                 }
833
834                 return result;
835                 
836         }
837         // [Kol] (5.12)
838         if (cln::abs(value) > 1) {
839                 
840                 cln::cl_N result;
841
842                 for (int s=0; s<p; s++) {
843                         for (int r=0; r<=s; r++) {
844                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
845                                         / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
846                                         * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
847                         }
848                 }
849                 result = result * cln::expt(cln::cl_I(-1),n);
850
851                 cln::cl_N res2;
852                 for (int r=0; r<n; r++) {
853                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
854                 }
855                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
856
857                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
858
859                 return result;
860         }
861         else {
862                 return S_projection(n, p, value, prec);
863         }
864 }
865
866
867 } // end of anonymous namespace
868
869
870 //////////////////////////////////////////////////////////////////////
871 //
872 // Nielsen's generalized polylogarithm  S(n,p,x)
873 //
874 // GiNaC function
875 //
876 //////////////////////////////////////////////////////////////////////
877
878
879 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
880 {
881         if (x2 == 1) {
882                 return Li(x1+1,x3);
883         }
884         if (x3.info(info_flags::numeric) && (!x3.info(info_flags::crational)) && 
885                         x1.info(info_flags::posint) && x2.info(info_flags::posint)) {
886                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
887         }
888         return S(x1,x2,x3).hold();
889 }
890
891
892 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
893 {
894         if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
895                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
896         }
897         return S(x1,x2,x3).hold();
898 }
899
900
901 static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options)
902 {
903         epvector seq;
904         seq.push_back(expair(S(x1,x2,x3), 0));
905         return pseries(rel,seq);
906 }
907
908
909 static ex S_deriv(const ex& x1, const ex& x2, const ex& x3, unsigned deriv_param)
910 {
911         GINAC_ASSERT(deriv_param < 3);
912         if (deriv_param < 2) {
913                 return _ex0;
914         }
915         if (x1 > 0) {
916                 return S(x1-1, x2, x3) / x3;
917         } else {
918                 return S(x1, x2-1, x3) / (1-x3);
919         }
920 }
921
922
923 REGISTER_FUNCTION(S,
924                 eval_func(S_eval).
925                 evalf_func(S_evalf).
926                 do_not_evalf_params().
927                 series_func(S_series).
928                 derivative_func(S_deriv));
929
930
931 //////////////////////////////////////////////////////////////////////
932 //
933 // Harmonic polylogarithm  H(m,x)
934 //
935 // helper function
936 //
937 //////////////////////////////////////////////////////////////////////
938
939
940 // anonymous namespace for helper functions
941 namespace {
942
943
944 // forward declaration
945 ex convert_from_RV(const lst& parameterlst, const ex& arg);
946
947
948 // multiplies an one-dimensional H with another H
949 // [ReV] (18)
950 ex trafo_H_mult(const ex& h1, const ex& h2)
951 {
952         ex res;
953         ex hshort;
954         lst hlong;
955         ex h1nops = h1.op(0).nops();
956         ex h2nops = h2.op(0).nops();
957         if (h1nops > 1) {
958                 hshort = h2.op(0).op(0);
959                 hlong = ex_to<lst>(h1.op(0));
960         } else {
961                 hshort = h1.op(0).op(0);
962                 if (h2nops > 1) {
963                         hlong = ex_to<lst>(h2.op(0));
964                 } else {
965                         hlong = h2.op(0).op(0);
966                 }
967         }
968         for (int i=0; i<=hlong.nops(); i++) {
969                 lst newparameter;
970                 int j=0;
971                 for (; j<i; j++) {
972                         newparameter.append(hlong[j]);
973                 }
974                 newparameter.append(hshort);
975                 for (; j<hlong.nops(); j++) {
976                         newparameter.append(hlong[j]);
977                 }
978                 res += H(newparameter, h1.op(1)).hold();
979         }
980         return res;
981 }
982
983
984 // applies trafo_H_mult recursively on expressions
985 struct map_trafo_H_mult : public map_function
986 {
987         ex operator()(const ex& e)
988         {
989                 if (is_a<add>(e)) {
990                         return e.map(*this);
991                 }
992
993                 if (is_a<mul>(e)) {
994
995                         ex result = 1;
996                         ex firstH;
997                         lst Hlst;
998                         for (int pos=0; pos<e.nops(); pos++) {
999                                 if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
1000                                         std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
1001                                         if (name == "H") {
1002                                                 for (ex i=0; i<e.op(pos).op(1); i++) {
1003                                                         Hlst.append(e.op(pos).op(0));
1004                                                 }
1005                                                 continue;
1006                                         }
1007                                 } else if (is_a<function>(e.op(pos))) {
1008                                         std::string name = ex_to<function>(e.op(pos)).get_name();
1009                                         if (name == "H") {
1010                                                 if (e.op(pos).op(0).nops() > 1) {
1011                                                         firstH = e.op(pos);
1012                                                 } else {
1013                                                         Hlst.append(e.op(pos));
1014                                                 }
1015                                                 continue;
1016                                         }
1017                                 }
1018                                 result *= e.op(pos);
1019                         }
1020                         if (firstH == 0) {
1021                                 if (Hlst.nops() > 0) {
1022                                         firstH = Hlst[Hlst.nops()-1];
1023                                         Hlst.remove_last();
1024                                 } else {
1025                                         return e;
1026                                 }
1027                         }
1028
1029                         if (Hlst.nops() > 0) {
1030                                 ex buffer = trafo_H_mult(firstH, Hlst.op(0));
1031                                 result *= buffer;
1032                                 for (int i=1; i<Hlst.nops(); i++) {
1033                                         result *= Hlst.op(i);
1034                                 }
1035                                 result = result.expand();
1036                                 map_trafo_H_mult recursion;
1037                                 return recursion(result);
1038                         } else {
1039                                 return e;
1040                         }
1041
1042                 }
1043                 return e;
1044         }
1045 };
1046
1047
1048 // do integration [ReV] (49)
1049 // put parameter 1 in front of existing parameters
1050 ex trafo_H_prepend_one(const ex& e, const ex& arg)
1051 {
1052         ex h;
1053         std::string name;
1054         if (is_a<function>(e)) {
1055                 name = ex_to<function>(e).get_name();
1056         }
1057         if (name == "H") {
1058                 h = e;
1059         } else {
1060                 for (int i=0; i<e.nops(); i++) {
1061                         if (is_a<function>(e.op(i))) {
1062                                 std::string name = ex_to<function>(e.op(i)).get_name();
1063                                 if (name == "H") {
1064                                         h = e.op(i);
1065                                 }
1066                         }
1067                 }
1068         }
1069         if (h != 0) {
1070                 lst newparameter = ex_to<lst>(h.op(0));
1071                 newparameter.prepend(1);
1072                 return e.subs(h == H(newparameter, h.op(1)).hold());
1073         } else {
1074                 return e * H(lst(1),1-arg).hold();
1075         }
1076 }
1077
1078
1079 // do integration [ReV] (55)
1080 // put parameter 0 in front of existing parameters
1081 ex trafo_H_prepend_zero(const ex& e, const ex& arg)
1082 {
1083         ex h;
1084         std::string name;
1085         if (is_a<function>(e)) {
1086                 name = ex_to<function>(e).get_name();
1087         }
1088         if (name == "H") {
1089                 h = e;
1090         } else {
1091                 for (int i=0; i<e.nops(); i++) {
1092                         if (is_a<function>(e.op(i))) {
1093                                 std::string name = ex_to<function>(e.op(i)).get_name();
1094                                 if (name == "H") {
1095                                         h = e.op(i);
1096                                 }
1097                         }
1098                 }
1099         }
1100         if (h != 0) {
1101                 lst newparameter = ex_to<lst>(h.op(0));
1102                 newparameter.prepend(0);
1103                 ex addzeta = convert_from_RV(newparameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1104                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
1105         } else {
1106                 return e * (-H(lst(0),1/arg).hold());
1107         }
1108 }
1109
1110
1111 // do x -> 1-x transformation
1112 struct map_trafo_H_1mx : public map_function
1113 {
1114         ex operator()(const ex& e)
1115         {
1116                 if (is_a<add>(e) || is_a<mul>(e)) {
1117                         return e.map(*this);
1118                 }
1119                 
1120                 if (is_a<function>(e)) {
1121                         std::string name = ex_to<function>(e).get_name();
1122                         if (name == "H") {
1123
1124                                 lst parameter = ex_to<lst>(e.op(0));
1125                                 ex arg = e.op(1);
1126
1127                                 // if all parameters are either zero or one return the transformed function
1128                                 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1129                                         lst newparameter;
1130                                         for (int i=parameter.nops(); i>0; i--) {
1131                                                 newparameter.append(0);
1132                                         }
1133                                         return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1134                                 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1135                                         lst newparameter;
1136                                         for (int i=parameter.nops(); i>0; i--) {
1137                                                 newparameter.append(1);
1138                                         }
1139                                         return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1140                                 }
1141
1142                                 lst newparameter = parameter;
1143                                 newparameter.remove_first();
1144
1145                                 if (parameter.op(0) == 0) {
1146
1147                                         // leading zero
1148                                         ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1149                                         map_trafo_H_1mx recursion;
1150                                         ex buffer = recursion(H(newparameter, arg).hold());
1151                                         if (is_a<add>(buffer)) {
1152                                                 for (int i=0; i<buffer.nops(); i++) {
1153                                                         res -= trafo_H_prepend_one(buffer.op(i), arg);
1154                                                 }
1155                                         } else {
1156                                                 res -= trafo_H_prepend_one(buffer, arg);
1157                                         }
1158                                         return res;
1159
1160                                 } else {
1161
1162                                         // leading one
1163                                         map_trafo_H_1mx recursion;
1164                                         map_trafo_H_mult unify;
1165                                         ex res;
1166                                         int firstzero = 0;
1167                                         while (parameter.op(firstzero) == 1) {
1168                                                 firstzero++;
1169                                         }
1170                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1171                                                 lst newparameter;
1172                                                 int j=0;
1173                                                 for (; j<=i; j++) {
1174                                                         newparameter.append(parameter[j+1]);
1175                                                 }
1176                                                 newparameter.append(1);
1177                                                 for (; j<parameter.nops()-1; j++) {
1178                                                         newparameter.append(parameter[j+1]);
1179                                                 }
1180                                                 res -= H(newparameter, arg).hold();
1181                                         }
1182                                         return (unify((-H(lst(0), 1-arg).hold() * recursion(H(newparameter, arg).hold())).expand()) +
1183                                                         recursion(res)) / firstzero;
1184
1185                                 }
1186
1187                         }
1188                 }
1189                 return e;
1190         }
1191 };
1192
1193
1194 // do x -> 1/x transformation
1195 struct map_trafo_H_1overx : public map_function
1196 {
1197         ex operator()(const ex& e)
1198         {
1199                 if (is_a<add>(e) || is_a<mul>(e)) {
1200                         return e.map(*this);
1201                 }
1202
1203                 if (is_a<function>(e)) {
1204                         std::string name = ex_to<function>(e).get_name();
1205                         if (name == "H") {
1206
1207                                 lst parameter = ex_to<lst>(e.op(0));
1208                                 ex arg = e.op(1);
1209
1210                                 // if all parameters are either zero or one return the transformed function
1211                                 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1212                                         map_trafo_H_mult unify;
1213                                         return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() - I*Pi, parameter.nops()) / 
1214                                                                 factorial(parameter.nops())).expand());
1215                                 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1216                                         return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
1217                                 }
1218
1219                                 lst newparameter = parameter;
1220                                 newparameter.remove_first();
1221
1222                                 if (parameter.op(0) == 0) {
1223                                         
1224                                         // leading zero
1225                                         ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1226                                         map_trafo_H_1overx recursion;
1227                                         ex buffer = recursion(H(newparameter, arg).hold());
1228                                         if (is_a<add>(buffer)) {
1229                                                 for (int i=0; i<buffer.nops(); i++) {
1230                                                         res += trafo_H_prepend_zero(buffer.op(i), arg);
1231                                                 }
1232                                         } else {
1233                                                 res += trafo_H_prepend_zero(buffer, arg);
1234                                         }
1235                                         return res;
1236
1237                                 } else {
1238
1239                                         // leading one
1240                                         map_trafo_H_1overx recursion;
1241                                         map_trafo_H_mult unify;
1242                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
1243                                         int firstzero = 0;
1244                                         while (parameter.op(firstzero) == 1) {
1245                                                 firstzero++;
1246                                         }
1247                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1248                                                 lst newparameter;
1249                                                 int j=0;
1250                                                 for (; j<=i; j++) {
1251                                                         newparameter.append(parameter[j+1]);
1252                                                 }
1253                                                 newparameter.append(1);
1254                                                 for (; j<parameter.nops()-1; j++) {
1255                                                         newparameter.append(parameter[j+1]);
1256                                                 }
1257                                                 res -= H(newparameter, arg).hold();
1258                                         }
1259                                         res = recursion(res).expand() / firstzero;
1260                                         return unify(res);
1261
1262                                 }
1263
1264                         }
1265                 }
1266                 return e;
1267         }
1268 };
1269
1270
1271 // remove trailing zeros from H-parameters
1272 struct map_trafo_H_reduce_trailing_zeros : public map_function
1273 {
1274         ex operator()(const ex& e)
1275         {
1276                 if (is_a<add>(e) || is_a<mul>(e)) {
1277                         return e.map(*this);
1278                 }
1279                 if (is_a<function>(e)) {
1280                         std::string name = ex_to<function>(e).get_name();
1281                         if (name == "H") {
1282                                 lst parameter;
1283                                 if (is_a<lst>(e.op(0))) {
1284                                                 parameter = ex_to<lst>(e.op(0));
1285                                 } else {
1286                                         parameter = lst(e.op(0));
1287                                 }
1288                                 ex arg = e.op(1);
1289                                 if (parameter.op(parameter.nops()-1) == 0) {
1290                                         
1291                                         //
1292                                         if (parameter.nops() == 1) {
1293                                                 return log(arg);
1294                                         }
1295                                         
1296                                         //
1297                                         lst::const_iterator it = parameter.begin();
1298                                         while ((it != parameter.end()) && (*it == 0)) {
1299                                                 it++;
1300                                         }
1301                                         if (it == parameter.end()) {
1302                                                 return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
1303                                         }
1304                                         
1305                                         //
1306                                         parameter.remove_last();
1307                                         int lastentry = parameter.nops();
1308                                         while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
1309                                                 lastentry--;
1310                                         }
1311                                         
1312                                         //
1313                                         ex result = log(arg) * H(parameter,arg).hold();
1314                                         for (ex i=0; i<lastentry; i++) {
1315                                                 parameter[i]++;
1316                                                 result -= (parameter[i]-1) * H(parameter, arg).hold();
1317                                                 parameter[i]--;
1318                                         }
1319                                         
1320                                         if (lastentry < parameter.nops()) {
1321                                                 result = result / (parameter.nops()-lastentry+1);
1322                                                 return result.map(*this);
1323                                         } else {
1324                                                 return result;
1325                                         }
1326                                 }
1327                         }
1328                 }
1329                 return e;
1330         }
1331 };
1332
1333
1334 // recursively call convert_from_RV on expression
1335 struct map_trafo_H_convert : public map_function
1336 {
1337         ex operator()(const ex& e)
1338         {
1339                 if (is_a<add>(e) || is_a<mul>(e) || is_a<power>(e)) {
1340                         return e.map(*this);
1341                 }
1342                 if (is_a<function>(e)) {
1343                         std::string name = ex_to<function>(e).get_name();
1344                         if (name == "H") {
1345                                 lst parameter = ex_to<lst>(e.op(0));
1346                                 ex arg = e.op(1);
1347                                 return convert_from_RV(parameter, arg);
1348                         }
1349                 }
1350                 return e;
1351         }
1352 };
1353
1354
1355 // translate notation from nested sums to Remiddi/Vermaseren
1356 lst convert_to_RV(const lst& o)
1357 {
1358         lst res;
1359         for (lst::const_iterator it = o.begin(); it != o.end(); it++) {
1360                 for (ex i=0; i<(*it)-1; i++) {
1361                         res.append(0);
1362                 }
1363                 res.append(1);
1364         }
1365         return res;
1366 }
1367
1368
1369 // translate notation from Remiddi/Vermaseren to nested sums
1370 ex convert_from_RV(const lst& parameterlst, const ex& arg)
1371 {
1372         lst newparameterlst;
1373
1374         lst::const_iterator it = parameterlst.begin();
1375         int count = 1;
1376         while (it != parameterlst.end()) {
1377                 if (*it == 0) {
1378                         count++;
1379                 } else {
1380                         newparameterlst.append((*it>0) ? count : -count);
1381                         count = 1;
1382                 }
1383                 it++;
1384         }
1385         for (int i=1; i<count; i++) {
1386                 newparameterlst.append(0);
1387         }
1388         
1389         map_trafo_H_reduce_trailing_zeros filter;
1390         return filter(H(newparameterlst, arg).hold());
1391 }
1392
1393
1394 // do the actual summation.
1395 cln::cl_N H_do_sum(const std::vector<int>& s, const cln::cl_N& x)
1396 {
1397         const int j = s.size();
1398
1399         std::vector<cln::cl_N> t(j);
1400
1401         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1402         cln::cl_N factor = cln::expt(x, j) * one;
1403         cln::cl_N t0buf;
1404         int q = 0;
1405         do {
1406                 t0buf = t[0];
1407                 q++;
1408                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1409                 for (int k=j-2; k>=1; k--) {
1410                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1411                 }
1412                 t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), s[0]);
1413                 factor = factor * x;
1414         } while (t[0] != t0buf);
1415         
1416         return t[0];
1417 }
1418
1419
1420 } // end of anonymous namespace
1421
1422
1423 //////////////////////////////////////////////////////////////////////
1424 //
1425 // Harmonic polylogarithm  H(m,x)
1426 //
1427 // GiNaC function
1428 //
1429 //////////////////////////////////////////////////////////////////////
1430
1431
1432 static ex H_eval(const ex& x1, const ex& x2)
1433 {
1434         if (x2 == 0) {
1435                 return 0;
1436         }
1437         if (x2 == 1) {
1438                 return zeta(x1);
1439         }
1440         if (x1.nops() == 1) {
1441                 return Li(x1.op(0), x2);
1442         }
1443         if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational))) {
1444                 return H(x1,x2).evalf();
1445         }
1446         return H(x1,x2).hold();
1447 }
1448
1449
1450 static ex H_evalf(const ex& x1, const ex& x2)
1451 {
1452         if (is_a<lst>(x1) && is_a<numeric>(x2)) {
1453                 for (int i=0; i<x1.nops(); i++) {
1454                         if (!x1.op(i).info(info_flags::posint)) {
1455                                 return H(x1,x2).hold();
1456                         }
1457                 }
1458                 if (x1.nops() < 1) {
1459                         return _ex1;
1460                 }
1461                 if (x1.nops() == 1) {
1462                         return Li(x1.op(0), x2).evalf();
1463                 }
1464                 cln::cl_N x = ex_to<numeric>(x2).to_cl_N();
1465                 if (x == 1) {
1466                         return zeta(x1).evalf();
1467                 }
1468
1469                 // choose trafo
1470                 if (cln::abs(x) > 1) {
1471                         symbol xtemp("xtemp");
1472                         map_trafo_H_1overx trafo;
1473                         ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1474                         map_trafo_H_convert converter;
1475                         res = converter(res);
1476                         return res.subs(xtemp==x2).evalf();
1477                 }
1478
1479                 // since the x->1-x transformation produces a lot of terms, it is only
1480                 // efficient for argument near one.
1481                 if (cln::realpart(x) > 0.95) {
1482                         symbol xtemp("xtemp");
1483                         map_trafo_H_1mx trafo;
1484                         ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1485                         map_trafo_H_convert converter;
1486                         res = converter(res);
1487                         return res.subs(xtemp==x2).evalf();
1488                 }
1489
1490                 // no trafo -> do summation
1491                 int count = x1.nops();
1492                 std::vector<int> r(count);
1493                 for (int i=0; i<count; i++) {
1494                         r[i] = ex_to<numeric>(x1.op(i)).to_int();
1495                 }
1496
1497                 return numeric(H_do_sum(r,x));
1498         }
1499
1500         return H(x1,x2).hold();
1501 }
1502
1503
1504 static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
1505 {
1506         epvector seq;
1507         seq.push_back(expair(H(x1,x2), 0));
1508         return pseries(rel,seq);
1509 }
1510
1511
1512 static ex H_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
1513 {
1514         GINAC_ASSERT(deriv_param < 2);
1515         if (deriv_param == 0) {
1516                 return _ex0;
1517         }
1518         if (is_a<lst>(x1)) {
1519                 lst newparameter = ex_to<lst>(x1);
1520                 if (x1.op(0) == 1) {
1521                         newparameter.remove_first();
1522                         return 1/(1-x2) * H(newparameter, x2);
1523                 } else {
1524                         newparameter[0]--;
1525                         return H(newparameter, x2).hold() / x2;
1526                 }
1527         } else {
1528                 if (x1 == 1) {
1529                         return 1/(1-x2);
1530                 } else {
1531                         return H(x1-1, x2).hold() / x2;
1532                 }
1533         }
1534 }
1535
1536
1537 REGISTER_FUNCTION(H,
1538                 eval_func(H_eval).
1539                 evalf_func(H_evalf).
1540                 do_not_evalf_params().
1541                 series_func(H_series).
1542                 derivative_func(H_deriv));
1543
1544
1545 //////////////////////////////////////////////////////////////////////
1546 //
1547 // Multiple zeta values  zeta(x) and zeta(x,s)
1548 //
1549 // helper functions
1550 //
1551 //////////////////////////////////////////////////////////////////////
1552
1553
1554 // anonymous namespace for helper functions
1555 namespace {
1556
1557
1558 // parameters and data for [Cra] algorithm
1559 const cln::cl_N lambda = cln::cl_N("319/320");
1560 int L1;
1561 int L2;
1562 std::vector<std::vector<cln::cl_N> > f_kj;
1563 std::vector<cln::cl_N> crB;
1564 std::vector<std::vector<cln::cl_N> > crG;
1565 std::vector<cln::cl_N> crX;
1566
1567
1568 void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
1569 {
1570         const int size = a.size();
1571         for (int n=0; n<size; n++) {
1572                 c[n] = 0;
1573                 for (int m=0; m<=n; m++) {
1574                         c[n] = c[n] + a[m]*b[n-m];
1575                 }
1576         }
1577 }
1578
1579
1580 // [Cra] section 4
1581 void initcX(const std::vector<int>& s)
1582 {
1583         const int k = s.size();
1584
1585         crX.clear();
1586         crG.clear();
1587         crB.clear();
1588
1589         for (int i=0; i<=L2; i++) {
1590                 crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
1591         }
1592
1593         int Sm = 0;
1594         int Smp1 = 0;
1595         for (int m=0; m<k-1; m++) {
1596                 std::vector<cln::cl_N> crGbuf;
1597                 Sm = Sm + s[m];
1598                 Smp1 = Sm + s[m+1];
1599                 for (int i=0; i<=L2; i++) {
1600                         crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
1601                 }
1602                 crG.push_back(crGbuf);
1603         }
1604
1605         crX = crB;
1606
1607         for (int m=0; m<k-1; m++) {
1608                 std::vector<cln::cl_N> Xbuf;
1609                 for (int i=0; i<=L2; i++) {
1610                         Xbuf.push_back(crX[i] * crG[m][i]);
1611                 }
1612                 halfcyclic_convolute(Xbuf, crB, crX);
1613         }
1614 }
1615
1616
1617 // [Cra] section 4
1618 cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
1619 {
1620         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1621         cln::cl_N factor = cln::expt(lambda, Sqk);
1622         cln::cl_N res = factor / Sqk * crX[0] * one;
1623         cln::cl_N resbuf;
1624         int N = 0;
1625         do {
1626                 resbuf = res;
1627                 factor = factor * lambda;
1628                 N++;
1629                 res = res + crX[N] * factor / (N+Sqk);
1630         } while ((res != resbuf) || cln::zerop(crX[N]));
1631         return res;
1632 }
1633
1634
1635 // [Cra] section 4
1636 void calc_f(int maxr)
1637 {
1638         f_kj.clear();
1639         f_kj.resize(L1);
1640         
1641         cln::cl_N t0, t1, t2, t3, t4;
1642         int i, j, k;
1643         std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
1644         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1645         
1646         t0 = cln::exp(-lambda);
1647         t2 = 1;
1648         for (k=1; k<=L1; k++) {
1649                 t1 = k * lambda;
1650                 t2 = t0 * t2;
1651                 for (j=1; j<=maxr; j++) {
1652                         t3 = 1;
1653                         t4 = 1;
1654                         for (i=2; i<=j; i++) {
1655                                 t4 = t4 * (j-i+1);
1656                                 t3 = t1 * t3 + t4;
1657                         }
1658                         (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one);
1659                 }
1660                 it++;
1661         }
1662 }
1663
1664
1665 // [Cra] (3.1)
1666 cln::cl_N crandall_Z(const std::vector<int>& s)
1667 {
1668         const int j = s.size();
1669
1670         if (j == 1) {   
1671                 cln::cl_N t0;
1672                 cln::cl_N t0buf;
1673                 int q = 0;
1674                 do {
1675                         t0buf = t0;
1676                         q++;
1677                         t0 = t0 + f_kj[q+j-2][s[0]-1];
1678                 } while (t0 != t0buf);
1679                 
1680                 return t0 / cln::factorial(s[0]-1);
1681         }
1682
1683         std::vector<cln::cl_N> t(j);
1684
1685         cln::cl_N t0buf;
1686         int q = 0;
1687         do {
1688                 t0buf = t[0];
1689                 q++;
1690                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1691                 for (int k=j-2; k>=1; k--) {
1692                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1693                 }
1694                 t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
1695         } while (t[0] != t0buf);
1696         
1697         return t[0] / cln::factorial(s[0]-1);
1698 }
1699
1700
1701 // [Cra] (2.4)
1702 cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
1703 {
1704         std::vector<int> r = s;
1705         const int j = r.size();
1706
1707         // decide on maximal size of f_kj for crandall_Z
1708         if (Digits < 50) {
1709                 L1 = 150;
1710         } else {
1711                 L1 = Digits * 3 + j*2;
1712         }
1713
1714         // decide on maximal size of crX for crandall_Y
1715         if (Digits < 38) {
1716                 L2 = 63;
1717         } else if (Digits < 86) {
1718                 L2 = 127;
1719         } else if (Digits < 192) {
1720                 L2 = 255;
1721         } else if (Digits < 394) {
1722                 L2 = 511;
1723         } else if (Digits < 808) {
1724                 L2 = 1023;
1725         } else {
1726                 L2 = 2047;
1727         }
1728
1729         cln::cl_N res;
1730
1731         int maxr = 0;
1732         int S = 0;
1733         for (int i=0; i<j; i++) {
1734                 S += r[i];
1735                 if (r[i] > maxr) {
1736                         maxr = r[i];
1737                 }
1738         }
1739
1740         calc_f(maxr);
1741
1742         const cln::cl_N r0factorial = cln::factorial(r[0]-1);
1743
1744         std::vector<int> rz;
1745         int skp1buf;
1746         int Srun = S;
1747         for (int k=r.size()-1; k>0; k--) {
1748
1749                 rz.insert(rz.begin(), r.back());
1750                 skp1buf = rz.front();
1751                 Srun -= skp1buf;
1752                 r.pop_back();
1753
1754                 initcX(r);
1755                 
1756                 for (int q=0; q<skp1buf; q++) {
1757                         
1758                         cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
1759                         cln::cl_N pp2 = crandall_Z(rz);
1760
1761                         rz.front()--;
1762                         
1763                         if (q & 1) {
1764                                 res = res - pp1 * pp2 / cln::factorial(q);
1765                         } else {
1766                                 res = res + pp1 * pp2 / cln::factorial(q);
1767                         }
1768                 }
1769                 rz.front() = skp1buf;
1770         }
1771         rz.insert(rz.begin(), r.back());
1772
1773         initcX(rz);
1774
1775         res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
1776
1777         return res;
1778 }
1779
1780
1781 cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
1782 {
1783         const int j = r.size();
1784
1785         // buffer for subsums
1786         std::vector<cln::cl_N> t(j);
1787         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1788
1789         cln::cl_N t0buf;
1790         int q = 0;
1791         do {
1792                 t0buf = t[0];
1793                 q++;
1794                 t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
1795                 for (int k=j-2; k>=0; k--) {
1796                         t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
1797                 }
1798         } while (t[0] != t0buf);
1799
1800         return t[0];
1801 }
1802
1803
1804 // does Hoelder convolution. see [BBB] (7.0)
1805 cln::cl_N zeta_do_Hoelder_convolution(const std::vector<int>& m_, const std::vector<int>& s_)
1806 {
1807         // prepare parameters
1808         // holds Li arguments in [BBB] notation
1809         std::vector<int> s = s_;
1810         std::vector<int> m_p = m_;
1811         std::vector<int> m_q;
1812         // holds Li arguments in nested sums notation
1813         std::vector<cln::cl_N> s_p(s.size(), cln::cl_N(1));
1814         s_p[0] = s_p[0] * cln::cl_N("1/2");
1815         // convert notations
1816         int sig = 1;
1817         for (int i=0; i<s_.size(); i++) {
1818                 if (s_[i] < 0) {
1819                         sig = -sig;
1820                         s_p[i] = -s_p[i];
1821                 }
1822                 s[i] = sig * std::abs(s[i]);
1823         }
1824         std::vector<cln::cl_N> s_q;
1825         cln::cl_N signum = 1;
1826
1827         // first term
1828         cln::cl_N res = multipleLi_do_sum(m_p, s_p);
1829
1830         // middle terms
1831         do {
1832
1833                 // change parameters
1834                 if (s.front() > 0) {
1835                         if (m_p.front() == 1) {
1836                                 m_p.erase(m_p.begin());
1837                                 s_p.erase(s_p.begin());
1838                                 if (s_p.size() > 0) {
1839                                         s_p.front() = s_p.front() * cln::cl_N("1/2");
1840                                 }
1841                                 s.erase(s.begin());
1842                                 m_q.front()++;
1843                         } else {
1844                                 m_p.front()--;
1845                                 m_q.insert(m_q.begin(), 1);
1846                                 if (s_q.size() > 0) {
1847                                         s_q.front() = s_q.front() * 2;
1848                                 }
1849                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
1850                         }
1851                 } else {
1852                         if (m_p.front() == 1) {
1853                                 m_p.erase(m_p.begin());
1854                                 s_p.erase(s_p.begin());
1855                                 if (s_p.size() > 0) {
1856                                         s_p.front() = s_p.front() * cln::cl_N("1/2");
1857                                 }
1858                                 s.erase(s.begin());
1859                                 for (int i=0; i<s.size(); i++) {
1860                                         s_p[i] = -s_p[i];
1861                                         if (s[i] > 0) break;
1862                                         
1863                                 }
1864                                 m_q.insert(m_q.begin(), 1);
1865                                 if (s_q.size() > 0) {
1866                                         s_q.front() = s_q.front() * 4;
1867                                 }
1868                                 s_q.insert(s_q.begin(), cln::cl_N("1/4"));
1869                                 signum = -signum;
1870                         } else {
1871                                 m_p.front()--;
1872                                 m_q.insert(m_q.begin(), 1);
1873                                 if (s_q.size() > 0) {
1874                                         s_q.front() = s_q.front() * 2;
1875                                 }
1876                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
1877                         }
1878                 }
1879
1880                 // exiting the loop
1881                 if (m_p.size() == 0) break;
1882
1883                 res = res + signum * multipleLi_do_sum(m_p, s_p) * multipleLi_do_sum(m_q, s_q);
1884                 
1885         } while (true);
1886
1887         // last term
1888         res = res + signum * multipleLi_do_sum(m_q, s_q);
1889         
1890         return res;
1891 }
1892
1893
1894 } // end of anonymous namespace
1895
1896
1897 //////////////////////////////////////////////////////////////////////
1898 //
1899 // Multiple zeta values  zeta(x)
1900 //
1901 // GiNaC function
1902 //
1903 //////////////////////////////////////////////////////////////////////
1904
1905
1906 static ex zeta1_evalf(const ex& x)
1907 {
1908         if (is_exactly_a<lst>(x) && (x.nops()>1)) {
1909
1910                 // multiple zeta value
1911                 const int count = x.nops();
1912                 const lst& xlst = ex_to<lst>(x);
1913                 std::vector<int> r(count);
1914
1915                 // check parameters and convert them
1916                 lst::const_iterator it1 = xlst.begin();
1917                 std::vector<int>::iterator it2 = r.begin();
1918                 do {
1919                         if (!(*it1).info(info_flags::posint)) {
1920                                 return zeta(x).hold();
1921                         }
1922                         *it2 = ex_to<numeric>(*it1).to_int();
1923                         it1++;
1924                         it2++;
1925                 } while (it2 != r.end());
1926
1927                 // check for divergence
1928                 if (r[0] == 1) {
1929                         return zeta(x).hold();
1930                 }
1931
1932                 // decide on summation algorithm
1933                 // this is still a bit clumsy
1934                 int limit = (Digits>17) ? 10 : 6;
1935                 if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) {
1936                         return numeric(zeta_do_sum_Crandall(r));
1937                 } else {
1938                         return numeric(zeta_do_sum_simple(r));
1939                 }
1940         }
1941                 
1942         // single zeta value
1943         if (is_exactly_a<numeric>(x) && (x != 1)) {
1944                 try {
1945                         return zeta(ex_to<numeric>(x));
1946                 } catch (const dunno &e) { }
1947         }
1948
1949         return zeta(x).hold();
1950 }
1951
1952
1953 static ex zeta1_eval(const ex& x)
1954 {
1955         if (is_exactly_a<lst>(x)) {
1956                 if (x.nops() == 1) {
1957                         return zeta(x.op(0));
1958                 }
1959                 return zeta(x).hold();
1960         }
1961
1962         if (x.info(info_flags::numeric)) {
1963                 const numeric& y = ex_to<numeric>(x);
1964                 // trap integer arguments:
1965                 if (y.is_integer()) {
1966                         if (y.is_zero()) {
1967                                 return _ex_1_2;
1968                         }
1969                         if (y.is_equal(_num1)) {
1970                                 return zeta(x).hold();
1971                         }
1972                         if (y.info(info_flags::posint)) {
1973                                 if (y.info(info_flags::odd)) {
1974                                         return zeta(x).hold();
1975                                 } else {
1976                                         return abs(bernoulli(y)) * pow(Pi, y) * pow(_num2, y-_num1) / factorial(y);
1977                                 }
1978                         } else {
1979                                 if (y.info(info_flags::odd)) {
1980                                         return -bernoulli(_num1-y) / (_num1-y);
1981                                 } else {
1982                                         return _ex0;
1983                                 }
1984                         }
1985                 }
1986                 // zeta(float)
1987                 if (y.info(info_flags::numeric) && !y.info(info_flags::crational))
1988                         return zeta1_evalf(x);
1989         }
1990         return zeta(x).hold();
1991 }
1992
1993
1994 static ex zeta1_deriv(const ex& x, unsigned deriv_param)
1995 {
1996         GINAC_ASSERT(deriv_param==0);
1997
1998         if (is_exactly_a<lst>(x)) {
1999                 return _ex0;
2000         } else {
2001                 return zeta(_ex1, x);
2002         }
2003 }
2004
2005
2006 unsigned zeta1_SERIAL::serial =
2007                         function::register_new(function_options("zeta").
2008                                                 eval_func(zeta1_eval).
2009                                                 evalf_func(zeta1_evalf).
2010                                                 do_not_evalf_params().
2011                                                 derivative_func(zeta1_deriv).
2012                                                 latex_name("\\zeta").
2013                                                 overloaded(2));
2014
2015
2016 //////////////////////////////////////////////////////////////////////
2017 //
2018 // Alternating Euler sum  zeta(x,s)
2019 //
2020 // GiNaC function
2021 //
2022 //////////////////////////////////////////////////////////////////////
2023
2024
2025 static ex zeta2_evalf(const ex& x, const ex& s)
2026 {
2027         if (is_exactly_a<lst>(x)) {
2028
2029                 // alternating Euler sum
2030                 const int count = x.nops();
2031                 const lst& xlst = ex_to<lst>(x);
2032                 const lst& slst = ex_to<lst>(s);
2033                 std::vector<int> xi(count);
2034                 std::vector<int> si(count);
2035
2036                 // check parameters and convert them
2037                 lst::const_iterator it_xread = xlst.begin();
2038                 lst::const_iterator it_sread = slst.begin();
2039                 std::vector<int>::iterator it_xwrite = xi.begin();
2040                 std::vector<int>::iterator it_swrite = si.begin();
2041                 do {
2042                         if (!(*it_xread).info(info_flags::posint)) {
2043                                 return zeta(x, s).hold();
2044                         }
2045                         *it_xwrite = ex_to<numeric>(*it_xread).to_int();
2046                         if (*it_sread > 0) {
2047                                 *it_swrite = 1;
2048                         } else {
2049                                 *it_swrite = -1;
2050                         }
2051                         it_xread++;
2052                         it_sread++;
2053                         it_xwrite++;
2054                         it_swrite++;
2055                 } while (it_xwrite != xi.end());
2056
2057                 // check for divergence
2058                 if ((xi[0] == 1) && (si[0] == 1)) {
2059                         return zeta(x, s).hold();
2060                 }
2061
2062                 // use Hoelder convolution
2063                 return numeric(zeta_do_Hoelder_convolution(xi, si));
2064         }
2065                 
2066         return zeta(x, s).hold();
2067 }
2068
2069
2070 static ex zeta2_eval(const ex& x, const ex& s)
2071 {
2072         if (is_exactly_a<lst>(s)) {
2073                 const lst& l = ex_to<lst>(s);
2074                 lst::const_iterator it = l.begin();
2075                 while (it != l.end()) {
2076                         if ((*it).info(info_flags::negative)) {
2077                                 return zeta(x, s).hold();
2078                         }
2079                         it++;
2080                 }
2081                 return zeta(x);
2082         } else {
2083                 if (s.info(info_flags::positive)) {
2084                         return zeta(x);
2085                 }
2086         }
2087
2088         return zeta(x, s).hold();
2089 }
2090
2091
2092 static ex zeta2_deriv(const ex& x, const ex& s, unsigned deriv_param)
2093 {
2094         GINAC_ASSERT(deriv_param==0);
2095
2096         if (is_exactly_a<lst>(x)) {
2097                 return _ex0;
2098         } else {
2099                 if ((is_exactly_a<lst>(s) && (s.op(0) > 0)) || (s > 0)) {
2100                         return zeta(_ex1, x);
2101                 }
2102                 return _ex0;
2103         }
2104 }
2105
2106
2107 unsigned zeta2_SERIAL::serial =
2108                         function::register_new(function_options("zeta").
2109                                                 eval_func(zeta2_eval).
2110                                                 evalf_func(zeta2_evalf).
2111                                                 do_not_evalf_params().
2112                                                 derivative_func(zeta2_deriv).
2113                                                 latex_name("\\zeta").
2114                                                 overloaded(2));
2115
2116
2117 } // namespace GiNaC
2118