* Fixed bug in zeta(m,s).
[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                 // do it once ...
411                 q++;
412                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
413                 for (int k=j-2; k>=0; k--) {
414                         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]);
415                 }
416                 // ... and do it again (to avoid premature drop out due to special arguments)
417                 q++;
418                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
419                 for (int k=j-2; k>=0; k--) {
420                         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]);
421                 }
422         } while (t[0] != t0buf);
423         
424         return t[0];
425 }
426
427
428 } // end of anonymous namespace
429
430
431 //////////////////////////////////////////////////////////////////////
432 //
433 // Classical polylogarithm and multiple polylogarithm  Li(n,x)
434 //
435 // GiNaC function
436 //
437 //////////////////////////////////////////////////////////////////////
438
439
440 static ex Li_eval(const ex& x1, const ex& x2)
441 {
442         if (x2.is_zero()) {
443                 return _ex0;
444         }
445         else {
446                 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational)))
447                         return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
448                 if (is_a<lst>(x2)) {
449                         for (int i=0; i<x2.nops(); i++) {
450                                 if (!is_a<numeric>(x2.op(i))) {
451                                         return Li(x1,x2).hold();
452                                 }
453                         }
454                         return Li(x1,x2).evalf();
455                 }
456                 return Li(x1,x2).hold();
457         }
458 }
459
460
461 static ex Li_evalf(const ex& x1, const ex& x2)
462 {
463         // classical polylogs
464         if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
465                 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
466         }
467         // multiple polylogs
468         else if (is_a<lst>(x1) && is_a<lst>(x2)) {
469                 ex conv = 1;
470                 for (int i=0; i<x1.nops(); i++) {
471                         if (!x1.op(i).info(info_flags::posint)) {
472                                 return Li(x1,x2).hold();
473                         }
474                         if (!is_a<numeric>(x2.op(i))) {
475                                 return Li(x1,x2).hold();
476                         }
477                         conv *= x2.op(i);
478                         if ((conv > 1) || ((conv == 1) && (x1.op(0) == 1))) {
479                                 return Li(x1,x2).hold();
480                         }
481                 }
482
483                 std::vector<int> m;
484                 std::vector<cln::cl_N> x;
485                 for (int i=0; i<ex_to<numeric>(x1.nops()).to_int(); i++) {
486                         m.push_back(ex_to<numeric>(x1.op(i)).to_int());
487                         x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
488                 }
489
490                 return numeric(multipleLi_do_sum(m, x));
491         }
492
493         return Li(x1,x2).hold();
494 }
495
496
497 static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
498 {
499         epvector seq;
500         seq.push_back(expair(Li(x1,x2), 0));
501         return pseries(rel,seq);
502 }
503
504
505 static ex Li_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
506 {
507         GINAC_ASSERT(deriv_param < 2);
508         if (deriv_param == 0) {
509                 return _ex0;
510         }
511         if (x1 > 0) {
512                 return Li(x1-1, x2) / x2;
513         } else {
514                 return 1/(1-x2);
515         }
516 }
517
518
519 REGISTER_FUNCTION(Li,
520                 eval_func(Li_eval).
521                 evalf_func(Li_evalf).
522                 do_not_evalf_params().
523                 series_func(Li_series).
524                 derivative_func(Li_deriv));
525
526
527 //////////////////////////////////////////////////////////////////////
528 //
529 // Nielsen's generalized polylogarithm  S(n,p,x)
530 //
531 // helper functions
532 //
533 //////////////////////////////////////////////////////////////////////
534
535
536 // anonymous namespace for helper functions
537 namespace {
538
539
540 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
541 // see fill_Yn()
542 std::vector<std::vector<cln::cl_N> > Yn;
543 int ynsize = 0; // number of Yn[]
544 int ynlength = 100; // initial length of all Yn[i]
545
546
547 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
548 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
549 // representing S_{n,p}(x).
550 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
551 // equivalent Z-sum.
552 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
553 // representing S_{n,p}(x).
554 // The calculation of Y_n uses the values from Y_{n-1}.
555 void fill_Yn(int n, const cln::float_format_t& prec)
556 {
557         const int initsize = ynlength;
558         //const int initsize = initsize_Yn;
559         cln::cl_N one = cln::cl_float(1, prec);
560
561         if (n) {
562                 std::vector<cln::cl_N> buf(initsize);
563                 std::vector<cln::cl_N>::iterator it = buf.begin();
564                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
565                 *it = (*itprev) / cln::cl_N(n+1) * one;
566                 it++;
567                 itprev++;
568                 // sums with an index smaller than the depth are zero and need not to be calculated.
569                 // calculation starts with depth, which is n+2)
570                 for (int i=n+2; i<=initsize+n; i++) {
571                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
572                         it++;
573                         itprev++;
574                 }
575                 Yn.push_back(buf);
576         } else {
577                 std::vector<cln::cl_N> buf(initsize);
578                 std::vector<cln::cl_N>::iterator it = buf.begin();
579                 *it = 1 * one;
580                 it++;
581                 for (int i=2; i<=initsize; i++) {
582                         *it = *(it-1) + 1 / cln::cl_N(i) * one;
583                         it++;
584                 }
585                 Yn.push_back(buf);
586         }
587         ynsize++;
588 }
589
590
591 // make Yn longer ... 
592 void make_Yn_longer(int newsize, const cln::float_format_t& prec)
593 {
594
595         cln::cl_N one = cln::cl_float(1, prec);
596
597         Yn[0].resize(newsize);
598         std::vector<cln::cl_N>::iterator it = Yn[0].begin();
599         it += ynlength;
600         for (int i=ynlength+1; i<=newsize; i++) {
601                 *it = *(it-1) + 1 / cln::cl_N(i) * one;
602                 it++;
603         }
604
605         for (int n=1; n<ynsize; n++) {
606                 Yn[n].resize(newsize);
607                 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
608                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
609                 it += ynlength;
610                 itprev += ynlength;
611                 for (int i=ynlength+n+1; i<=newsize+n; i++) {
612                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
613                         it++;
614                         itprev++;
615                 }
616         }
617         
618         ynlength = newsize;
619 }
620
621
622 // helper function for S(n,p,x)
623 // [Kol] (7.2)
624 cln::cl_N C(int n, int p)
625 {
626         cln::cl_N result;
627
628         for (int k=0; k<p; k++) {
629                 for (int j=0; j<=(n+k-1)/2; j++) {
630                         if (k == 0) {
631                                 if (n & 1) {
632                                         if (j & 1) {
633                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
634                                         }
635                                         else {
636                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
637                                         }
638                                 }
639                         }
640                         else {
641                                 if (k & 1) {
642                                         if (j & 1) {
643                                                 result = result + cln::factorial(n+k-1)
644                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
645                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
646                                         }
647                                         else {
648                                                 result = result - cln::factorial(n+k-1)
649                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
650                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
651                                         }
652                                 }
653                                 else {
654                                         if (j & 1) {
655                                                 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()
656                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
657                                         }
658                                         else {
659                                                 result = result + cln::factorial(n+k-1)
660                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
661                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
662                                         }
663                                 }
664                         }
665                 }
666         }
667         int np = n+p;
668         if ((np-1) & 1) {
669                 if (((np)/2+n) & 1) {
670                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
671                 }
672                 else {
673                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
674                 }
675         }
676
677         return result;
678 }
679
680
681 // helper function for S(n,p,x)
682 // [Kol] remark to (9.1)
683 cln::cl_N a_k(int k)
684 {
685         cln::cl_N result;
686
687         if (k == 0) {
688                 return 1;
689         }
690
691         result = result;
692         for (int m=2; m<=k; m++) {
693                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
694         }
695
696         return -result / k;
697 }
698
699
700 // helper function for S(n,p,x)
701 // [Kol] remark to (9.1)
702 cln::cl_N b_k(int k)
703 {
704         cln::cl_N result;
705
706         if (k == 0) {
707                 return 1;
708         }
709
710         result = result;
711         for (int m=2; m<=k; m++) {
712                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
713         }
714
715         return result / k;
716 }
717
718
719 // helper function for S(n,p,x)
720 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
721 {
722         if (p==1) {
723                 return Li_projection(n+1, x, prec);
724         }
725         
726         // check if precalculated values are sufficient
727         if (p > ynsize+1) {
728                 for (int i=ynsize; i<p-1; i++) {
729                         fill_Yn(i, prec);
730                 }
731         }
732
733         // should be done otherwise
734         cln::cl_N xf = x * cln::cl_float(1, prec);
735
736         cln::cl_N res;
737         cln::cl_N resbuf;
738         cln::cl_N factor = cln::expt(xf, p);
739         int i = p;
740         do {
741                 resbuf = res;
742                 if (i-p >= ynlength) {
743                         // make Yn longer
744                         make_Yn_longer(ynlength*2, prec);
745                 }
746                 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? ...
747                 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
748                 factor = factor * xf;
749                 i++;
750         } while (res != resbuf);
751         
752         return res;
753 }
754
755
756 // helper function for S(n,p,x)
757 cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
758 {
759         // [Kol] (5.3)
760         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
761
762                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
763                         * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
764
765                 for (int s=0; s<n; s++) {
766                         cln::cl_N res2;
767                         for (int r=0; r<p; r++) {
768                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
769                                         * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
770                         }
771                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
772                 }
773
774                 return result;
775         }
776         
777         return S_do_sum(n, p, x, prec);
778 }
779
780
781 // helper function for S(n,p,x)
782 numeric S_num(int n, int p, const numeric& x)
783 {
784         if (x == 1) {
785                 if (n == 1) {
786                     // [Kol] (2.22) with (2.21)
787                         return cln::zeta(p+1);
788                 }
789
790                 if (p == 1) {
791                     // [Kol] (2.22)
792                         return cln::zeta(n+1);
793                 }
794
795                 // [Kol] (9.1)
796                 cln::cl_N result;
797                 for (int nu=0; nu<n; nu++) {
798                         for (int rho=0; rho<=p; rho++) {
799                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
800                                         * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
801                         }
802                 }
803                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
804
805                 return result;
806         }
807         else if (x == -1) {
808                 // [Kol] (2.22)
809                 if (p == 1) {
810                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
811                 }
812 //              throw std::runtime_error("don't know how to evaluate this function!");
813         }
814
815         // what is the desired float format?
816         // first guess: default format
817         cln::float_format_t prec = cln::default_float_format;
818         const cln::cl_N value = x.to_cl_N();
819         // second guess: the argument's format
820         if (!x.real().is_rational())
821                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
822         else if (!x.imag().is_rational())
823                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
824
825
826         // [Kol] (5.3)
827         if (cln::realpart(value) < -0.5) {
828
829                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
830                         * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
831
832                 for (int s=0; s<n; s++) {
833                         cln::cl_N res2;
834                         for (int r=0; r<p; r++) {
835                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
836                                         * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
837                         }
838                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
839                 }
840
841                 return result;
842                 
843         }
844         // [Kol] (5.12)
845         if (cln::abs(value) > 1) {
846                 
847                 cln::cl_N result;
848
849                 for (int s=0; s<p; s++) {
850                         for (int r=0; r<=s; r++) {
851                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
852                                         / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
853                                         * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
854                         }
855                 }
856                 result = result * cln::expt(cln::cl_I(-1),n);
857
858                 cln::cl_N res2;
859                 for (int r=0; r<n; r++) {
860                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
861                 }
862                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
863
864                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
865
866                 return result;
867         }
868         else {
869                 return S_projection(n, p, value, prec);
870         }
871 }
872
873
874 } // end of anonymous namespace
875
876
877 //////////////////////////////////////////////////////////////////////
878 //
879 // Nielsen's generalized polylogarithm  S(n,p,x)
880 //
881 // GiNaC function
882 //
883 //////////////////////////////////////////////////////////////////////
884
885
886 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
887 {
888         if (x2 == 1) {
889                 return Li(x1+1,x3);
890         }
891         if (x3.info(info_flags::numeric) && (!x3.info(info_flags::crational)) && 
892                         x1.info(info_flags::posint) && x2.info(info_flags::posint)) {
893                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
894         }
895         return S(x1,x2,x3).hold();
896 }
897
898
899 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
900 {
901         if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
902                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
903         }
904         return S(x1,x2,x3).hold();
905 }
906
907
908 static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options)
909 {
910         epvector seq;
911         seq.push_back(expair(S(x1,x2,x3), 0));
912         return pseries(rel,seq);
913 }
914
915
916 static ex S_deriv(const ex& x1, const ex& x2, const ex& x3, unsigned deriv_param)
917 {
918         GINAC_ASSERT(deriv_param < 3);
919         if (deriv_param < 2) {
920                 return _ex0;
921         }
922         if (x1 > 0) {
923                 return S(x1-1, x2, x3) / x3;
924         } else {
925                 return S(x1, x2-1, x3) / (1-x3);
926         }
927 }
928
929
930 REGISTER_FUNCTION(S,
931                 eval_func(S_eval).
932                 evalf_func(S_evalf).
933                 do_not_evalf_params().
934                 series_func(S_series).
935                 derivative_func(S_deriv));
936
937
938 //////////////////////////////////////////////////////////////////////
939 //
940 // Harmonic polylogarithm  H(m,x)
941 //
942 // helper function
943 //
944 //////////////////////////////////////////////////////////////////////
945
946
947 // anonymous namespace for helper functions
948 namespace {
949
950
951 // forward declaration
952 ex convert_from_RV(const lst& parameterlst, const ex& arg);
953
954
955 // multiplies an one-dimensional H with another H
956 // [ReV] (18)
957 ex trafo_H_mult(const ex& h1, const ex& h2)
958 {
959         ex res;
960         ex hshort;
961         lst hlong;
962         ex h1nops = h1.op(0).nops();
963         ex h2nops = h2.op(0).nops();
964         if (h1nops > 1) {
965                 hshort = h2.op(0).op(0);
966                 hlong = ex_to<lst>(h1.op(0));
967         } else {
968                 hshort = h1.op(0).op(0);
969                 if (h2nops > 1) {
970                         hlong = ex_to<lst>(h2.op(0));
971                 } else {
972                         hlong = h2.op(0).op(0);
973                 }
974         }
975         for (int i=0; i<=hlong.nops(); i++) {
976                 lst newparameter;
977                 int j=0;
978                 for (; j<i; j++) {
979                         newparameter.append(hlong[j]);
980                 }
981                 newparameter.append(hshort);
982                 for (; j<hlong.nops(); j++) {
983                         newparameter.append(hlong[j]);
984                 }
985                 res += H(newparameter, h1.op(1)).hold();
986         }
987         return res;
988 }
989
990
991 // applies trafo_H_mult recursively on expressions
992 struct map_trafo_H_mult : public map_function
993 {
994         ex operator()(const ex& e)
995         {
996                 if (is_a<add>(e)) {
997                         return e.map(*this);
998                 }
999
1000                 if (is_a<mul>(e)) {
1001
1002                         ex result = 1;
1003                         ex firstH;
1004                         lst Hlst;
1005                         for (int pos=0; pos<e.nops(); pos++) {
1006                                 if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
1007                                         std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
1008                                         if (name == "H") {
1009                                                 for (ex i=0; i<e.op(pos).op(1); i++) {
1010                                                         Hlst.append(e.op(pos).op(0));
1011                                                 }
1012                                                 continue;
1013                                         }
1014                                 } else if (is_a<function>(e.op(pos))) {
1015                                         std::string name = ex_to<function>(e.op(pos)).get_name();
1016                                         if (name == "H") {
1017                                                 if (e.op(pos).op(0).nops() > 1) {
1018                                                         firstH = e.op(pos);
1019                                                 } else {
1020                                                         Hlst.append(e.op(pos));
1021                                                 }
1022                                                 continue;
1023                                         }
1024                                 }
1025                                 result *= e.op(pos);
1026                         }
1027                         if (firstH == 0) {
1028                                 if (Hlst.nops() > 0) {
1029                                         firstH = Hlst[Hlst.nops()-1];
1030                                         Hlst.remove_last();
1031                                 } else {
1032                                         return e;
1033                                 }
1034                         }
1035
1036                         if (Hlst.nops() > 0) {
1037                                 ex buffer = trafo_H_mult(firstH, Hlst.op(0));
1038                                 result *= buffer;
1039                                 for (int i=1; i<Hlst.nops(); i++) {
1040                                         result *= Hlst.op(i);
1041                                 }
1042                                 result = result.expand();
1043                                 map_trafo_H_mult recursion;
1044                                 return recursion(result);
1045                         } else {
1046                                 return e;
1047                         }
1048
1049                 }
1050                 return e;
1051         }
1052 };
1053
1054
1055 // do integration [ReV] (49)
1056 // put parameter 1 in front of existing parameters
1057 ex trafo_H_prepend_one(const ex& e, const ex& arg)
1058 {
1059         ex h;
1060         std::string name;
1061         if (is_a<function>(e)) {
1062                 name = ex_to<function>(e).get_name();
1063         }
1064         if (name == "H") {
1065                 h = e;
1066         } else {
1067                 for (int i=0; i<e.nops(); i++) {
1068                         if (is_a<function>(e.op(i))) {
1069                                 std::string name = ex_to<function>(e.op(i)).get_name();
1070                                 if (name == "H") {
1071                                         h = e.op(i);
1072                                 }
1073                         }
1074                 }
1075         }
1076         if (h != 0) {
1077                 lst newparameter = ex_to<lst>(h.op(0));
1078                 newparameter.prepend(1);
1079                 return e.subs(h == H(newparameter, h.op(1)).hold());
1080         } else {
1081                 return e * H(lst(1),1-arg).hold();
1082         }
1083 }
1084
1085
1086 // do integration [ReV] (55)
1087 // put parameter 0 in front of existing parameters
1088 ex trafo_H_prepend_zero(const ex& e, const ex& arg)
1089 {
1090         ex h;
1091         std::string name;
1092         if (is_a<function>(e)) {
1093                 name = ex_to<function>(e).get_name();
1094         }
1095         if (name == "H") {
1096                 h = e;
1097         } else {
1098                 for (int i=0; i<e.nops(); i++) {
1099                         if (is_a<function>(e.op(i))) {
1100                                 std::string name = ex_to<function>(e.op(i)).get_name();
1101                                 if (name == "H") {
1102                                         h = e.op(i);
1103                                 }
1104                         }
1105                 }
1106         }
1107         if (h != 0) {
1108                 lst newparameter = ex_to<lst>(h.op(0));
1109                 newparameter.prepend(0);
1110                 ex addzeta = convert_from_RV(newparameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1111                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
1112         } else {
1113                 return e * (-H(lst(0),1/arg).hold());
1114         }
1115 }
1116
1117
1118 // do x -> 1-x transformation
1119 struct map_trafo_H_1mx : public map_function
1120 {
1121         ex operator()(const ex& e)
1122         {
1123                 if (is_a<add>(e) || is_a<mul>(e)) {
1124                         return e.map(*this);
1125                 }
1126                 
1127                 if (is_a<function>(e)) {
1128                         std::string name = ex_to<function>(e).get_name();
1129                         if (name == "H") {
1130
1131                                 lst parameter = ex_to<lst>(e.op(0));
1132                                 ex arg = e.op(1);
1133
1134                                 // if all parameters are either zero or one return the transformed function
1135                                 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1136                                         lst newparameter;
1137                                         for (int i=parameter.nops(); i>0; i--) {
1138                                                 newparameter.append(0);
1139                                         }
1140                                         return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1141                                 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1142                                         lst newparameter;
1143                                         for (int i=parameter.nops(); i>0; i--) {
1144                                                 newparameter.append(1);
1145                                         }
1146                                         return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1147                                 }
1148
1149                                 lst newparameter = parameter;
1150                                 newparameter.remove_first();
1151
1152                                 if (parameter.op(0) == 0) {
1153
1154                                         // leading zero
1155                                         ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1156                                         map_trafo_H_1mx recursion;
1157                                         ex buffer = recursion(H(newparameter, arg).hold());
1158                                         if (is_a<add>(buffer)) {
1159                                                 for (int i=0; i<buffer.nops(); i++) {
1160                                                         res -= trafo_H_prepend_one(buffer.op(i), arg);
1161                                                 }
1162                                         } else {
1163                                                 res -= trafo_H_prepend_one(buffer, arg);
1164                                         }
1165                                         return res;
1166
1167                                 } else {
1168
1169                                         // leading one
1170                                         map_trafo_H_1mx recursion;
1171                                         map_trafo_H_mult unify;
1172                                         ex res;
1173                                         int firstzero = 0;
1174                                         while (parameter.op(firstzero) == 1) {
1175                                                 firstzero++;
1176                                         }
1177                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1178                                                 lst newparameter;
1179                                                 int j=0;
1180                                                 for (; j<=i; j++) {
1181                                                         newparameter.append(parameter[j+1]);
1182                                                 }
1183                                                 newparameter.append(1);
1184                                                 for (; j<parameter.nops()-1; j++) {
1185                                                         newparameter.append(parameter[j+1]);
1186                                                 }
1187                                                 res -= H(newparameter, arg).hold();
1188                                         }
1189                                         return (unify((-H(lst(0), 1-arg).hold() * recursion(H(newparameter, arg).hold())).expand()) +
1190                                                         recursion(res)) / firstzero;
1191
1192                                 }
1193
1194                         }
1195                 }
1196                 return e;
1197         }
1198 };
1199
1200
1201 // do x -> 1/x transformation
1202 struct map_trafo_H_1overx : public map_function
1203 {
1204         ex operator()(const ex& e)
1205         {
1206                 if (is_a<add>(e) || is_a<mul>(e)) {
1207                         return e.map(*this);
1208                 }
1209
1210                 if (is_a<function>(e)) {
1211                         std::string name = ex_to<function>(e).get_name();
1212                         if (name == "H") {
1213
1214                                 lst parameter = ex_to<lst>(e.op(0));
1215                                 ex arg = e.op(1);
1216
1217                                 // if all parameters are either zero or one return the transformed function
1218                                 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1219                                         map_trafo_H_mult unify;
1220                                         return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() - I*Pi, parameter.nops()) / 
1221                                                                 factorial(parameter.nops())).expand());
1222                                 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1223                                         return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
1224                                 }
1225
1226                                 lst newparameter = parameter;
1227                                 newparameter.remove_first();
1228
1229                                 if (parameter.op(0) == 0) {
1230                                         
1231                                         // leading zero
1232                                         ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1233                                         map_trafo_H_1overx recursion;
1234                                         ex buffer = recursion(H(newparameter, arg).hold());
1235                                         if (is_a<add>(buffer)) {
1236                                                 for (int i=0; i<buffer.nops(); i++) {
1237                                                         res += trafo_H_prepend_zero(buffer.op(i), arg);
1238                                                 }
1239                                         } else {
1240                                                 res += trafo_H_prepend_zero(buffer, arg);
1241                                         }
1242                                         return res;
1243
1244                                 } else {
1245
1246                                         // leading one
1247                                         map_trafo_H_1overx recursion;
1248                                         map_trafo_H_mult unify;
1249                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
1250                                         int firstzero = 0;
1251                                         while (parameter.op(firstzero) == 1) {
1252                                                 firstzero++;
1253                                         }
1254                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1255                                                 lst newparameter;
1256                                                 int j=0;
1257                                                 for (; j<=i; j++) {
1258                                                         newparameter.append(parameter[j+1]);
1259                                                 }
1260                                                 newparameter.append(1);
1261                                                 for (; j<parameter.nops()-1; j++) {
1262                                                         newparameter.append(parameter[j+1]);
1263                                                 }
1264                                                 res -= H(newparameter, arg).hold();
1265                                         }
1266                                         res = recursion(res).expand() / firstzero;
1267                                         return unify(res);
1268
1269                                 }
1270
1271                         }
1272                 }
1273                 return e;
1274         }
1275 };
1276
1277
1278 // remove trailing zeros from H-parameters
1279 struct map_trafo_H_reduce_trailing_zeros : public map_function
1280 {
1281         ex operator()(const ex& e)
1282         {
1283                 if (is_a<add>(e) || is_a<mul>(e)) {
1284                         return e.map(*this);
1285                 }
1286                 if (is_a<function>(e)) {
1287                         std::string name = ex_to<function>(e).get_name();
1288                         if (name == "H") {
1289                                 lst parameter;
1290                                 if (is_a<lst>(e.op(0))) {
1291                                                 parameter = ex_to<lst>(e.op(0));
1292                                 } else {
1293                                         parameter = lst(e.op(0));
1294                                 }
1295                                 ex arg = e.op(1);
1296                                 if (parameter.op(parameter.nops()-1) == 0) {
1297                                         
1298                                         //
1299                                         if (parameter.nops() == 1) {
1300                                                 return log(arg);
1301                                         }
1302                                         
1303                                         //
1304                                         lst::const_iterator it = parameter.begin();
1305                                         while ((it != parameter.end()) && (*it == 0)) {
1306                                                 it++;
1307                                         }
1308                                         if (it == parameter.end()) {
1309                                                 return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
1310                                         }
1311                                         
1312                                         //
1313                                         parameter.remove_last();
1314                                         int lastentry = parameter.nops();
1315                                         while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
1316                                                 lastentry--;
1317                                         }
1318                                         
1319                                         //
1320                                         ex result = log(arg) * H(parameter,arg).hold();
1321                                         for (ex i=0; i<lastentry; i++) {
1322                                                 parameter[i]++;
1323                                                 result -= (parameter[i]-1) * H(parameter, arg).hold();
1324                                                 parameter[i]--;
1325                                         }
1326                                         
1327                                         if (lastentry < parameter.nops()) {
1328                                                 result = result / (parameter.nops()-lastentry+1);
1329                                                 return result.map(*this);
1330                                         } else {
1331                                                 return result;
1332                                         }
1333                                 }
1334                         }
1335                 }
1336                 return e;
1337         }
1338 };
1339
1340
1341 // recursively call convert_from_RV on expression
1342 struct map_trafo_H_convert : public map_function
1343 {
1344         ex operator()(const ex& e)
1345         {
1346                 if (is_a<add>(e) || is_a<mul>(e) || is_a<power>(e)) {
1347                         return e.map(*this);
1348                 }
1349                 if (is_a<function>(e)) {
1350                         std::string name = ex_to<function>(e).get_name();
1351                         if (name == "H") {
1352                                 lst parameter = ex_to<lst>(e.op(0));
1353                                 ex arg = e.op(1);
1354                                 return convert_from_RV(parameter, arg);
1355                         }
1356                 }
1357                 return e;
1358         }
1359 };
1360
1361
1362 // translate notation from nested sums to Remiddi/Vermaseren
1363 lst convert_to_RV(const lst& o)
1364 {
1365         lst res;
1366         for (lst::const_iterator it = o.begin(); it != o.end(); it++) {
1367                 for (ex i=0; i<(*it)-1; i++) {
1368                         res.append(0);
1369                 }
1370                 res.append(1);
1371         }
1372         return res;
1373 }
1374
1375
1376 // translate notation from Remiddi/Vermaseren to nested sums
1377 ex convert_from_RV(const lst& parameterlst, const ex& arg)
1378 {
1379         lst newparameterlst;
1380
1381         lst::const_iterator it = parameterlst.begin();
1382         int count = 1;
1383         while (it != parameterlst.end()) {
1384                 if (*it == 0) {
1385                         count++;
1386                 } else {
1387                         newparameterlst.append((*it>0) ? count : -count);
1388                         count = 1;
1389                 }
1390                 it++;
1391         }
1392         for (int i=1; i<count; i++) {
1393                 newparameterlst.append(0);
1394         }
1395         
1396         map_trafo_H_reduce_trailing_zeros filter;
1397         return filter(H(newparameterlst, arg).hold());
1398 }
1399
1400
1401 // do the actual summation.
1402 cln::cl_N H_do_sum(const std::vector<int>& s, const cln::cl_N& x)
1403 {
1404         const int j = s.size();
1405
1406         std::vector<cln::cl_N> t(j);
1407
1408         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1409         cln::cl_N factor = cln::expt(x, j) * one;
1410         cln::cl_N t0buf;
1411         int q = 0;
1412         do {
1413                 t0buf = t[0];
1414                 q++;
1415                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1416                 for (int k=j-2; k>=1; k--) {
1417                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1418                 }
1419                 t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), s[0]);
1420                 factor = factor * x;
1421         } while (t[0] != t0buf);
1422         
1423         return t[0];
1424 }
1425
1426
1427 } // end of anonymous namespace
1428
1429
1430 //////////////////////////////////////////////////////////////////////
1431 //
1432 // Harmonic polylogarithm  H(m,x)
1433 //
1434 // GiNaC function
1435 //
1436 //////////////////////////////////////////////////////////////////////
1437
1438
1439 static ex H_eval(const ex& x1, const ex& x2)
1440 {
1441         if (x2 == 0) {
1442                 return 0;
1443         }
1444         if (x2 == 1) {
1445                 return zeta(x1);
1446         }
1447         if (x1.nops() == 1) {
1448                 return Li(x1.op(0), x2);
1449         }
1450         if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational))) {
1451                 return H(x1,x2).evalf();
1452         }
1453         return H(x1,x2).hold();
1454 }
1455
1456
1457 static ex H_evalf(const ex& x1, const ex& x2)
1458 {
1459         if (is_a<lst>(x1) && is_a<numeric>(x2)) {
1460                 for (int i=0; i<x1.nops(); i++) {
1461                         if (!x1.op(i).info(info_flags::posint)) {
1462                                 return H(x1,x2).hold();
1463                         }
1464                 }
1465                 if (x1.nops() < 1) {
1466                         return _ex1;
1467                 }
1468                 if (x1.nops() == 1) {
1469                         return Li(x1.op(0), x2).evalf();
1470                 }
1471                 cln::cl_N x = ex_to<numeric>(x2).to_cl_N();
1472                 if (x == 1) {
1473                         return zeta(x1).evalf();
1474                 }
1475
1476                 // choose trafo
1477                 if (cln::abs(x) > 1) {
1478                         symbol xtemp("xtemp");
1479                         map_trafo_H_1overx trafo;
1480                         ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1481                         map_trafo_H_convert converter;
1482                         res = converter(res);
1483                         return res.subs(xtemp==x2).evalf();
1484                 }
1485
1486                 // since the x->1-x transformation produces a lot of terms, it is only
1487                 // efficient for argument near one.
1488                 if (cln::realpart(x) > 0.95) {
1489                         symbol xtemp("xtemp");
1490                         map_trafo_H_1mx trafo;
1491                         ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1492                         map_trafo_H_convert converter;
1493                         res = converter(res);
1494                         return res.subs(xtemp==x2).evalf();
1495                 }
1496
1497                 // no trafo -> do summation
1498                 int count = x1.nops();
1499                 std::vector<int> r(count);
1500                 for (int i=0; i<count; i++) {
1501                         r[i] = ex_to<numeric>(x1.op(i)).to_int();
1502                 }
1503
1504                 return numeric(H_do_sum(r,x));
1505         }
1506
1507         return H(x1,x2).hold();
1508 }
1509
1510
1511 static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
1512 {
1513         epvector seq;
1514         seq.push_back(expair(H(x1,x2), 0));
1515         return pseries(rel,seq);
1516 }
1517
1518
1519 static ex H_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
1520 {
1521         GINAC_ASSERT(deriv_param < 2);
1522         if (deriv_param == 0) {
1523                 return _ex0;
1524         }
1525         if (is_a<lst>(x1)) {
1526                 lst newparameter = ex_to<lst>(x1);
1527                 if (x1.op(0) == 1) {
1528                         newparameter.remove_first();
1529                         return 1/(1-x2) * H(newparameter, x2);
1530                 } else {
1531                         newparameter[0]--;
1532                         return H(newparameter, x2).hold() / x2;
1533                 }
1534         } else {
1535                 if (x1 == 1) {
1536                         return 1/(1-x2);
1537                 } else {
1538                         return H(x1-1, x2).hold() / x2;
1539                 }
1540         }
1541 }
1542
1543
1544 REGISTER_FUNCTION(H,
1545                 eval_func(H_eval).
1546                 evalf_func(H_evalf).
1547                 do_not_evalf_params().
1548                 series_func(H_series).
1549                 derivative_func(H_deriv));
1550
1551
1552 //////////////////////////////////////////////////////////////////////
1553 //
1554 // Multiple zeta values  zeta(x) and zeta(x,s)
1555 //
1556 // helper functions
1557 //
1558 //////////////////////////////////////////////////////////////////////
1559
1560
1561 // anonymous namespace for helper functions
1562 namespace {
1563
1564
1565 // parameters and data for [Cra] algorithm
1566 const cln::cl_N lambda = cln::cl_N("319/320");
1567 int L1;
1568 int L2;
1569 std::vector<std::vector<cln::cl_N> > f_kj;
1570 std::vector<cln::cl_N> crB;
1571 std::vector<std::vector<cln::cl_N> > crG;
1572 std::vector<cln::cl_N> crX;
1573
1574
1575 void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
1576 {
1577         const int size = a.size();
1578         for (int n=0; n<size; n++) {
1579                 c[n] = 0;
1580                 for (int m=0; m<=n; m++) {
1581                         c[n] = c[n] + a[m]*b[n-m];
1582                 }
1583         }
1584 }
1585
1586
1587 // [Cra] section 4
1588 void initcX(const std::vector<int>& s)
1589 {
1590         const int k = s.size();
1591
1592         crX.clear();
1593         crG.clear();
1594         crB.clear();
1595
1596         for (int i=0; i<=L2; i++) {
1597                 crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
1598         }
1599
1600         int Sm = 0;
1601         int Smp1 = 0;
1602         for (int m=0; m<k-1; m++) {
1603                 std::vector<cln::cl_N> crGbuf;
1604                 Sm = Sm + s[m];
1605                 Smp1 = Sm + s[m+1];
1606                 for (int i=0; i<=L2; i++) {
1607                         crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
1608                 }
1609                 crG.push_back(crGbuf);
1610         }
1611
1612         crX = crB;
1613
1614         for (int m=0; m<k-1; m++) {
1615                 std::vector<cln::cl_N> Xbuf;
1616                 for (int i=0; i<=L2; i++) {
1617                         Xbuf.push_back(crX[i] * crG[m][i]);
1618                 }
1619                 halfcyclic_convolute(Xbuf, crB, crX);
1620         }
1621 }
1622
1623
1624 // [Cra] section 4
1625 cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
1626 {
1627         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1628         cln::cl_N factor = cln::expt(lambda, Sqk);
1629         cln::cl_N res = factor / Sqk * crX[0] * one;
1630         cln::cl_N resbuf;
1631         int N = 0;
1632         do {
1633                 resbuf = res;
1634                 factor = factor * lambda;
1635                 N++;
1636                 res = res + crX[N] * factor / (N+Sqk);
1637         } while ((res != resbuf) || cln::zerop(crX[N]));
1638         return res;
1639 }
1640
1641
1642 // [Cra] section 4
1643 void calc_f(int maxr)
1644 {
1645         f_kj.clear();
1646         f_kj.resize(L1);
1647         
1648         cln::cl_N t0, t1, t2, t3, t4;
1649         int i, j, k;
1650         std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
1651         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1652         
1653         t0 = cln::exp(-lambda);
1654         t2 = 1;
1655         for (k=1; k<=L1; k++) {
1656                 t1 = k * lambda;
1657                 t2 = t0 * t2;
1658                 for (j=1; j<=maxr; j++) {
1659                         t3 = 1;
1660                         t4 = 1;
1661                         for (i=2; i<=j; i++) {
1662                                 t4 = t4 * (j-i+1);
1663                                 t3 = t1 * t3 + t4;
1664                         }
1665                         (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one);
1666                 }
1667                 it++;
1668         }
1669 }
1670
1671
1672 // [Cra] (3.1)
1673 cln::cl_N crandall_Z(const std::vector<int>& s)
1674 {
1675         const int j = s.size();
1676
1677         if (j == 1) {   
1678                 cln::cl_N t0;
1679                 cln::cl_N t0buf;
1680                 int q = 0;
1681                 do {
1682                         t0buf = t0;
1683                         q++;
1684                         t0 = t0 + f_kj[q+j-2][s[0]-1];
1685                 } while (t0 != t0buf);
1686                 
1687                 return t0 / cln::factorial(s[0]-1);
1688         }
1689
1690         std::vector<cln::cl_N> t(j);
1691
1692         cln::cl_N t0buf;
1693         int q = 0;
1694         do {
1695                 t0buf = t[0];
1696                 q++;
1697                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1698                 for (int k=j-2; k>=1; k--) {
1699                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1700                 }
1701                 t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
1702         } while (t[0] != t0buf);
1703         
1704         return t[0] / cln::factorial(s[0]-1);
1705 }
1706
1707
1708 // [Cra] (2.4)
1709 cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
1710 {
1711         std::vector<int> r = s;
1712         const int j = r.size();
1713
1714         // decide on maximal size of f_kj for crandall_Z
1715         if (Digits < 50) {
1716                 L1 = 150;
1717         } else {
1718                 L1 = Digits * 3 + j*2;
1719         }
1720
1721         // decide on maximal size of crX for crandall_Y
1722         if (Digits < 38) {
1723                 L2 = 63;
1724         } else if (Digits < 86) {
1725                 L2 = 127;
1726         } else if (Digits < 192) {
1727                 L2 = 255;
1728         } else if (Digits < 394) {
1729                 L2 = 511;
1730         } else if (Digits < 808) {
1731                 L2 = 1023;
1732         } else {
1733                 L2 = 2047;
1734         }
1735
1736         cln::cl_N res;
1737
1738         int maxr = 0;
1739         int S = 0;
1740         for (int i=0; i<j; i++) {
1741                 S += r[i];
1742                 if (r[i] > maxr) {
1743                         maxr = r[i];
1744                 }
1745         }
1746
1747         calc_f(maxr);
1748
1749         const cln::cl_N r0factorial = cln::factorial(r[0]-1);
1750
1751         std::vector<int> rz;
1752         int skp1buf;
1753         int Srun = S;
1754         for (int k=r.size()-1; k>0; k--) {
1755
1756                 rz.insert(rz.begin(), r.back());
1757                 skp1buf = rz.front();
1758                 Srun -= skp1buf;
1759                 r.pop_back();
1760
1761                 initcX(r);
1762                 
1763                 for (int q=0; q<skp1buf; q++) {
1764                         
1765                         cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
1766                         cln::cl_N pp2 = crandall_Z(rz);
1767
1768                         rz.front()--;
1769                         
1770                         if (q & 1) {
1771                                 res = res - pp1 * pp2 / cln::factorial(q);
1772                         } else {
1773                                 res = res + pp1 * pp2 / cln::factorial(q);
1774                         }
1775                 }
1776                 rz.front() = skp1buf;
1777         }
1778         rz.insert(rz.begin(), r.back());
1779
1780         initcX(rz);
1781
1782         res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
1783
1784         return res;
1785 }
1786
1787
1788 cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
1789 {
1790         const int j = r.size();
1791
1792         // buffer for subsums
1793         std::vector<cln::cl_N> t(j);
1794         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1795
1796         cln::cl_N t0buf;
1797         int q = 0;
1798         do {
1799                 t0buf = t[0];
1800                 q++;
1801                 t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
1802                 for (int k=j-2; k>=0; k--) {
1803                         t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
1804                 }
1805         } while (t[0] != t0buf);
1806
1807         return t[0];
1808 }
1809
1810
1811 // does Hoelder convolution. see [BBB] (7.0)
1812 cln::cl_N zeta_do_Hoelder_convolution(const std::vector<int>& m_, const std::vector<int>& s_)
1813 {
1814         // prepare parameters
1815         // holds Li arguments in [BBB] notation
1816         std::vector<int> s = s_;
1817         std::vector<int> m_p = m_;
1818         std::vector<int> m_q;
1819         // holds Li arguments in nested sums notation
1820         std::vector<cln::cl_N> s_p(s.size(), cln::cl_N(1));
1821         s_p[0] = s_p[0] * cln::cl_N("1/2");
1822         // convert notations
1823         int sig = 1;
1824         for (int i=0; i<s_.size(); i++) {
1825                 if (s_[i] < 0) {
1826                         sig = -sig;
1827                         s_p[i] = -s_p[i];
1828                 }
1829                 s[i] = sig * std::abs(s[i]);
1830         }
1831         std::vector<cln::cl_N> s_q;
1832         cln::cl_N signum = 1;
1833
1834         // first term
1835         cln::cl_N res = multipleLi_do_sum(m_p, s_p);
1836
1837         // middle terms
1838         do {
1839
1840                 // change parameters
1841                 if (s.front() > 0) {
1842                         if (m_p.front() == 1) {
1843                                 m_p.erase(m_p.begin());
1844                                 s_p.erase(s_p.begin());
1845                                 if (s_p.size() > 0) {
1846                                         s_p.front() = s_p.front() * cln::cl_N("1/2");
1847                                 }
1848                                 s.erase(s.begin());
1849                                 m_q.front()++;
1850                         } else {
1851                                 m_p.front()--;
1852                                 m_q.insert(m_q.begin(), 1);
1853                                 if (s_q.size() > 0) {
1854                                         s_q.front() = s_q.front() * 2;
1855                                 }
1856                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
1857                         }
1858                 } else {
1859                         if (m_p.front() == 1) {
1860                                 m_p.erase(m_p.begin());
1861                                 cln::cl_N spbuf = s_p.front();
1862                                 s_p.erase(s_p.begin());
1863                                 if (s_p.size() > 0) {
1864                                         s_p.front() = s_p.front() * spbuf;
1865                                 }
1866                                 s.erase(s.begin());
1867                                 m_q.insert(m_q.begin(), 1);
1868                                 if (s_q.size() > 0) {
1869                                         s_q.front() = s_q.front() * 4;
1870                                 }
1871                                 s_q.insert(s_q.begin(), cln::cl_N("1/4"));
1872                                 signum = -signum;
1873                         } else {
1874                                 m_p.front()--;
1875                                 m_q.insert(m_q.begin(), 1);
1876                                 if (s_q.size() > 0) {
1877                                         s_q.front() = s_q.front() * 2;
1878                                 }
1879                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
1880                         }
1881                 }
1882
1883                 // exiting the loop
1884                 if (m_p.size() == 0) break;
1885
1886                 res = res + signum * multipleLi_do_sum(m_p, s_p) * multipleLi_do_sum(m_q, s_q);
1887                 
1888         } while (true);
1889
1890         // last term
1891         res = res + signum * multipleLi_do_sum(m_q, s_q);
1892         
1893         return res;
1894 }
1895
1896
1897 } // end of anonymous namespace
1898
1899
1900 //////////////////////////////////////////////////////////////////////
1901 //
1902 // Multiple zeta values  zeta(x)
1903 //
1904 // GiNaC function
1905 //
1906 //////////////////////////////////////////////////////////////////////
1907
1908
1909 static ex zeta1_evalf(const ex& x)
1910 {
1911         if (is_exactly_a<lst>(x) && (x.nops()>1)) {
1912
1913                 // multiple zeta value
1914                 const int count = x.nops();
1915                 const lst& xlst = ex_to<lst>(x);
1916                 std::vector<int> r(count);
1917
1918                 // check parameters and convert them
1919                 lst::const_iterator it1 = xlst.begin();
1920                 std::vector<int>::iterator it2 = r.begin();
1921                 do {
1922                         if (!(*it1).info(info_flags::posint)) {
1923                                 return zeta(x).hold();
1924                         }
1925                         *it2 = ex_to<numeric>(*it1).to_int();
1926                         it1++;
1927                         it2++;
1928                 } while (it2 != r.end());
1929
1930                 // check for divergence
1931                 if (r[0] == 1) {
1932                         return zeta(x).hold();
1933                 }
1934
1935                 // decide on summation algorithm
1936                 // this is still a bit clumsy
1937                 int limit = (Digits>17) ? 10 : 6;
1938                 if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) {
1939                         return numeric(zeta_do_sum_Crandall(r));
1940                 } else {
1941                         return numeric(zeta_do_sum_simple(r));
1942                 }
1943         }
1944                 
1945         // single zeta value
1946         if (is_exactly_a<numeric>(x) && (x != 1)) {
1947                 try {
1948                         return zeta(ex_to<numeric>(x));
1949                 } catch (const dunno &e) { }
1950         }
1951
1952         return zeta(x).hold();
1953 }
1954
1955
1956 static ex zeta1_eval(const ex& x)
1957 {
1958         if (is_exactly_a<lst>(x)) {
1959                 if (x.nops() == 1) {
1960                         return zeta(x.op(0));
1961                 }
1962                 return zeta(x).hold();
1963         }
1964
1965         if (x.info(info_flags::numeric)) {
1966                 const numeric& y = ex_to<numeric>(x);
1967                 // trap integer arguments:
1968                 if (y.is_integer()) {
1969                         if (y.is_zero()) {
1970                                 return _ex_1_2;
1971                         }
1972                         if (y.is_equal(_num1)) {
1973                                 return zeta(x).hold();
1974                         }
1975                         if (y.info(info_flags::posint)) {
1976                                 if (y.info(info_flags::odd)) {
1977                                         return zeta(x).hold();
1978                                 } else {
1979                                         return abs(bernoulli(y)) * pow(Pi, y) * pow(_num2, y-_num1) / factorial(y);
1980                                 }
1981                         } else {
1982                                 if (y.info(info_flags::odd)) {
1983                                         return -bernoulli(_num1-y) / (_num1-y);
1984                                 } else {
1985                                         return _ex0;
1986                                 }
1987                         }
1988                 }
1989                 // zeta(float)
1990                 if (y.info(info_flags::numeric) && !y.info(info_flags::crational))
1991                         return zeta1_evalf(x);
1992         }
1993         return zeta(x).hold();
1994 }
1995
1996
1997 static ex zeta1_deriv(const ex& x, unsigned deriv_param)
1998 {
1999         GINAC_ASSERT(deriv_param==0);
2000
2001         if (is_exactly_a<lst>(x)) {
2002                 return _ex0;
2003         } else {
2004                 return zeta(_ex1, x);
2005         }
2006 }
2007
2008
2009 unsigned zeta1_SERIAL::serial =
2010                         function::register_new(function_options("zeta").
2011                                                 eval_func(zeta1_eval).
2012                                                 evalf_func(zeta1_evalf).
2013                                                 do_not_evalf_params().
2014                                                 derivative_func(zeta1_deriv).
2015                                                 latex_name("\\zeta").
2016                                                 overloaded(2));
2017
2018
2019 //////////////////////////////////////////////////////////////////////
2020 //
2021 // Alternating Euler sum  zeta(x,s)
2022 //
2023 // GiNaC function
2024 //
2025 //////////////////////////////////////////////////////////////////////
2026
2027
2028 static ex zeta2_evalf(const ex& x, const ex& s)
2029 {
2030         if (is_exactly_a<lst>(x)) {
2031
2032                 // alternating Euler sum
2033                 const int count = x.nops();
2034                 const lst& xlst = ex_to<lst>(x);
2035                 const lst& slst = ex_to<lst>(s);
2036                 std::vector<int> xi(count);
2037                 std::vector<int> si(count);
2038
2039                 // check parameters and convert them
2040                 lst::const_iterator it_xread = xlst.begin();
2041                 lst::const_iterator it_sread = slst.begin();
2042                 std::vector<int>::iterator it_xwrite = xi.begin();
2043                 std::vector<int>::iterator it_swrite = si.begin();
2044                 do {
2045                         if (!(*it_xread).info(info_flags::posint)) {
2046                                 return zeta(x, s).hold();
2047                         }
2048                         *it_xwrite = ex_to<numeric>(*it_xread).to_int();
2049                         if (*it_sread > 0) {
2050                                 *it_swrite = 1;
2051                         } else {
2052                                 *it_swrite = -1;
2053                         }
2054                         it_xread++;
2055                         it_sread++;
2056                         it_xwrite++;
2057                         it_swrite++;
2058                 } while (it_xwrite != xi.end());
2059
2060                 // check for divergence
2061                 if ((xi[0] == 1) && (si[0] == 1)) {
2062                         return zeta(x, s).hold();
2063                 }
2064
2065                 // use Hoelder convolution
2066                 return numeric(zeta_do_Hoelder_convolution(xi, si));
2067         }
2068                 
2069         return zeta(x, s).hold();
2070 }
2071
2072
2073 static ex zeta2_eval(const ex& x, const ex& s)
2074 {
2075         if (is_exactly_a<lst>(s)) {
2076                 const lst& l = ex_to<lst>(s);
2077                 lst::const_iterator it = l.begin();
2078                 while (it != l.end()) {
2079                         if ((*it).info(info_flags::negative)) {
2080                                 return zeta(x, s).hold();
2081                         }
2082                         it++;
2083                 }
2084                 return zeta(x);
2085         } else {
2086                 if (s.info(info_flags::positive)) {
2087                         return zeta(x);
2088                 }
2089         }
2090
2091         return zeta(x, s).hold();
2092 }
2093
2094
2095 static ex zeta2_deriv(const ex& x, const ex& s, unsigned deriv_param)
2096 {
2097         GINAC_ASSERT(deriv_param==0);
2098
2099         if (is_exactly_a<lst>(x)) {
2100                 return _ex0;
2101         } else {
2102                 if ((is_exactly_a<lst>(s) && (s.op(0) > 0)) || (s > 0)) {
2103                         return zeta(_ex1, x);
2104                 }
2105                 return _ex0;
2106         }
2107 }
2108
2109
2110 unsigned zeta2_SERIAL::serial =
2111                         function::register_new(function_options("zeta").
2112                                                 eval_func(zeta2_eval).
2113                                                 evalf_func(zeta2_evalf).
2114                                                 do_not_evalf_params().
2115                                                 derivative_func(zeta2_deriv).
2116                                                 latex_name("\\zeta").
2117                                                 overloaded(2));
2118
2119
2120 } // namespace GiNaC
2121