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