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