inifcns_nstdsums.cpp: Lin_numeric takes cl_N as an argument instead of numeric.
[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  *                                         G(lst(a_1,...,a_k),y) or G(lst(a_1,...,a_k),lst(s_1,...,s_k),y)
9  *    Nielsen's generalized polylogarithm  S(n,p,x)
10  *    harmonic polylogarithm               H(m,x) or H(lst(m_1,...,m_k),x)
11  *    multiple zeta value                  zeta(m) or zeta(lst(m_1,...,m_k))
12  *    alternating Euler sum                zeta(m,s) or zeta(lst(m_1,...,m_k),lst(s_1,...,s_k))
13  *
14  *  Some remarks:
15  *
16  *    - All formulae used can be looked up in the following publications:
17  *      [Kol] Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
18  *      [Cra] Fast Evaluation of Multiple Zeta Sums, R.E.Crandall, Math.Comp. 67 (1998), pp. 1163-1172.
19  *      [ReV] Harmonic Polylogarithms, E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
20  *      [BBB] Special Values of Multiple Polylogarithms, J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001), pp. 907-941
21  *      [VSW] Numerical evaluation of multiple polylogarithms, J.Vollinga, S.Weinzierl, hep-ph/0410259
22  *
23  *    - The order of parameters and arguments of Li and zeta is defined according to the nested sums
24  *      representation. The parameters for H are understood as in [ReV]. They can be in expanded --- only
25  *      0, 1 and -1 --- or in compactified --- a string with zeros in front of 1 or -1 is written as a single
26  *      number --- notation.
27  *
28  *    - All functions can be nummerically evaluated with arguments in the whole complex plane. The parameters
29  *      for Li, zeta and S must be positive integers. If you want to have an alternating Euler sum, you have
30  *      to give the signs of the parameters as a second argument s to zeta(m,s) containing 1 and -1.
31  *
32  *    - The calculation of classical polylogarithms is speeded up by using Bernoulli numbers and 
33  *      look-up tables. S uses look-up tables as well. The zeta function applies the algorithms in
34  *      [Cra] and [BBB] for speed up. Multiple polylogarithms use Hoelder convolution [BBB].
35  *
36  *    - The functions have no means to do a series expansion into nested sums. To do this, you have to convert
37  *      these functions into the appropriate objects from the nestedsums library, do the expansion and convert
38  *      the result back.
39  *
40  *    - Numerical testing of this implementation has been performed by doing a comparison of results
41  *      between this software and the commercial M.......... 4.1. Multiple zeta values have been checked
42  *      by means of evaluations into simple zeta values. Harmonic polylogarithms have been checked by
43  *      comparison to S(n,p,x) for corresponding parameter combinations and by continuity checks
44  *      around |x|=1 along with comparisons to corresponding zeta functions. Multiple polylogarithms were
45  *      checked against H and zeta and by means of shuffle and quasi-shuffle relations.
46  *
47  */
48
49 /*
50  *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
51  *
52  *  This program is free software; you can redistribute it and/or modify
53  *  it under the terms of the GNU General Public License as published by
54  *  the Free Software Foundation; either version 2 of the License, or
55  *  (at your option) any later version.
56  *
57  *  This program is distributed in the hope that it will be useful,
58  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
59  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
60  *  GNU General Public License for more details.
61  *
62  *  You should have received a copy of the GNU General Public License
63  *  along with this program; if not, write to the Free Software
64  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
65  */
66
67 #include <sstream>
68 #include <stdexcept>
69 #include <vector>
70 #include <cln/cln.h>
71
72 #include "inifcns.h"
73
74 #include "add.h"
75 #include "constant.h"
76 #include "lst.h"
77 #include "mul.h"
78 #include "numeric.h"
79 #include "operators.h"
80 #include "power.h"
81 #include "pseries.h"
82 #include "relational.h"
83 #include "symbol.h"
84 #include "utils.h"
85 #include "wildcard.h"
86
87
88 namespace GiNaC {
89
90
91 //////////////////////////////////////////////////////////////////////
92 //
93 // Classical polylogarithm  Li(n,x)
94 //
95 // helper functions
96 //
97 //////////////////////////////////////////////////////////////////////
98
99
100 // anonymous namespace for helper functions
101 namespace {
102
103
104 // lookup table for factors built from Bernoulli numbers
105 // see fill_Xn()
106 std::vector<std::vector<cln::cl_N> > Xn;
107 // initial size of Xn that should suffice for 32bit machines (must be even)
108 const int xninitsizestep = 26;
109 int xninitsize = xninitsizestep;
110 int xnsize = 0;
111
112
113 // This function calculates the X_n. The X_n are needed for speed up of classical polylogarithms.
114 // With these numbers the polylogs can be calculated as follows:
115 //   Li_p (x)  =  \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with  u = -log(1-x)
116 //   X_0(n) = B_n (Bernoulli numbers)
117 //   X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
118 // The calculation of Xn depends on X0 and X{n-1}.
119 // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
120 // This results in a slightly more complicated algorithm for the X_n.
121 // The first index in Xn corresponds to the index of the polylog minus 2.
122 // The second index in Xn corresponds to the index from the actual sum.
123 void fill_Xn(int n)
124 {
125         if (n>1) {
126                 // calculate X_2 and higher (corresponding to Li_4 and higher)
127                 std::vector<cln::cl_N> buf(xninitsize);
128                 std::vector<cln::cl_N>::iterator it = buf.begin();
129                 cln::cl_N result;
130                 *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
131                 it++;
132                 for (int i=2; i<=xninitsize; i++) {
133                         if (i&1) {
134                                 result = 0; // k == 0
135                         } else {
136                                 result = Xn[0][i/2-1]; // k == 0
137                         }
138                         for (int k=1; k<i-1; k++) {
139                                 if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
140                                         result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
141                                 }
142                         }
143                         result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
144                         result = result + Xn[n-1][i-1] / (i+1); // k == i
145                         
146                         *it = result;
147                         it++;
148                 }
149                 Xn.push_back(buf);
150         } else if (n==1) {
151                 // special case to handle the X_0 correct
152                 std::vector<cln::cl_N> buf(xninitsize);
153                 std::vector<cln::cl_N>::iterator it = buf.begin();
154                 cln::cl_N result;
155                 *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
156                 it++;
157                 *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
158                 it++;
159                 for (int i=3; i<=xninitsize; i++) {
160                         if (i & 1) {
161                                 result = -Xn[0][(i-3)/2]/2;
162                                 *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
163                                 it++;
164                         } else {
165                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
166                                 for (int k=1; k<i/2; k++) {
167                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
168                                 }
169                                 *it = result;
170                                 it++;
171                         }
172                 }
173                 Xn.push_back(buf);
174         } else {
175                 // calculate X_0
176                 std::vector<cln::cl_N> buf(xninitsize/2);
177                 std::vector<cln::cl_N>::iterator it = buf.begin();
178                 for (int i=1; i<=xninitsize/2; i++) {
179                         *it = bernoulli(i*2).to_cl_N();
180                         it++;
181                 }
182                 Xn.push_back(buf);
183         }
184
185         xnsize++;
186 }
187
188 // doubles the number of entries in each Xn[]
189 void double_Xn()
190 {
191         const int pos0 = xninitsize / 2;
192         // X_0
193         for (int i=1; i<=xninitsizestep/2; ++i) {
194                 Xn[0].push_back(bernoulli((i+pos0)*2).to_cl_N());
195         }
196         if (Xn.size() > 1) {
197                 int xend = xninitsize + xninitsizestep;
198                 cln::cl_N result;
199                 // X_1
200                 for (int i=xninitsize+1; i<=xend; ++i) {
201                         if (i & 1) {
202                                 result = -Xn[0][(i-3)/2]/2;
203                                 Xn[1].push_back((cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result);
204                         } else {
205                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
206                                 for (int k=1; k<i/2; k++) {
207                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
208                                 }
209                                 Xn[1].push_back(result);
210                         }
211                 }
212                 // X_n
213                 for (int n=2; n<Xn.size(); ++n) {
214                         for (int i=xninitsize+1; i<=xend; ++i) {
215                                 if (i & 1) {
216                                         result = 0; // k == 0
217                                 } else {
218                                         result = Xn[0][i/2-1]; // k == 0
219                                 }
220                                 for (int k=1; k<i-1; ++k) {
221                                         if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
222                                                 result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
223                                         }
224                                 }
225                                 result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
226                                 result = result + Xn[n-1][i-1] / (i+1); // k == i
227                                 Xn[n].push_back(result);
228                         }
229                 }
230         }
231         xninitsize += xninitsizestep;
232 }
233
234
235 // calculates Li(2,x) without Xn
236 cln::cl_N Li2_do_sum(const cln::cl_N& x)
237 {
238         cln::cl_N res = x;
239         cln::cl_N resbuf;
240         cln::cl_N num = x * cln::cl_float(1, cln::float_format(Digits));
241         cln::cl_I den = 1; // n^2 = 1
242         unsigned i = 3;
243         do {
244                 resbuf = res;
245                 num = num * x;
246                 den = den + i;  // n^2 = 4, 9, 16, ...
247                 i += 2;
248                 res = res + num / den;
249         } while (res != resbuf);
250         return res;
251 }
252
253
254 // calculates Li(2,x) with Xn
255 cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x)
256 {
257         std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
258         std::vector<cln::cl_N>::const_iterator xend = Xn[0].end();
259         cln::cl_N u = -cln::log(1-x);
260         cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
261         cln::cl_N uu = cln::square(u);
262         cln::cl_N res = u - uu/4;
263         cln::cl_N resbuf;
264         unsigned i = 1;
265         do {
266                 resbuf = res;
267                 factor = factor * uu / (2*i * (2*i+1));
268                 res = res + (*it) * factor;
269                 i++;
270                 if (++it == xend) {
271                         double_Xn();
272                         it = Xn[0].begin() + (i-1);
273                         xend = Xn[0].end();
274                 }
275         } while (res != resbuf);
276         return res;
277 }
278
279
280 // calculates Li(n,x), n>2 without Xn
281 cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
282 {
283         cln::cl_N factor = x * cln::cl_float(1, cln::float_format(Digits));
284         cln::cl_N res = x;
285         cln::cl_N resbuf;
286         int i=2;
287         do {
288                 resbuf = res;
289                 factor = factor * x;
290                 res = res + factor / cln::expt(cln::cl_I(i),n);
291                 i++;
292         } while (res != resbuf);
293         return res;
294 }
295
296
297 // calculates Li(n,x), n>2 with Xn
298 cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
299 {
300         std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
301         std::vector<cln::cl_N>::const_iterator xend = Xn[n-2].end();
302         cln::cl_N u = -cln::log(1-x);
303         cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
304         cln::cl_N res = u;
305         cln::cl_N resbuf;
306         unsigned i=2;
307         do {
308                 resbuf = res;
309                 factor = factor * u / i;
310                 res = res + (*it) * factor;
311                 i++;
312                 if (++it == xend) {
313                         double_Xn();
314                         it = Xn[n-2].begin() + (i-2);
315                         xend = Xn[n-2].end();
316                 }
317         } while (res != resbuf);
318         return res;
319 }
320
321
322 // forward declaration needed by function Li_projection and C below
323 const cln::cl_N S_num(int n, int p, const cln::cl_N& x);
324
325
326 // helper function for classical polylog Li
327 cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
328 {
329         // treat n=2 as special case
330         if (n == 2) {
331                 // check if precalculated X0 exists
332                 if (xnsize == 0) {
333                         fill_Xn(0);
334                 }
335
336                 if (cln::realpart(x) < 0.5) {
337                         // choose the faster algorithm
338                         // the switching point was empirically determined. the optimal point
339                         // depends on hardware, Digits, ... so an approx value is okay.
340                         // it solves also the problem with precision due to the u=-log(1-x) transformation
341                         if (cln::abs(cln::realpart(x)) < 0.25) {
342                                 
343                                 return Li2_do_sum(x);
344                         } else {
345                                 return Li2_do_sum_Xn(x);
346                         }
347                 } else {
348                         // choose the faster algorithm
349                         if (cln::abs(cln::realpart(x)) > 0.75) {
350                                 return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
351                         } else {
352                                 return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
353                         }
354                 }
355         } else {
356                 // check if precalculated Xn exist
357                 if (n > xnsize+1) {
358                         for (int i=xnsize; i<n-1; i++) {
359                                 fill_Xn(i);
360                         }
361                 }
362
363                 if (cln::realpart(x) < 0.5) {
364                         // choose the faster algorithm
365                         // with n>=12 the "normal" summation always wins against the method with Xn
366                         if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
367                                 return Lin_do_sum(n, x);
368                         } else {
369                                 return Lin_do_sum_Xn(n, x);
370                         }
371                 } else {
372                         cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
373                         for (int j=0; j<n-1; j++) {
374                                 result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
375                                                   * cln::expt(cln::log(x), j) / cln::factorial(j);
376                         }
377                         return result;
378                 }
379         }
380 }
381
382 // helper function for classical polylog Li
383 const cln::cl_N Lin_numeric(const int n, const cln::cl_N& x)
384 {
385         if (n == 1) {
386                 // just a log
387                 return -cln::log(1-x);
388         }
389         if (zerop(x)) {
390                 return 0;
391         }
392         if (x == 1) {
393                 // [Kol] (2.22)
394                 return cln::zeta(n);
395         }
396         else if (x == -1) {
397                 // [Kol] (2.22)
398                 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
399         }
400         if (cln::abs(realpart(x)) < 0.4 && cln::abs(cln::abs(x)-1) < 0.01) {
401                 cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
402                 for (int j=0; j<n-1; j++) {
403                         result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
404                                 * cln::expt(cln::log(x), j) / cln::factorial(j);
405                 }
406                 return result;
407         }
408
409         // what is the desired float format?
410         // first guess: default format
411         cln::float_format_t prec = cln::default_float_format;
412         const cln::cl_N value = x;
413         // second guess: the argument's format
414         if (!instanceof(realpart(x), cln::cl_RA_ring))
415                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
416         else if (!instanceof(imagpart(x), cln::cl_RA_ring))
417                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
418         
419         // [Kol] (5.15)
420         if (cln::abs(value) > 1) {
421                 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
422                 // check if argument is complex. if it is real, the new polylog has to be conjugated.
423                 if (cln::zerop(cln::imagpart(value))) {
424                         if (n & 1) {
425                                 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
426                         }
427                         else {
428                                 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
429                         }
430                 }
431                 else {
432                         if (n & 1) {
433                                 result = result + Li_projection(n, cln::recip(value), prec);
434                         }
435                         else {
436                                 result = result - Li_projection(n, cln::recip(value), prec);
437                         }
438                 }
439                 cln::cl_N add;
440                 for (int j=0; j<n-1; j++) {
441                         add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
442                                     * Lin_numeric(n-j,1) * cln::expt(cln::log(-value),j) / cln::factorial(j);
443                 }
444                 result = result - add;
445                 return result;
446         }
447         else {
448                 return Li_projection(n, value, prec);
449         }
450 }
451
452
453 } // end of anonymous namespace
454
455
456 //////////////////////////////////////////////////////////////////////
457 //
458 // Multiple polylogarithm  Li(n,x)
459 //
460 // helper function
461 //
462 //////////////////////////////////////////////////////////////////////
463
464
465 // anonymous namespace for helper function
466 namespace {
467
468
469 // performs the actual series summation for multiple polylogarithms
470 cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
471 {
472         // ensure all x <> 0.
473         for (std::vector<cln::cl_N>::const_iterator it = x.begin(); it != x.end(); ++it) {
474                 if ( *it == 0 ) return cln::cl_float(0, cln::float_format(Digits));
475         }
476
477         const int j = s.size();
478         bool flag_accidental_zero = false;
479
480         std::vector<cln::cl_N> t(j);
481         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
482
483         cln::cl_N t0buf;
484         int q = 0;
485         do {
486                 t0buf = t[0];
487                 q++;
488                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
489                 for (int k=j-2; k>=0; k--) {
490                         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]);
491                 }
492                 q++;
493                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
494                 for (int k=j-2; k>=0; k--) {
495                         flag_accidental_zero = cln::zerop(t[k+1]);
496                         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]);
497                 }
498         } while ( (t[0] != t0buf) || cln::zerop(t[0]) || flag_accidental_zero );
499
500         return t[0];
501 }
502
503
504 // converts parameter types and calls multipleLi_do_sum (convenience function for G_numeric)
505 cln::cl_N mLi_do_summation(const lst& m, const lst& x)
506 {
507         std::vector<int> m_int;
508         std::vector<cln::cl_N> x_cln;
509         for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
510                 m_int.push_back(ex_to<numeric>(*itm).to_int());
511                 x_cln.push_back(ex_to<numeric>(*itx).to_cl_N());
512         }
513         return multipleLi_do_sum(m_int, x_cln);
514 }
515
516
517 // forward declaration for Li_eval()
518 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf);
519
520
521 // holding dummy-symbols for the G/Li transformations
522 std::vector<ex> gsyms;
523
524
525 // type used by the transformation functions for G
526 typedef std::vector<int> Gparameter;
527
528
529 // G_eval1-function for G transformations
530 ex G_eval1(int a, int scale)
531 {
532         if (a != 0) {
533                 const ex& scs = gsyms[std::abs(scale)];
534                 const ex& as = gsyms[std::abs(a)];
535                 if (as != scs) {
536                         return -log(1 - scs/as);
537                 } else {
538                         return -zeta(1);
539                 }
540         } else {
541                 return log(gsyms[std::abs(scale)]);
542         }
543 }
544
545
546 // G_eval-function for G transformations
547 ex G_eval(const Gparameter& a, int scale)
548 {
549         // check for properties of G
550         ex sc = gsyms[std::abs(scale)];
551         lst newa;
552         bool all_zero = true;
553         bool all_ones = true;
554         int count_ones = 0;
555         for (Gparameter::const_iterator it = a.begin(); it != a.end(); ++it) {
556                 if (*it != 0) {
557                         const ex sym = gsyms[std::abs(*it)];
558                         newa.append(sym);
559                         all_zero = false;
560                         if (sym != sc) {
561                                 all_ones = false;
562                         }
563                         if (all_ones) {
564                                 ++count_ones;
565                         }
566                 } else {
567                         all_ones = false;
568                 }
569         }
570
571         // care about divergent G: shuffle to separate divergencies that will be canceled
572         // later on in the transformation
573         if (newa.nops() > 1 && newa.op(0) == sc && !all_ones && a.front()!=0) {
574                 // do shuffle
575                 Gparameter short_a;
576                 Gparameter::const_iterator it = a.begin();
577                 ++it;
578                 for (; it != a.end(); ++it) {
579                         short_a.push_back(*it);
580                 }
581                 ex result = G_eval1(a.front(), scale) * G_eval(short_a, scale);
582                 it = short_a.begin();
583                 for (int i=1; i<count_ones; ++i) {
584                         ++it;
585                 }
586                 for (; it != short_a.end(); ++it) {
587
588                         Gparameter newa;
589                         Gparameter::const_iterator it2 = short_a.begin();
590                         for (--it2; it2 != it;) {
591                                 ++it2;
592                                 newa.push_back(*it2);
593                         }
594                         newa.push_back(a[0]);
595                         ++it2;
596                         for (; it2 != short_a.end(); ++it2) {
597                                 newa.push_back(*it2);   
598                         }
599                         result -= G_eval(newa, scale);
600                 }
601                 return result / count_ones;
602         }
603
604         // G({1,...,1};y) -> G({1};y)^k / k!
605         if (all_ones && a.size() > 1) {
606                 return pow(G_eval1(a.front(),scale), count_ones) / factorial(count_ones);
607         }
608
609         // G({0,...,0};y) -> log(y)^k / k!
610         if (all_zero) {
611                 return pow(log(gsyms[std::abs(scale)]), a.size()) / factorial(a.size());
612         }
613
614         // no special cases anymore -> convert it into Li
615         lst m;
616         lst x;
617         ex argbuf = gsyms[std::abs(scale)];
618         ex mval = _ex1;
619         for (Gparameter::const_iterator it=a.begin(); it!=a.end(); ++it) {
620                 if (*it != 0) {
621                         const ex& sym = gsyms[std::abs(*it)];
622                         x.append(argbuf / sym);
623                         m.append(mval);
624                         mval = _ex1;
625                         argbuf = sym;
626                 } else {
627                         ++mval;
628                 }
629         }
630         return pow(-1, x.nops()) * Li(m, x);
631 }
632
633
634 // converts data for G: pending_integrals -> a
635 Gparameter convert_pending_integrals_G(const Gparameter& pending_integrals)
636 {
637         GINAC_ASSERT(pending_integrals.size() != 1);
638
639         if (pending_integrals.size() > 0) {
640                 // get rid of the first element, which would stand for the new upper limit
641                 Gparameter new_a(pending_integrals.begin()+1, pending_integrals.end());
642                 return new_a;
643         } else {
644                 // just return empty parameter list
645                 Gparameter new_a;
646                 return new_a;
647         }
648 }
649
650
651 // check the parameters a and scale for G and return information about convergence, depth, etc.
652 // convergent     : true if G(a,scale) is convergent
653 // depth          : depth of G(a,scale)
654 // trailing_zeros : number of trailing zeros of a
655 // min_it         : iterator of a pointing on the smallest element in a
656 Gparameter::const_iterator check_parameter_G(const Gparameter& a, int scale,
657                 bool& convergent, int& depth, int& trailing_zeros, Gparameter::const_iterator& min_it)
658 {
659         convergent = true;
660         depth = 0;
661         trailing_zeros = 0;
662         min_it = a.end();
663         Gparameter::const_iterator lastnonzero = a.end();
664         for (Gparameter::const_iterator it = a.begin(); it != a.end(); ++it) {
665                 if (std::abs(*it) > 0) {
666                         ++depth;
667                         trailing_zeros = 0;
668                         lastnonzero = it;
669                         if (std::abs(*it) < scale) {
670                                 convergent = false;
671                                 if ((min_it == a.end()) || (std::abs(*it) < std::abs(*min_it))) {
672                                         min_it = it;
673                                 }
674                         }
675                 } else {
676                         ++trailing_zeros;
677                 }
678         }
679         return ++lastnonzero;
680 }
681
682
683 // add scale to pending_integrals if pending_integrals is empty
684 Gparameter prepare_pending_integrals(const Gparameter& pending_integrals, int scale)
685 {
686         GINAC_ASSERT(pending_integrals.size() != 1);
687
688         if (pending_integrals.size() > 0) {
689                 return pending_integrals;
690         } else {
691                 Gparameter new_pending_integrals;
692                 new_pending_integrals.push_back(scale);
693                 return new_pending_integrals;
694         }
695 }
696
697
698 // handles trailing zeroes for an otherwise convergent integral
699 ex trailing_zeros_G(const Gparameter& a, int scale)
700 {
701         bool convergent;
702         int depth, trailing_zeros;
703         Gparameter::const_iterator last, dummyit;
704         last = check_parameter_G(a, scale, convergent, depth, trailing_zeros, dummyit);
705
706         GINAC_ASSERT(convergent);
707
708         if ((trailing_zeros > 0) && (depth > 0)) {
709                 ex result;
710                 Gparameter new_a(a.begin(), a.end()-1);
711                 result += G_eval1(0, scale) * trailing_zeros_G(new_a, scale);
712                 for (Gparameter::const_iterator it = a.begin(); it != last; ++it) {
713                         Gparameter new_a(a.begin(), it);
714                         new_a.push_back(0);
715                         new_a.insert(new_a.end(), it, a.end()-1);
716                         result -= trailing_zeros_G(new_a, scale);
717                 }
718
719                 return result / trailing_zeros;
720         } else {
721                 return G_eval(a, scale);
722         }
723 }
724
725
726 // G transformation [VSW] (57),(58)
727 ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, int scale)
728 {
729         // pendint = ( y1, b1, ..., br )
730         //       a = ( 0, ..., 0, amin )
731         //   scale = y2
732         //
733         // int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(0, ..., 0, sr; y2)
734         // where sr replaces amin
735
736         GINAC_ASSERT(a.back() != 0);
737         GINAC_ASSERT(a.size() > 0);
738
739         ex result;
740         Gparameter new_pending_integrals = prepare_pending_integrals(pending_integrals, std::abs(a.back()));
741         const int psize = pending_integrals.size();
742
743         // length == 1
744         // G(sr_{+-}; y2 ) = G(y2_{-+}; sr) - G(0; sr) + ln(-y2_{-+})
745
746         if (a.size() == 1) {
747
748           // ln(-y2_{-+})
749           result += log(gsyms[ex_to<numeric>(scale).to_int()]);
750                 if (a.back() > 0) {
751                         new_pending_integrals.push_back(-scale);
752                         result += I*Pi;
753                 } else {
754                         new_pending_integrals.push_back(scale);
755                         result -= I*Pi;
756                 }
757                 if (psize) {
758                         result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front());
759                 }
760                 
761                 // G(y2_{-+}; sr)
762                 result += trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front());
763                 
764                 // G(0; sr)
765                 new_pending_integrals.back() = 0;
766                 result -= trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front());
767
768                 return result;
769         }
770
771         // length > 1
772         // G_m(sr_{+-}; y2) = -zeta_m + int_0^y2 dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
773         //                            - int_0^sr dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
774
775         //term zeta_m
776         result -= zeta(a.size());
777         if (psize) {
778                 result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front());
779         }
780         
781         // term int_0^sr dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
782         //    = int_0^sr dt/t G_{m-1}( t_{+-}; y2 )
783         Gparameter new_a(a.begin()+1, a.end());
784         new_pending_integrals.push_back(0);
785         result -= depth_one_trafo_G(new_pending_integrals, new_a, scale);
786         
787         // term int_0^y2 dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
788         //    = int_0^y2 dt/t G_{m-1}( t_{+-}; y2 )
789         Gparameter new_pending_integrals_2;
790         new_pending_integrals_2.push_back(scale);
791         new_pending_integrals_2.push_back(0);
792         if (psize) {
793                 result += trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front())
794                           * depth_one_trafo_G(new_pending_integrals_2, new_a, scale);
795         } else {
796                 result += depth_one_trafo_G(new_pending_integrals_2, new_a, scale);
797         }
798
799         return result;
800 }
801
802
803 // forward declaration
804 ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2,
805              const Gparameter& pendint, const Gparameter& a_old, int scale);
806
807
808 // G transformation [VSW]
809 ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
810 {
811         // main recursion routine
812         //
813         // pendint = ( y1, b1, ..., br )
814         //       a = ( a1, ..., amin, ..., aw )
815         //   scale = y2
816         //
817         // int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(a1,...,sr,...,aw,y2)
818         // where sr replaces amin
819
820         // find smallest alpha, determine depth and trailing zeros, and check for convergence
821         bool convergent;
822         int depth, trailing_zeros;
823         Gparameter::const_iterator min_it;
824         Gparameter::const_iterator firstzero = 
825                 check_parameter_G(a, scale, convergent, depth, trailing_zeros, min_it);
826         int min_it_pos = min_it - a.begin();
827
828         // special case: all a's are zero
829         if (depth == 0) {
830                 ex result;
831
832                 if (a.size() == 0) {
833                   result = 1;
834                 } else {
835                   result = G_eval(a, scale);
836                 }
837                 if (pendint.size() > 0) {
838                   result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front());
839                 } 
840                 return result;
841         }
842
843         // handle trailing zeros
844         if (trailing_zeros > 0) {
845                 ex result;
846                 Gparameter new_a(a.begin(), a.end()-1);
847                 result += G_eval1(0, scale) * G_transform(pendint, new_a, scale);
848                 for (Gparameter::const_iterator it = a.begin(); it != firstzero; ++it) {
849                         Gparameter new_a(a.begin(), it);
850                         new_a.push_back(0);
851                         new_a.insert(new_a.end(), it, a.end()-1);
852                         result -= G_transform(pendint, new_a, scale);
853                 }
854                 return result / trailing_zeros;
855         }
856
857         // convergence case
858         if (convergent) {
859                 if (pendint.size() > 0) {
860                         return G_eval(convert_pending_integrals_G(pendint), pendint.front()) * G_eval(a, scale);
861                 } else {
862                         return G_eval(a, scale);
863                 }
864         }
865
866         // call basic transformation for depth equal one
867         if (depth == 1) {
868                 return depth_one_trafo_G(pendint, a, scale);
869         }
870
871         // do recursion
872         // int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(a1,...,sr,...,aw,y2)
873         //  =  int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(a1,...,0,...,aw,y2)
874         //   + int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) int_0^{sr} ds_{r+1} d/ds_{r+1} G(a1,...,s_{r+1},...,aw,y2)
875
876         // smallest element in last place
877         if (min_it + 1 == a.end()) {
878                 do { --min_it; } while (*min_it == 0);
879                 Gparameter empty;
880                 Gparameter a1(a.begin(),min_it+1);
881                 Gparameter a2(min_it+1,a.end());
882
883                 ex result = G_transform(pendint,a2,scale)*G_transform(empty,a1,scale);
884
885                 result -= shuffle_G(empty,a1,a2,pendint,a,scale);
886                 return result;
887         }
888
889         Gparameter empty;
890         Gparameter::iterator changeit;
891
892         // first term G(a_1,..,0,...,a_w;a_0)
893         Gparameter new_pendint = prepare_pending_integrals(pendint, a[min_it_pos]);
894         Gparameter new_a = a;
895         new_a[min_it_pos] = 0;
896         ex result = G_transform(empty, new_a, scale);
897         if (pendint.size() > 0) {
898                 result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front());
899         }
900
901         // other terms
902         changeit = new_a.begin() + min_it_pos;
903         changeit = new_a.erase(changeit);
904         if (changeit != new_a.begin()) {
905                 // smallest in the middle
906                 new_pendint.push_back(*changeit);
907                 result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
908                         * G_transform(empty, new_a, scale);
909                 int buffer = *changeit;
910                 *changeit = *min_it;
911                 result += G_transform(new_pendint, new_a, scale);
912                 *changeit = buffer;
913                 new_pendint.pop_back();
914                 --changeit;
915                 new_pendint.push_back(*changeit);
916                 result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
917                         * G_transform(empty, new_a, scale);
918                 *changeit = *min_it;
919                 result -= G_transform(new_pendint, new_a, scale);
920         } else {
921                 // smallest at the front
922                 new_pendint.push_back(scale);
923                 result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
924                         * G_transform(empty, new_a, scale);
925                 new_pendint.back() =  *changeit;
926                 result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
927                         * G_transform(empty, new_a, scale);
928                 *changeit = *min_it;
929                 result += G_transform(new_pendint, new_a, scale);
930         }
931         return result;
932 }
933
934
935 // shuffles the two parameter list a1 and a2 and calls G_transform for every term except
936 // for the one that is equal to a_old
937 ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2,
938              const Gparameter& pendint, const Gparameter& a_old, int scale) 
939 {
940         if (a1.size()==0 && a2.size()==0) {
941                 // veto the one configuration we don't want
942                 if ( a0 == a_old ) return 0;
943
944                 return G_transform(pendint,a0,scale);
945         }
946
947         if (a2.size()==0) {
948                 Gparameter empty;
949                 Gparameter aa0 = a0;
950                 aa0.insert(aa0.end(),a1.begin(),a1.end());
951                 return shuffle_G(aa0,empty,empty,pendint,a_old,scale);
952         }
953
954         if (a1.size()==0) {
955                 Gparameter empty;
956                 Gparameter aa0 = a0;
957                 aa0.insert(aa0.end(),a2.begin(),a2.end());
958                 return shuffle_G(aa0,empty,empty,pendint,a_old,scale);
959         }
960
961         Gparameter a1_removed(a1.begin()+1,a1.end());
962         Gparameter a2_removed(a2.begin()+1,a2.end());
963
964         Gparameter a01 = a0;
965         Gparameter a02 = a0;
966
967         a01.push_back( a1[0] );
968         a02.push_back( a2[0] );
969
970         return shuffle_G(a01,a1_removed,a2,pendint,a_old,scale)
971              + shuffle_G(a02,a1,a2_removed,pendint,a_old,scale);
972 }
973
974
975 // handles the transformations and the numerical evaluation of G
976 // the parameter x, s and y must only contain numerics
977 ex G_numeric(const lst& x, const lst& s, const ex& y)
978 {
979         // check for convergence and necessary accelerations
980         bool need_trafo = false;
981         bool need_hoelder = false;
982         int depth = 0;
983         for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
984                 if (!(*it).is_zero()) {
985                         ++depth;
986                         if (abs(*it) - y < -pow(10,-Digits+1)) {
987                                 need_trafo = true;
988                         }
989                         if (abs((abs(*it) - y)/y) < 0.01) {
990                                 need_hoelder = true;
991                         }
992                 }
993         }
994         if (x.op(x.nops()-1).is_zero()) {
995                 need_trafo = true;
996         }
997         if (depth == 1 && x.nops() == 2 && !need_trafo) {
998                 return -Li(x.nops(), y / x.op(x.nops()-1)).evalf();
999         }
1000         
1001         // do acceleration transformation (hoelder convolution [BBB])
1002         if (need_hoelder) {
1003                 
1004                 ex result;
1005                 const int size = x.nops();
1006                 lst newx;
1007                 for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
1008                         newx.append(*it / y);
1009                 }
1010                 
1011                 for (int r=0; r<=size; ++r) {
1012                         ex buffer = pow(-1, r);
1013                         ex p = 2;
1014                         bool adjustp;
1015                         do {
1016                                 adjustp = false;
1017                                 for (lst::const_iterator it = newx.begin(); it != newx.end(); ++it) {
1018                                         if (*it == 1/p) {
1019                                                 p += (3-p)/2; 
1020                                                 adjustp = true;
1021                                                 continue;
1022                                         }
1023                                 }
1024                         } while (adjustp);
1025                         ex q = p / (p-1);
1026                         lst qlstx;
1027                         lst qlsts;
1028                         for (int j=r; j>=1; --j) {
1029                                 qlstx.append(1-newx.op(j-1));
1030                                 if (newx.op(j-1).info(info_flags::real) && newx.op(j-1) > 1 && newx.op(j-1) <= 2) {
1031                                         qlsts.append( s.op(j-1));
1032                                 } else {
1033                                         qlsts.append( -s.op(j-1));
1034                                 }
1035                         }
1036                         if (qlstx.nops() > 0) {
1037                                 buffer *= G_numeric(qlstx, qlsts, 1/q);
1038                         }
1039                         lst plstx;
1040                         lst plsts;
1041                         for (int j=r+1; j<=size; ++j) {
1042                                 plstx.append(newx.op(j-1));
1043                                 plsts.append(s.op(j-1));
1044                         }
1045                         if (plstx.nops() > 0) {
1046                                 buffer *= G_numeric(plstx, plsts, 1/p);
1047                         }
1048                         result += buffer;
1049                 }
1050                 return result;
1051         }
1052         
1053         // convergence transformation
1054         if (need_trafo) {
1055
1056                 // sort (|x|<->position) to determine indices
1057                 std::multimap<ex,int> sortmap;
1058                 int size = 0;
1059                 for (int i=0; i<x.nops(); ++i) {
1060                         if (!x[i].is_zero()) {
1061                                 sortmap.insert(std::pair<ex,int>(abs(x[i]), i));
1062                                 ++size;
1063                         }
1064                 }
1065                 // include upper limit (scale)
1066                 sortmap.insert(std::pair<ex,int>(abs(y), x.nops()));
1067
1068                 // generate missing dummy-symbols
1069                 int i = 1;
1070                 gsyms.clear();
1071                 gsyms.push_back(symbol("GSYMS_ERROR"));
1072                 ex lastentry;
1073                 for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
1074                         if (it != sortmap.begin()) {
1075                                 if (it->second < x.nops()) {
1076                                         if (x[it->second] == lastentry) {
1077                                                 gsyms.push_back(gsyms.back());
1078                                                 continue;
1079                                         }
1080                                 } else {
1081                                         if (y == lastentry) {
1082                                                 gsyms.push_back(gsyms.back());
1083                                                 continue;
1084                                         }
1085                                 }
1086                         }
1087                         std::ostringstream os;
1088                         os << "a" << i;
1089                         gsyms.push_back(symbol(os.str()));
1090                         ++i;
1091                         if (it->second < x.nops()) {
1092                                 lastentry = x[it->second];
1093                         } else {
1094                                 lastentry = y;
1095                         }
1096                 }
1097
1098                 // fill position data according to sorted indices and prepare substitution list
1099                 Gparameter a(x.nops());
1100                 lst subslst;
1101                 int pos = 1;
1102                 int scale;
1103                 for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
1104                         if (it->second < x.nops()) {
1105                                 if (s[it->second] > 0) {
1106                                         a[it->second] = pos;
1107                                 } else {
1108                                         a[it->second] = -pos;
1109                                 }
1110                                 subslst.append(gsyms[pos] == x[it->second]);
1111                         } else {
1112                                 scale = pos;
1113                                 subslst.append(gsyms[pos] == y);
1114                         }
1115                         ++pos;
1116                 }
1117
1118                 // do transformation
1119                 Gparameter pendint;
1120                 ex result = G_transform(pendint, a, scale);
1121                 // replace dummy symbols with their values
1122                 result = result.eval().expand();
1123                 result = result.subs(subslst).evalf();
1124                 
1125                 return result;
1126         }
1127
1128         // do summation
1129         lst newx;
1130         lst m;
1131         int mcount = 1;
1132         ex sign = 1;
1133         ex factor = y;
1134         for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
1135                 if ((*it).is_zero()) {
1136                         ++mcount;
1137                 } else {
1138                         newx.append(factor / (*it));
1139                         factor = *it;
1140                         m.append(mcount);
1141                         mcount = 1;
1142                         sign = -sign;
1143                 }
1144         }
1145
1146         return sign * numeric(mLi_do_summation(m, newx));
1147 }
1148
1149
1150 ex mLi_numeric(const lst& m, const lst& x)
1151 {
1152         // let G_numeric do the transformation
1153         lst newx;
1154         lst s;
1155         ex factor = 1;
1156         for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
1157                 for (int i = 1; i < *itm; ++i) {
1158                         newx.append(0);
1159                         s.append(1);
1160                 }
1161                 newx.append(factor / *itx);
1162                 factor /= *itx;
1163                 s.append(1);
1164         }
1165         return pow(-1, m.nops()) * G_numeric(newx, s, _ex1);
1166 }
1167
1168
1169 } // end of anonymous namespace
1170
1171
1172 //////////////////////////////////////////////////////////////////////
1173 //
1174 // Generalized multiple polylogarithm  G(x, y) and G(x, s, y)
1175 //
1176 // GiNaC function
1177 //
1178 //////////////////////////////////////////////////////////////////////
1179
1180
1181 static ex G2_evalf(const ex& x_, const ex& y)
1182 {
1183         if (!y.info(info_flags::positive)) {
1184                 return G(x_, y).hold();
1185         }
1186         lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
1187         if (x.nops() == 0) {
1188                 return _ex1;
1189         }
1190         if (x.op(0) == y) {
1191                 return G(x_, y).hold();
1192         }
1193         lst s;
1194         bool all_zero = true;
1195         for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
1196                 if (!(*it).info(info_flags::numeric)) {
1197                         return G(x_, y).hold();
1198                 }
1199                 if (*it != _ex0) {
1200                         all_zero = false;
1201                 }
1202                 s.append(+1);
1203         }
1204         if (all_zero) {
1205                 return pow(log(y), x.nops()) / factorial(x.nops());
1206         }
1207         return G_numeric(x, s, y);
1208 }
1209
1210
1211 static ex G2_eval(const ex& x_, const ex& y)
1212 {
1213         //TODO eval to MZV or H or S or Lin
1214
1215         if (!y.info(info_flags::positive)) {
1216                 return G(x_, y).hold();
1217         }
1218         lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
1219         if (x.nops() == 0) {
1220                 return _ex1;
1221         }
1222         if (x.op(0) == y) {
1223                 return G(x_, y).hold();
1224         }
1225         lst s;
1226         bool all_zero = true;
1227         bool crational = true;
1228         for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
1229                 if (!(*it).info(info_flags::numeric)) {
1230                         return G(x_, y).hold();
1231                 }
1232                 if (!(*it).info(info_flags::crational)) {
1233                         crational = false;
1234                 }
1235                 if (*it != _ex0) {
1236                         all_zero = false;
1237                 }
1238                 s.append(+1);
1239         }
1240         if (all_zero) {
1241                 return pow(log(y), x.nops()) / factorial(x.nops());
1242         }
1243         if (!y.info(info_flags::crational)) {
1244                 crational = false;
1245         }
1246         if (crational) {
1247                 return G(x_, y).hold();
1248         }
1249         return G_numeric(x, s, y);
1250 }
1251
1252
1253 unsigned G2_SERIAL::serial = function::register_new(function_options("G", 2).
1254                                 evalf_func(G2_evalf).
1255                                 eval_func(G2_eval).
1256                                 do_not_evalf_params().
1257                                 overloaded(2));
1258 //TODO
1259 //                                derivative_func(G2_deriv).
1260 //                                print_func<print_latex>(G2_print_latex).
1261
1262
1263 static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
1264 {
1265         if (!y.info(info_flags::positive)) {
1266                 return G(x_, s_, y).hold();
1267         }
1268         lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
1269         lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst(s_);
1270         if (x.nops() != s.nops()) {
1271                 return G(x_, s_, y).hold();
1272         }
1273         if (x.nops() == 0) {
1274                 return _ex1;
1275         }
1276         if (x.op(0) == y) {
1277                 return G(x_, s_, y).hold();
1278         }
1279         lst sn;
1280         bool all_zero = true;
1281         for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
1282                 if (!(*itx).info(info_flags::numeric)) {
1283                         return G(x_, y).hold();
1284                 }
1285                 if (!(*its).info(info_flags::real)) {
1286                         return G(x_, y).hold();
1287                 }
1288                 if (*itx != _ex0) {
1289                         all_zero = false;
1290                 }
1291                 if (*its >= 0) {
1292                         sn.append(+1);
1293                 } else {
1294                         sn.append(-1);
1295                 }
1296         }
1297         if (all_zero) {
1298                 return pow(log(y), x.nops()) / factorial(x.nops());
1299         }
1300         return G_numeric(x, sn, y);
1301 }
1302
1303
1304 static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
1305 {
1306         //TODO eval to MZV or H or S or Lin
1307
1308         if (!y.info(info_flags::positive)) {
1309                 return G(x_, s_, y).hold();
1310         }
1311         lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
1312         lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst(s_);
1313         if (x.nops() != s.nops()) {
1314                 return G(x_, s_, y).hold();
1315         }
1316         if (x.nops() == 0) {
1317                 return _ex1;
1318         }
1319         if (x.op(0) == y) {
1320                 return G(x_, s_, y).hold();
1321         }
1322         lst sn;
1323         bool all_zero = true;
1324         bool crational = true;
1325         for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
1326                 if (!(*itx).info(info_flags::numeric)) {
1327                         return G(x_, s_, y).hold();
1328                 }
1329                 if (!(*its).info(info_flags::real)) {
1330                         return G(x_, s_, y).hold();
1331                 }
1332                 if (!(*itx).info(info_flags::crational)) {
1333                         crational = false;
1334                 }
1335                 if (*itx != _ex0) {
1336                         all_zero = false;
1337                 }
1338                 if (*its >= 0) {
1339                         sn.append(+1);
1340                 } else {
1341                         sn.append(-1);
1342                 }
1343         }
1344         if (all_zero) {
1345                 return pow(log(y), x.nops()) / factorial(x.nops());
1346         }
1347         if (!y.info(info_flags::crational)) {
1348                 crational = false;
1349         }
1350         if (crational) {
1351                 return G(x_, s_, y).hold();
1352         }
1353         return G_numeric(x, sn, y);
1354 }
1355
1356
1357 unsigned G3_SERIAL::serial = function::register_new(function_options("G", 3).
1358                                 evalf_func(G3_evalf).
1359                                 eval_func(G3_eval).
1360                                 do_not_evalf_params().
1361                                 overloaded(2));
1362 //TODO
1363 //                                derivative_func(G3_deriv).
1364 //                                print_func<print_latex>(G3_print_latex).
1365
1366
1367 //////////////////////////////////////////////////////////////////////
1368 //
1369 // Classical polylogarithm and multiple polylogarithm  Li(m,x)
1370 //
1371 // GiNaC function
1372 //
1373 //////////////////////////////////////////////////////////////////////
1374
1375
1376 static ex Li_evalf(const ex& m_, const ex& x_)
1377 {
1378         // classical polylogs
1379         if (m_.info(info_flags::posint)) {
1380                 if (x_.info(info_flags::numeric)) {
1381                         int m__ = ex_to<numeric>(m_).to_int();
1382                         const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
1383                         const cln::cl_N result = Lin_numeric(m__, x__);
1384                         return numeric(result);
1385                 } else {
1386                         // try to numerically evaluate second argument
1387                         ex x_val = x_.evalf();
1388                         if (x_val.info(info_flags::numeric)) {
1389                                 int m__ = ex_to<numeric>(m_).to_int();
1390                                 const cln::cl_N x__ = ex_to<numeric>(x_val).to_cl_N();
1391                                 const cln::cl_N result = Lin_numeric(m__, x__);
1392                                 return numeric(result);
1393                         }
1394                 }
1395         }
1396         // multiple polylogs
1397         if (is_a<lst>(m_) && is_a<lst>(x_)) {
1398
1399                 const lst& m = ex_to<lst>(m_);
1400                 const lst& x = ex_to<lst>(x_);
1401                 if (m.nops() != x.nops()) {
1402                         return Li(m_,x_).hold();
1403                 }
1404                 if (x.nops() == 0) {
1405                         return _ex1;
1406                 }
1407                 if ((m.op(0) == _ex1) && (x.op(0) == _ex1)) {
1408                         return Li(m_,x_).hold();
1409                 }
1410
1411                 for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
1412                         if (!(*itm).info(info_flags::posint)) {
1413                                 return Li(m_, x_).hold();
1414                         }
1415                         if (!(*itx).info(info_flags::numeric)) {
1416                                 return Li(m_, x_).hold();
1417                         }
1418                         if (*itx == _ex0) {
1419                                 return _ex0;
1420                         }
1421                 }
1422
1423                 return mLi_numeric(m, x);
1424         }
1425
1426         return Li(m_,x_).hold();
1427 }
1428
1429
1430 static ex Li_eval(const ex& m_, const ex& x_)
1431 {
1432         if (is_a<lst>(m_)) {
1433                 if (is_a<lst>(x_)) {
1434                         // multiple polylogs
1435                         const lst& m = ex_to<lst>(m_);
1436                         const lst& x = ex_to<lst>(x_);
1437                         if (m.nops() != x.nops()) {
1438                                 return Li(m_,x_).hold();
1439                         }
1440                         if (x.nops() == 0) {
1441                                 return _ex1;
1442                         }
1443                         bool is_H = true;
1444                         bool is_zeta = true;
1445                         bool do_evalf = true;
1446                         bool crational = true;
1447                         for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
1448                                 if (!(*itm).info(info_flags::posint)) {
1449                                         return Li(m_,x_).hold();
1450                                 }
1451                                 if ((*itx != _ex1) && (*itx != _ex_1)) {
1452                                         if (itx != x.begin()) {
1453                                                 is_H = false;
1454                                         }
1455                                         is_zeta = false;
1456                                 }
1457                                 if (*itx == _ex0) {
1458                                         return _ex0;
1459                                 }
1460                                 if (!(*itx).info(info_flags::numeric)) {
1461                                         do_evalf = false;
1462                                 }
1463                                 if (!(*itx).info(info_flags::crational)) {
1464                                         crational = false;
1465                                 }
1466                         }
1467                         if (is_zeta) {
1468                                 return zeta(m_,x_);
1469                         }
1470                         if (is_H) {
1471                                 ex prefactor;
1472                                 lst newm = convert_parameter_Li_to_H(m, x, prefactor);
1473                                 return prefactor * H(newm, x[0]);
1474                         }
1475                         if (do_evalf && !crational) {
1476                                 return mLi_numeric(m,x);
1477                         }
1478                 }
1479                 return Li(m_, x_).hold();
1480         } else if (is_a<lst>(x_)) {
1481                 return Li(m_, x_).hold();
1482         }
1483
1484         // classical polylogs
1485         if (x_ == _ex0) {
1486                 return _ex0;
1487         }
1488         if (x_ == _ex1) {
1489                 return zeta(m_);
1490         }
1491         if (x_ == _ex_1) {
1492                 return (pow(2,1-m_)-1) * zeta(m_);
1493         }
1494         if (m_ == _ex1) {
1495                 return -log(1-x_);
1496         }
1497         if (m_ == _ex2) {
1498                 if (x_.is_equal(I)) {
1499                         return power(Pi,_ex2)/_ex_48 + Catalan*I;
1500                 }
1501                 if (x_.is_equal(-I)) {
1502                         return power(Pi,_ex2)/_ex_48 - Catalan*I;
1503                 }
1504         }
1505         if (m_.info(info_flags::posint) && x_.info(info_flags::numeric) && !x_.info(info_flags::crational)) {
1506                 int m__ = ex_to<numeric>(m_).to_int();
1507                 const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
1508                 const cln::cl_N result = Lin_numeric(m__, x__);
1509                 return numeric(result);
1510         }
1511
1512         return Li(m_, x_).hold();
1513 }
1514
1515
1516 static ex Li_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
1517 {
1518         if (is_a<lst>(m) || is_a<lst>(x)) {
1519                 // multiple polylog
1520                 epvector seq;
1521                 seq.push_back(expair(Li(m, x), 0));
1522                 return pseries(rel, seq);
1523         }
1524         
1525         // classical polylog
1526         const ex x_pt = x.subs(rel, subs_options::no_pattern);
1527         if (m.info(info_flags::numeric) && x_pt.info(info_flags::numeric)) {
1528                 // First special case: x==0 (derivatives have poles)
1529                 if (x_pt.is_zero()) {
1530                         const symbol s;
1531                         ex ser;
1532                         // manually construct the primitive expansion
1533                         for (int i=1; i<order; ++i)
1534                                 ser += pow(s,i) / pow(numeric(i), m);
1535                         // substitute the argument's series expansion
1536                         ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
1537                         // maybe that was terminating, so add a proper order term
1538                         epvector nseq;
1539                         nseq.push_back(expair(Order(_ex1), order));
1540                         ser += pseries(rel, nseq);
1541                         // reexpanding it will collapse the series again
1542                         return ser.series(rel, order);
1543                 }
1544                 // TODO special cases: x==1 (branch point) and x real, >=1 (branch cut)
1545                 throw std::runtime_error("Li_series: don't know how to do the series expansion at this point!");
1546         }
1547         // all other cases should be safe, by now:
1548         throw do_taylor();  // caught by function::series()
1549 }
1550
1551
1552 static ex Li_deriv(const ex& m_, const ex& x_, unsigned deriv_param)
1553 {
1554         GINAC_ASSERT(deriv_param < 2);
1555         if (deriv_param == 0) {
1556                 return _ex0;
1557         }
1558         if (m_.nops() > 1) {
1559                 throw std::runtime_error("don't know how to derivate multiple polylogarithm!");
1560         }
1561         ex m;
1562         if (is_a<lst>(m_)) {
1563                 m = m_.op(0);
1564         } else {
1565                 m = m_;
1566         }
1567         ex x;
1568         if (is_a<lst>(x_)) {
1569                 x = x_.op(0);
1570         } else {
1571                 x = x_;
1572         }
1573         if (m > 0) {
1574                 return Li(m-1, x) / x;
1575         } else {
1576                 return 1/(1-x);
1577         }
1578 }
1579
1580
1581 static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c)
1582 {
1583         lst m;
1584         if (is_a<lst>(m_)) {
1585                 m = ex_to<lst>(m_);
1586         } else {
1587                 m = lst(m_);
1588         }
1589         lst x;
1590         if (is_a<lst>(x_)) {
1591                 x = ex_to<lst>(x_);
1592         } else {
1593                 x = lst(x_);
1594         }
1595         c.s << "\\mbox{Li}_{";
1596         lst::const_iterator itm = m.begin();
1597         (*itm).print(c);
1598         itm++;
1599         for (; itm != m.end(); itm++) {
1600                 c.s << ",";
1601                 (*itm).print(c);
1602         }
1603         c.s << "}(";
1604         lst::const_iterator itx = x.begin();
1605         (*itx).print(c);
1606         itx++;
1607         for (; itx != x.end(); itx++) {
1608                 c.s << ",";
1609                 (*itx).print(c);
1610         }
1611         c.s << ")";
1612 }
1613
1614
1615 REGISTER_FUNCTION(Li,
1616                   evalf_func(Li_evalf).
1617                   eval_func(Li_eval).
1618                   series_func(Li_series).
1619                   derivative_func(Li_deriv).
1620                   print_func<print_latex>(Li_print_latex).
1621                   do_not_evalf_params());
1622
1623
1624 //////////////////////////////////////////////////////////////////////
1625 //
1626 // Nielsen's generalized polylogarithm  S(n,p,x)
1627 //
1628 // helper functions
1629 //
1630 //////////////////////////////////////////////////////////////////////
1631
1632
1633 // anonymous namespace for helper functions
1634 namespace {
1635
1636
1637 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
1638 // see fill_Yn()
1639 std::vector<std::vector<cln::cl_N> > Yn;
1640 int ynsize = 0; // number of Yn[]
1641 int ynlength = 100; // initial length of all Yn[i]
1642
1643
1644 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
1645 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
1646 // representing S_{n,p}(x).
1647 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
1648 // equivalent Z-sum.
1649 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
1650 // representing S_{n,p}(x).
1651 // The calculation of Y_n uses the values from Y_{n-1}.
1652 void fill_Yn(int n, const cln::float_format_t& prec)
1653 {
1654         const int initsize = ynlength;
1655         //const int initsize = initsize_Yn;
1656         cln::cl_N one = cln::cl_float(1, prec);
1657
1658         if (n) {
1659                 std::vector<cln::cl_N> buf(initsize);
1660                 std::vector<cln::cl_N>::iterator it = buf.begin();
1661                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
1662                 *it = (*itprev) / cln::cl_N(n+1) * one;
1663                 it++;
1664                 itprev++;
1665                 // sums with an index smaller than the depth are zero and need not to be calculated.
1666                 // calculation starts with depth, which is n+2)
1667                 for (int i=n+2; i<=initsize+n; i++) {
1668                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
1669                         it++;
1670                         itprev++;
1671                 }
1672                 Yn.push_back(buf);
1673         } else {
1674                 std::vector<cln::cl_N> buf(initsize);
1675                 std::vector<cln::cl_N>::iterator it = buf.begin();
1676                 *it = 1 * one;
1677                 it++;
1678                 for (int i=2; i<=initsize; i++) {
1679                         *it = *(it-1) + 1 / cln::cl_N(i) * one;
1680                         it++;
1681                 }
1682                 Yn.push_back(buf);
1683         }
1684         ynsize++;
1685 }
1686
1687
1688 // make Yn longer ... 
1689 void make_Yn_longer(int newsize, const cln::float_format_t& prec)
1690 {
1691
1692         cln::cl_N one = cln::cl_float(1, prec);
1693
1694         Yn[0].resize(newsize);
1695         std::vector<cln::cl_N>::iterator it = Yn[0].begin();
1696         it += ynlength;
1697         for (int i=ynlength+1; i<=newsize; i++) {
1698                 *it = *(it-1) + 1 / cln::cl_N(i) * one;
1699                 it++;
1700         }
1701
1702         for (int n=1; n<ynsize; n++) {
1703                 Yn[n].resize(newsize);
1704                 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
1705                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
1706                 it += ynlength;
1707                 itprev += ynlength;
1708                 for (int i=ynlength+n+1; i<=newsize+n; i++) {
1709                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
1710                         it++;
1711                         itprev++;
1712                 }
1713         }
1714         
1715         ynlength = newsize;
1716 }
1717
1718
1719 // helper function for S(n,p,x)
1720 // [Kol] (7.2)
1721 cln::cl_N C(int n, int p)
1722 {
1723         cln::cl_N result;
1724
1725         for (int k=0; k<p; k++) {
1726                 for (int j=0; j<=(n+k-1)/2; j++) {
1727                         if (k == 0) {
1728                                 if (n & 1) {
1729                                         if (j & 1) {
1730                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
1731                                         }
1732                                         else {
1733                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
1734                                         }
1735                                 }
1736                         }
1737                         else {
1738                                 if (k & 1) {
1739                                         if (j & 1) {
1740                                                 result = result + cln::factorial(n+k-1)
1741                                                                   * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1742                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1743                                         }
1744                                         else {
1745                                                 result = result - cln::factorial(n+k-1)
1746                                                                   * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1747                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1748                                         }
1749                                 }
1750                                 else {
1751                                         if (j & 1) {
1752                                                 result = result - cln::factorial(n+k-1) * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1753                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1754                                         }
1755                                         else {
1756                                                 result = result + cln::factorial(n+k-1)
1757                                                                   * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1758                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1759                                         }
1760                                 }
1761                         }
1762                 }
1763         }
1764         int np = n+p;
1765         if ((np-1) & 1) {
1766                 if (((np)/2+n) & 1) {
1767                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
1768                 }
1769                 else {
1770                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
1771                 }
1772         }
1773
1774         return result;
1775 }
1776
1777
1778 // helper function for S(n,p,x)
1779 // [Kol] remark to (9.1)
1780 cln::cl_N a_k(int k)
1781 {
1782         cln::cl_N result;
1783
1784         if (k == 0) {
1785                 return 1;
1786         }
1787
1788         result = result;
1789         for (int m=2; m<=k; m++) {
1790                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
1791         }
1792
1793         return -result / k;
1794 }
1795
1796
1797 // helper function for S(n,p,x)
1798 // [Kol] remark to (9.1)
1799 cln::cl_N b_k(int k)
1800 {
1801         cln::cl_N result;
1802
1803         if (k == 0) {
1804                 return 1;
1805         }
1806
1807         result = result;
1808         for (int m=2; m<=k; m++) {
1809                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
1810         }
1811
1812         return result / k;
1813 }
1814
1815
1816 // helper function for S(n,p,x)
1817 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
1818 {
1819         static cln::float_format_t oldprec = cln::default_float_format;
1820
1821         if (p==1) {
1822                 return Li_projection(n+1, x, prec);
1823         }
1824
1825         // precision has changed, we need to clear lookup table Yn
1826         if ( oldprec != prec ) {
1827                 Yn.clear();
1828                 ynsize = 0;
1829                 ynlength = 100;
1830                 oldprec = prec;
1831         }
1832                 
1833         // check if precalculated values are sufficient
1834         if (p > ynsize+1) {
1835                 for (int i=ynsize; i<p-1; i++) {
1836                         fill_Yn(i, prec);
1837                 }
1838         }
1839
1840         // should be done otherwise
1841         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1842         cln::cl_N xf = x * one;
1843         //cln::cl_N xf = x * cln::cl_float(1, prec);
1844
1845         cln::cl_N res;
1846         cln::cl_N resbuf;
1847         cln::cl_N factor = cln::expt(xf, p);
1848         int i = p;
1849         do {
1850                 resbuf = res;
1851                 if (i-p >= ynlength) {
1852                         // make Yn longer
1853                         make_Yn_longer(ynlength*2, prec);
1854                 }
1855                 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? ...
1856                 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
1857                 factor = factor * xf;
1858                 i++;
1859         } while (res != resbuf);
1860         
1861         return res;
1862 }
1863
1864
1865 // helper function for S(n,p,x)
1866 cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
1867 {
1868         // [Kol] (5.3)
1869         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
1870
1871                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
1872                                    * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
1873
1874                 for (int s=0; s<n; s++) {
1875                         cln::cl_N res2;
1876                         for (int r=0; r<p; r++) {
1877                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
1878                                               * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
1879                         }
1880                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
1881                 }
1882
1883                 return result;
1884         }
1885         
1886         return S_do_sum(n, p, x, prec);
1887 }
1888
1889
1890 // helper function for S(n,p,x)
1891 const cln::cl_N S_num(int n, int p, const cln::cl_N& x)
1892 {
1893         if (x == 1) {
1894                 if (n == 1) {
1895                     // [Kol] (2.22) with (2.21)
1896                         return cln::zeta(p+1);
1897                 }
1898
1899                 if (p == 1) {
1900                     // [Kol] (2.22)
1901                         return cln::zeta(n+1);
1902                 }
1903
1904                 // [Kol] (9.1)
1905                 cln::cl_N result;
1906                 for (int nu=0; nu<n; nu++) {
1907                         for (int rho=0; rho<=p; rho++) {
1908                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
1909                                                   * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
1910                         }
1911                 }
1912                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
1913
1914                 return result;
1915         }
1916         else if (x == -1) {
1917                 // [Kol] (2.22)
1918                 if (p == 1) {
1919                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
1920                 }
1921 //              throw std::runtime_error("don't know how to evaluate this function!");
1922         }
1923
1924         // what is the desired float format?
1925         // first guess: default format
1926         cln::float_format_t prec = cln::default_float_format;
1927         const cln::cl_N value = x;
1928         // second guess: the argument's format
1929         if (!instanceof(realpart(value), cln::cl_RA_ring))
1930                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
1931         else if (!instanceof(imagpart(value), cln::cl_RA_ring))
1932                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
1933
1934         // [Kol] (5.3)
1935         if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95))) {
1936
1937                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
1938                                    * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
1939
1940                 for (int s=0; s<n; s++) {
1941                         cln::cl_N res2;
1942                         for (int r=0; r<p; r++) {
1943                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
1944                                               * S_num(p-r,n-s,1-value) / cln::factorial(r);
1945                         }
1946                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
1947                 }
1948
1949                 return result;
1950                 
1951         }
1952         // [Kol] (5.12)
1953         if (cln::abs(value) > 1) {
1954                 
1955                 cln::cl_N result;
1956
1957                 for (int s=0; s<p; s++) {
1958                         for (int r=0; r<=s; r++) {
1959                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
1960                                                   / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
1961                                                   * S_num(n+s-r,p-s,cln::recip(value));
1962                         }
1963                 }
1964                 result = result * cln::expt(cln::cl_I(-1),n);
1965
1966                 cln::cl_N res2;
1967                 for (int r=0; r<n; r++) {
1968                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
1969                 }
1970                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
1971
1972                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
1973
1974                 return result;
1975         }
1976         else {
1977                 return S_projection(n, p, value, prec);
1978         }
1979 }
1980
1981
1982 } // end of anonymous namespace
1983
1984
1985 //////////////////////////////////////////////////////////////////////
1986 //
1987 // Nielsen's generalized polylogarithm  S(n,p,x)
1988 //
1989 // GiNaC function
1990 //
1991 //////////////////////////////////////////////////////////////////////
1992
1993
1994 static ex S_evalf(const ex& n, const ex& p, const ex& x)
1995 {
1996         if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
1997                 const int n_ = ex_to<numeric>(n).to_int();
1998                 const int p_ = ex_to<numeric>(p).to_int();
1999                 if (is_a<numeric>(x)) {
2000                         const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
2001                         const cln::cl_N result = S_num(n_, p_, x_);
2002                         return numeric(result);
2003                 } else {
2004                         ex x_val = x.evalf();
2005                         if (is_a<numeric>(x_val)) {
2006                                 const cln::cl_N x_val_ = ex_to<numeric>(x_val).to_cl_N();
2007                                 const cln::cl_N result = S_num(n_, p_, x_val_);
2008                                 return numeric(result);
2009                         }
2010                 }
2011         }
2012         return S(n, p, x).hold();
2013 }
2014
2015
2016 static ex S_eval(const ex& n, const ex& p, const ex& x)
2017 {
2018         if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
2019                 if (x == 0) {
2020                         return _ex0;
2021                 }
2022                 if (x == 1) {
2023                         lst m(n+1);
2024                         for (int i=ex_to<numeric>(p).to_int()-1; i>0; i--) {
2025                                 m.append(1);
2026                         }
2027                         return zeta(m);
2028                 }
2029                 if (p == 1) {
2030                         return Li(n+1, x);
2031                 }
2032                 if (x.info(info_flags::numeric) && (!x.info(info_flags::crational))) {
2033                         int n_ = ex_to<numeric>(n).to_int();
2034                         int p_ = ex_to<numeric>(p).to_int();
2035                         const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
2036                         const cln::cl_N result = S_num(n_, p_, x_);
2037                         return numeric(result);
2038                 }
2039         }
2040         if (n.is_zero()) {
2041                 // [Kol] (5.3)
2042                 return pow(-log(1-x), p) / factorial(p);
2043         }
2044         return S(n, p, x).hold();
2045 }
2046
2047
2048 static ex S_series(const ex& n, const ex& p, const ex& x, const relational& rel, int order, unsigned options)
2049 {
2050         if (p == _ex1) {
2051                 return Li(n+1, x).series(rel, order, options);
2052         }
2053
2054         const ex x_pt = x.subs(rel, subs_options::no_pattern);
2055         if (n.info(info_flags::posint) && p.info(info_flags::posint) && x_pt.info(info_flags::numeric)) {
2056                 // First special case: x==0 (derivatives have poles)
2057                 if (x_pt.is_zero()) {
2058                         const symbol s;
2059                         ex ser;
2060                         // manually construct the primitive expansion
2061                         // subsum = Euler-Zagier-Sum is needed
2062                         // dirty hack (slow ...) calculation of subsum:
2063                         std::vector<ex> presubsum, subsum;
2064                         subsum.push_back(0);
2065                         for (int i=1; i<order-1; ++i) {
2066                                 subsum.push_back(subsum[i-1] + numeric(1, i));
2067                         }
2068                         for (int depth=2; depth<p; ++depth) {
2069                                 presubsum = subsum;
2070                                 for (int i=1; i<order-1; ++i) {
2071                                         subsum[i] = subsum[i-1] + numeric(1, i) * presubsum[i-1];
2072                                 }
2073                         }
2074                                 
2075                         for (int i=1; i<order; ++i) {
2076                                 ser += pow(s,i) / pow(numeric(i), n+1) * subsum[i-1];
2077                         }
2078                         // substitute the argument's series expansion
2079                         ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
2080                         // maybe that was terminating, so add a proper order term
2081                         epvector nseq;
2082                         nseq.push_back(expair(Order(_ex1), order));
2083                         ser += pseries(rel, nseq);
2084                         // reexpanding it will collapse the series again
2085                         return ser.series(rel, order);
2086                 }
2087                 // TODO special cases: x==1 (branch point) and x real, >=1 (branch cut)
2088                 throw std::runtime_error("S_series: don't know how to do the series expansion at this point!");
2089         }
2090         // all other cases should be safe, by now:
2091         throw do_taylor();  // caught by function::series()
2092 }
2093
2094
2095 static ex S_deriv(const ex& n, const ex& p, const ex& x, unsigned deriv_param)
2096 {
2097         GINAC_ASSERT(deriv_param < 3);
2098         if (deriv_param < 2) {
2099                 return _ex0;
2100         }
2101         if (n > 0) {
2102                 return S(n-1, p, x) / x;
2103         } else {
2104                 return S(n, p-1, x) / (1-x);
2105         }
2106 }
2107
2108
2109 static void S_print_latex(const ex& n, const ex& p, const ex& x, const print_context& c)
2110 {
2111         c.s << "\\mbox{S}_{";
2112         n.print(c);
2113         c.s << ",";
2114         p.print(c);
2115         c.s << "}(";
2116         x.print(c);
2117         c.s << ")";
2118 }
2119
2120
2121 REGISTER_FUNCTION(S,
2122                   evalf_func(S_evalf).
2123                   eval_func(S_eval).
2124                   series_func(S_series).
2125                   derivative_func(S_deriv).
2126                   print_func<print_latex>(S_print_latex).
2127                   do_not_evalf_params());
2128
2129
2130 //////////////////////////////////////////////////////////////////////
2131 //
2132 // Harmonic polylogarithm  H(m,x)
2133 //
2134 // helper functions
2135 //
2136 //////////////////////////////////////////////////////////////////////
2137
2138
2139 // anonymous namespace for helper functions
2140 namespace {
2141
2142         
2143 // regulates the pole (used by 1/x-transformation)
2144 symbol H_polesign("IMSIGN");
2145
2146
2147 // convert parameters from H to Li representation
2148 // parameters are expected to be in expanded form, i.e. only 0, 1 and -1
2149 // returns true if some parameters are negative
2150 bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
2151 {
2152         // expand parameter list
2153         lst mexp;
2154         for (lst::const_iterator it = l.begin(); it != l.end(); it++) {
2155                 if (*it > 1) {
2156                         for (ex count=*it-1; count > 0; count--) {
2157                                 mexp.append(0);
2158                         }
2159                         mexp.append(1);
2160                 } else if (*it < -1) {
2161                         for (ex count=*it+1; count < 0; count++) {
2162                                 mexp.append(0);
2163                         }
2164                         mexp.append(-1);
2165                 } else {
2166                         mexp.append(*it);
2167                 }
2168         }
2169         
2170         ex signum = 1;
2171         pf = 1;
2172         bool has_negative_parameters = false;
2173         ex acc = 1;
2174         for (lst::const_iterator it = mexp.begin(); it != mexp.end(); it++) {
2175                 if (*it == 0) {
2176                         acc++;
2177                         continue;
2178                 }
2179                 if (*it > 0) {
2180                         m.append((*it+acc-1) * signum);
2181                 } else {
2182                         m.append((*it-acc+1) * signum);
2183                 }
2184                 acc = 1;
2185                 signum = *it;
2186                 pf *= *it;
2187                 if (pf < 0) {
2188                         has_negative_parameters = true;
2189                 }
2190         }
2191         if (has_negative_parameters) {
2192                 for (int i=0; i<m.nops(); i++) {
2193                         if (m.op(i) < 0) {
2194                                 m.let_op(i) = -m.op(i);
2195                                 s.append(-1);
2196                         } else {
2197                                 s.append(1);
2198                         }
2199                 }
2200         }
2201         
2202         return has_negative_parameters;
2203 }
2204
2205
2206 // recursivly transforms H to corresponding multiple polylogarithms
2207 struct map_trafo_H_convert_to_Li : public map_function
2208 {
2209         ex operator()(const ex& e)
2210         {
2211                 if (is_a<add>(e) || is_a<mul>(e)) {
2212                         return e.map(*this);
2213                 }
2214                 if (is_a<function>(e)) {
2215                         std::string name = ex_to<function>(e).get_name();
2216                         if (name == "H") {
2217                                 lst parameter;
2218                                 if (is_a<lst>(e.op(0))) {
2219                                                 parameter = ex_to<lst>(e.op(0));
2220                                 } else {
2221                                         parameter = lst(e.op(0));
2222                                 }
2223                                 ex arg = e.op(1);
2224
2225                                 lst m;
2226                                 lst s;
2227                                 ex pf;
2228                                 if (convert_parameter_H_to_Li(parameter, m, s, pf)) {
2229                                         s.let_op(0) = s.op(0) * arg;
2230                                         return pf * Li(m, s).hold();
2231                                 } else {
2232                                         for (int i=0; i<m.nops(); i++) {
2233                                                 s.append(1);
2234                                         }
2235                                         s.let_op(0) = s.op(0) * arg;
2236                                         return Li(m, s).hold();
2237                                 }
2238                         }
2239                 }
2240                 return e;
2241         }
2242 };
2243
2244
2245 // recursivly transforms H to corresponding zetas
2246 struct map_trafo_H_convert_to_zeta : public map_function
2247 {
2248         ex operator()(const ex& e)
2249         {
2250                 if (is_a<add>(e) || is_a<mul>(e)) {
2251                         return e.map(*this);
2252                 }
2253                 if (is_a<function>(e)) {
2254                         std::string name = ex_to<function>(e).get_name();
2255                         if (name == "H") {
2256                                 lst parameter;
2257                                 if (is_a<lst>(e.op(0))) {
2258                                                 parameter = ex_to<lst>(e.op(0));
2259                                 } else {
2260                                         parameter = lst(e.op(0));
2261                                 }
2262
2263                                 lst m;
2264                                 lst s;
2265                                 ex pf;
2266                                 if (convert_parameter_H_to_Li(parameter, m, s, pf)) {
2267                                         return pf * zeta(m, s);
2268                                 } else {
2269                                         return zeta(m);
2270                                 }
2271                         }
2272                 }
2273                 return e;
2274         }
2275 };
2276
2277
2278 // remove trailing zeros from H-parameters
2279 struct map_trafo_H_reduce_trailing_zeros : public map_function
2280 {
2281         ex operator()(const ex& e)
2282         {
2283                 if (is_a<add>(e) || is_a<mul>(e)) {
2284                         return e.map(*this);
2285                 }
2286                 if (is_a<function>(e)) {
2287                         std::string name = ex_to<function>(e).get_name();
2288                         if (name == "H") {
2289                                 lst parameter;
2290                                 if (is_a<lst>(e.op(0))) {
2291                                         parameter = ex_to<lst>(e.op(0));
2292                                 } else {
2293                                         parameter = lst(e.op(0));
2294                                 }
2295                                 ex arg = e.op(1);
2296                                 if (parameter.op(parameter.nops()-1) == 0) {
2297                                         
2298                                         //
2299                                         if (parameter.nops() == 1) {
2300                                                 return log(arg);
2301                                         }
2302                                         
2303                                         //
2304                                         lst::const_iterator it = parameter.begin();
2305                                         while ((it != parameter.end()) && (*it == 0)) {
2306                                                 it++;
2307                                         }
2308                                         if (it == parameter.end()) {
2309                                                 return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
2310                                         }
2311                                         
2312                                         //
2313                                         parameter.remove_last();
2314                                         int lastentry = parameter.nops();
2315                                         while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
2316                                                 lastentry--;
2317                                         }
2318                                         
2319                                         //
2320                                         ex result = log(arg) * H(parameter,arg).hold();
2321                                         ex acc = 0;
2322                                         for (ex i=0; i<lastentry; i++) {
2323                                                 if (parameter[i] > 0) {
2324                                                         parameter[i]++;
2325                                                         result -= (acc + parameter[i]-1) * H(parameter, arg).hold();
2326                                                         parameter[i]--;
2327                                                         acc = 0;
2328                                                 } else if (parameter[i] < 0) {
2329                                                         parameter[i]--;
2330                                                         result -= (acc + abs(parameter[i]+1)) * H(parameter, arg).hold();
2331                                                         parameter[i]++;
2332                                                         acc = 0;
2333                                                 } else {
2334                                                         acc++;
2335                                                 }
2336                                         }
2337                                         
2338                                         if (lastentry < parameter.nops()) {
2339                                                 result = result / (parameter.nops()-lastentry+1);
2340                                                 return result.map(*this);
2341                                         } else {
2342                                                 return result;
2343                                         }
2344                                 }
2345                         }
2346                 }
2347                 return e;
2348         }
2349 };
2350
2351
2352 // returns an expression with zeta functions corresponding to the parameter list for H
2353 ex convert_H_to_zeta(const lst& m)
2354 {
2355         symbol xtemp("xtemp");
2356         map_trafo_H_reduce_trailing_zeros filter;
2357         map_trafo_H_convert_to_zeta filter2;
2358         return filter2(filter(H(m, xtemp).hold())).subs(xtemp == 1);
2359 }
2360
2361
2362 // convert signs form Li to H representation
2363 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf)
2364 {
2365         lst res;
2366         lst::const_iterator itm = m.begin();
2367         lst::const_iterator itx = ++x.begin();
2368         int signum = 1;
2369         pf = _ex1;
2370         res.append(*itm);
2371         itm++;
2372         while (itx != x.end()) {
2373                 signum *= (*itx > 0) ? 1 : -1;
2374                 pf *= signum;
2375                 res.append((*itm) * signum);
2376                 itm++;
2377                 itx++;
2378         }
2379         return res;
2380 }
2381
2382
2383 // multiplies an one-dimensional H with another H
2384 // [ReV] (18)
2385 ex trafo_H_mult(const ex& h1, const ex& h2)
2386 {
2387         ex res;
2388         ex hshort;
2389         lst hlong;
2390         ex h1nops = h1.op(0).nops();
2391         ex h2nops = h2.op(0).nops();
2392         if (h1nops > 1) {
2393                 hshort = h2.op(0).op(0);
2394                 hlong = ex_to<lst>(h1.op(0));
2395         } else {
2396                 hshort = h1.op(0).op(0);
2397                 if (h2nops > 1) {
2398                         hlong = ex_to<lst>(h2.op(0));
2399                 } else {
2400                         hlong = h2.op(0).op(0);
2401                 }
2402         }
2403         for (int i=0; i<=hlong.nops(); i++) {
2404                 lst newparameter;
2405                 int j=0;
2406                 for (; j<i; j++) {
2407                         newparameter.append(hlong[j]);
2408                 }
2409                 newparameter.append(hshort);
2410                 for (; j<hlong.nops(); j++) {
2411                         newparameter.append(hlong[j]);
2412                 }
2413                 res += H(newparameter, h1.op(1)).hold();
2414         }
2415         return res;
2416 }
2417
2418
2419 // applies trafo_H_mult recursively on expressions
2420 struct map_trafo_H_mult : public map_function
2421 {
2422         ex operator()(const ex& e)
2423         {
2424                 if (is_a<add>(e)) {
2425                         return e.map(*this);
2426                 }
2427
2428                 if (is_a<mul>(e)) {
2429
2430                         ex result = 1;
2431                         ex firstH;
2432                         lst Hlst;
2433                         for (int pos=0; pos<e.nops(); pos++) {
2434                                 if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
2435                                         std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
2436                                         if (name == "H") {
2437                                                 for (ex i=0; i<e.op(pos).op(1); i++) {
2438                                                         Hlst.append(e.op(pos).op(0));
2439                                                 }
2440                                                 continue;
2441                                         }
2442                                 } else if (is_a<function>(e.op(pos))) {
2443                                         std::string name = ex_to<function>(e.op(pos)).get_name();
2444                                         if (name == "H") {
2445                                                 if (e.op(pos).op(0).nops() > 1) {
2446                                                         firstH = e.op(pos);
2447                                                 } else {
2448                                                         Hlst.append(e.op(pos));
2449                                                 }
2450                                                 continue;
2451                                         }
2452                                 }
2453                                 result *= e.op(pos);
2454                         }
2455                         if (firstH == 0) {
2456                                 if (Hlst.nops() > 0) {
2457                                         firstH = Hlst[Hlst.nops()-1];
2458                                         Hlst.remove_last();
2459                                 } else {
2460                                         return e;
2461                                 }
2462                         }
2463
2464                         if (Hlst.nops() > 0) {
2465                                 ex buffer = trafo_H_mult(firstH, Hlst.op(0));
2466                                 result *= buffer;
2467                                 for (int i=1; i<Hlst.nops(); i++) {
2468                                         result *= Hlst.op(i);
2469                                 }
2470                                 result = result.expand();
2471                                 map_trafo_H_mult recursion;
2472                                 return recursion(result);
2473                         } else {
2474                                 return e;
2475                         }
2476
2477                 }
2478                 return e;
2479         }
2480 };
2481
2482
2483 // do integration [ReV] (55)
2484 // put parameter 0 in front of existing parameters
2485 ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
2486 {
2487         ex h;
2488         std::string name;
2489         if (is_a<function>(e)) {
2490                 name = ex_to<function>(e).get_name();
2491         }
2492         if (name == "H") {
2493                 h = e;
2494         } else {
2495                 for (int i=0; i<e.nops(); i++) {
2496                         if (is_a<function>(e.op(i))) {
2497                                 std::string name = ex_to<function>(e.op(i)).get_name();
2498                                 if (name == "H") {
2499                                         h = e.op(i);
2500                                 }
2501                         }
2502                 }
2503         }
2504         if (h != 0) {
2505                 lst newparameter = ex_to<lst>(h.op(0));
2506                 newparameter.prepend(0);
2507                 ex addzeta = convert_H_to_zeta(newparameter);
2508                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
2509         } else {
2510                 return e * (-H(lst(0),1/arg).hold());
2511         }
2512 }
2513
2514
2515 // do integration [ReV] (49)
2516 // put parameter 1 in front of existing parameters
2517 ex trafo_H_prepend_one(const ex& e, const ex& arg)
2518 {
2519         ex h;
2520         std::string name;
2521         if (is_a<function>(e)) {
2522                 name = ex_to<function>(e).get_name();
2523         }
2524         if (name == "H") {
2525                 h = e;
2526         } else {
2527                 for (int i=0; i<e.nops(); i++) {
2528                         if (is_a<function>(e.op(i))) {
2529                                 std::string name = ex_to<function>(e.op(i)).get_name();
2530                                 if (name == "H") {
2531                                         h = e.op(i);
2532                                 }
2533                         }
2534                 }
2535         }
2536         if (h != 0) {
2537                 lst newparameter = ex_to<lst>(h.op(0));
2538                 newparameter.prepend(1);
2539                 return e.subs(h == H(newparameter, h.op(1)).hold());
2540         } else {
2541                 return e * H(lst(1),1-arg).hold();
2542         }
2543 }
2544
2545
2546 // do integration [ReV] (55)
2547 // put parameter -1 in front of existing parameters
2548 ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
2549 {
2550         ex h;
2551         std::string name;
2552         if (is_a<function>(e)) {
2553                 name = ex_to<function>(e).get_name();
2554         }
2555         if (name == "H") {
2556                 h = e;
2557         } else {
2558                 for (int i=0; i<e.nops(); i++) {
2559                         if (is_a<function>(e.op(i))) {
2560                                 std::string name = ex_to<function>(e.op(i)).get_name();
2561                                 if (name == "H") {
2562                                         h = e.op(i);
2563                                 }
2564                         }
2565                 }
2566         }
2567         if (h != 0) {
2568                 lst newparameter = ex_to<lst>(h.op(0));
2569                 newparameter.prepend(-1);
2570                 ex addzeta = convert_H_to_zeta(newparameter);
2571                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
2572         } else {
2573                 ex addzeta = convert_H_to_zeta(lst(-1));
2574                 return (e * (addzeta - H(lst(-1),1/arg).hold())).expand();
2575         }
2576 }
2577
2578
2579 // do integration [ReV] (55)
2580 // put parameter -1 in front of existing parameters
2581 ex trafo_H_1mxt1px_prepend_minusone(const ex& e, const ex& arg)
2582 {
2583         ex h;
2584         std::string name;
2585         if (is_a<function>(e)) {
2586                 name = ex_to<function>(e).get_name();
2587         }
2588         if (name == "H") {
2589                 h = e;
2590         } else {
2591                 for (int i=0; i<e.nops(); i++) {
2592                         if (is_a<function>(e.op(i))) {
2593                                 std::string name = ex_to<function>(e.op(i)).get_name();
2594                                 if (name == "H") {
2595                                         h = e.op(i);
2596                                 }
2597                         }
2598                 }
2599         }
2600         if (h != 0) {
2601                 lst newparameter = ex_to<lst>(h.op(0));
2602                 newparameter.prepend(-1);
2603                 return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
2604         } else {
2605                 return (e * H(lst(-1),(1-arg)/(1+arg)).hold()).expand();
2606         }
2607 }
2608
2609
2610 // do integration [ReV] (55)
2611 // put parameter 1 in front of existing parameters
2612 ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
2613 {
2614         ex h;
2615         std::string name;
2616         if (is_a<function>(e)) {
2617                 name = ex_to<function>(e).get_name();
2618         }
2619         if (name == "H") {
2620                 h = e;
2621         } else {
2622                 for (int i=0; i<e.nops(); i++) {
2623                         if (is_a<function>(e.op(i))) {
2624                                 std::string name = ex_to<function>(e.op(i)).get_name();
2625                                 if (name == "H") {
2626                                         h = e.op(i);
2627                                 }
2628                         }
2629                 }
2630         }
2631         if (h != 0) {
2632                 lst newparameter = ex_to<lst>(h.op(0));
2633                 newparameter.prepend(1);
2634                 return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
2635         } else {
2636                 return (e * H(lst(1),(1-arg)/(1+arg)).hold()).expand();
2637         }
2638 }
2639
2640
2641 // do x -> 1-x transformation
2642 struct map_trafo_H_1mx : public map_function
2643 {
2644         ex operator()(const ex& e)
2645         {
2646                 if (is_a<add>(e) || is_a<mul>(e)) {
2647                         return e.map(*this);
2648                 }
2649                 
2650                 if (is_a<function>(e)) {
2651                         std::string name = ex_to<function>(e).get_name();
2652                         if (name == "H") {
2653
2654                                 lst parameter = ex_to<lst>(e.op(0));
2655                                 ex arg = e.op(1);
2656
2657                                 // special cases if all parameters are either 0, 1 or -1
2658                                 bool allthesame = true;
2659                                 if (parameter.op(0) == 0) {
2660                                         for (int i=1; i<parameter.nops(); i++) {
2661                                                 if (parameter.op(i) != 0) {
2662                                                         allthesame = false;
2663                                                         break;
2664                                                 }
2665                                         }
2666                                         if (allthesame) {
2667                                                 lst newparameter;
2668                                                 for (int i=parameter.nops(); i>0; i--) {
2669                                                         newparameter.append(1);
2670                                                 }
2671                                                 return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
2672                                         }
2673                                 } else if (parameter.op(0) == -1) {
2674                                         throw std::runtime_error("map_trafo_H_1mx: cannot handle weights equal -1!");
2675                                 } else {
2676                                         for (int i=1; i<parameter.nops(); i++) {
2677                                                 if (parameter.op(i) != 1) {
2678                                                         allthesame = false;
2679                                                         break;
2680                                                 }
2681                                         }
2682                                         if (allthesame) {
2683                                                 lst newparameter;
2684                                                 for (int i=parameter.nops(); i>0; i--) {
2685                                                         newparameter.append(0);
2686                                                 }
2687                                                 return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
2688                                         }
2689                                 }
2690
2691                                 lst newparameter = parameter;
2692                                 newparameter.remove_first();
2693
2694                                 if (parameter.op(0) == 0) {
2695
2696                                         // leading zero
2697                                         ex res = convert_H_to_zeta(parameter);
2698                                         //ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
2699                                         map_trafo_H_1mx recursion;
2700                                         ex buffer = recursion(H(newparameter, arg).hold());
2701                                         if (is_a<add>(buffer)) {
2702                                                 for (int i=0; i<buffer.nops(); i++) {
2703                                                         res -= trafo_H_prepend_one(buffer.op(i), arg);
2704                                                 }
2705                                         } else {
2706                                                 res -= trafo_H_prepend_one(buffer, arg);
2707                                         }
2708                                         return res;
2709
2710                                 } else {
2711
2712                                         // leading one
2713                                         map_trafo_H_1mx recursion;
2714                                         map_trafo_H_mult unify;
2715                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
2716                                         int firstzero = 0;
2717                                         while (parameter.op(firstzero) == 1) {
2718                                                 firstzero++;
2719                                         }
2720                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
2721                                                 lst newparameter;
2722                                                 int j=0;
2723                                                 for (; j<=i; j++) {
2724                                                         newparameter.append(parameter[j+1]);
2725                                                 }
2726                                                 newparameter.append(1);
2727                                                 for (; j<parameter.nops()-1; j++) {
2728                                                         newparameter.append(parameter[j+1]);
2729                                                 }
2730                                                 res -= H(newparameter, arg).hold();
2731                                         }
2732                                         res = recursion(res).expand() / firstzero;
2733                                         return unify(res);
2734                                 }
2735                         }
2736                 }
2737                 return e;
2738         }
2739 };
2740
2741
2742 // do x -> 1/x transformation
2743 struct map_trafo_H_1overx : public map_function
2744 {
2745         ex operator()(const ex& e)
2746         {
2747                 if (is_a<add>(e) || is_a<mul>(e)) {
2748                         return e.map(*this);
2749                 }
2750
2751                 if (is_a<function>(e)) {
2752                         std::string name = ex_to<function>(e).get_name();
2753                         if (name == "H") {
2754
2755                                 lst parameter = ex_to<lst>(e.op(0));
2756                                 ex arg = e.op(1);
2757
2758                                 // special cases if all parameters are either 0, 1 or -1
2759                                 bool allthesame = true;
2760                                 if (parameter.op(0) == 0) {
2761                                         for (int i=1; i<parameter.nops(); i++) {
2762                                                 if (parameter.op(i) != 0) {
2763                                                         allthesame = false;
2764                                                         break;
2765                                                 }
2766                                         }
2767                                         if (allthesame) {
2768                                                 return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
2769                                         }
2770                                 } else if (parameter.op(0) == -1) {
2771                                         for (int i=1; i<parameter.nops(); i++) {
2772                                                 if (parameter.op(i) != -1) {
2773                                                         allthesame = false;
2774                                                         break;
2775                                                 }
2776                                         }
2777                                         if (allthesame) {
2778                                                 map_trafo_H_mult unify;
2779                                                 return unify((pow(H(lst(-1),1/arg).hold() - H(lst(0),1/arg).hold(), parameter.nops())
2780                                                        / factorial(parameter.nops())).expand());
2781                                         }
2782                                 } else {
2783                                         for (int i=1; i<parameter.nops(); i++) {
2784                                                 if (parameter.op(i) != 1) {
2785                                                         allthesame = false;
2786                                                         break;
2787                                                 }
2788                                         }
2789                                         if (allthesame) {
2790                                                 map_trafo_H_mult unify;
2791                                                 return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() + H_polesign, parameter.nops())
2792                                                        / factorial(parameter.nops())).expand());
2793                                         }
2794                                 }
2795
2796                                 lst newparameter = parameter;
2797                                 newparameter.remove_first();
2798
2799                                 if (parameter.op(0) == 0) {
2800                                         
2801                                         // leading zero
2802                                         ex res = convert_H_to_zeta(parameter);
2803                                         map_trafo_H_1overx recursion;
2804                                         ex buffer = recursion(H(newparameter, arg).hold());
2805                                         if (is_a<add>(buffer)) {
2806                                                 for (int i=0; i<buffer.nops(); i++) {
2807                                                         res += trafo_H_1tx_prepend_zero(buffer.op(i), arg);
2808                                                 }
2809                                         } else {
2810                                                 res += trafo_H_1tx_prepend_zero(buffer, arg);
2811                                         }
2812                                         return res;
2813
2814                                 } else if (parameter.op(0) == -1) {
2815
2816                                         // leading negative one
2817                                         ex res = convert_H_to_zeta(parameter);
2818                                         map_trafo_H_1overx recursion;
2819                                         ex buffer = recursion(H(newparameter, arg).hold());
2820                                         if (is_a<add>(buffer)) {
2821                                                 for (int i=0; i<buffer.nops(); i++) {
2822                                                         res += trafo_H_1tx_prepend_zero(buffer.op(i), arg) - trafo_H_1tx_prepend_minusone(buffer.op(i), arg);
2823                                                 }
2824                                         } else {
2825                                                 res += trafo_H_1tx_prepend_zero(buffer, arg) - trafo_H_1tx_prepend_minusone(buffer, arg);
2826                                         }
2827                                         return res;
2828
2829                                 } else {
2830
2831                                         // leading one
2832                                         map_trafo_H_1overx recursion;
2833                                         map_trafo_H_mult unify;
2834                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
2835                                         int firstzero = 0;
2836                                         while (parameter.op(firstzero) == 1) {
2837                                                 firstzero++;
2838                                         }
2839                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
2840                                                 lst newparameter;
2841                                                 int j=0;
2842                                                 for (; j<=i; j++) {
2843                                                         newparameter.append(parameter[j+1]);
2844                                                 }
2845                                                 newparameter.append(1);
2846                                                 for (; j<parameter.nops()-1; j++) {
2847                                                         newparameter.append(parameter[j+1]);
2848                                                 }
2849                                                 res -= H(newparameter, arg).hold();
2850                                         }
2851                                         res = recursion(res).expand() / firstzero;
2852                                         return unify(res);
2853
2854                                 }
2855
2856                         }
2857                 }
2858                 return e;
2859         }
2860 };
2861
2862
2863 // do x -> (1-x)/(1+x) transformation
2864 struct map_trafo_H_1mxt1px : public map_function
2865 {
2866         ex operator()(const ex& e)
2867         {
2868                 if (is_a<add>(e) || is_a<mul>(e)) {
2869                         return e.map(*this);
2870                 }
2871
2872                 if (is_a<function>(e)) {
2873                         std::string name = ex_to<function>(e).get_name();
2874                         if (name == "H") {
2875
2876                                 lst parameter = ex_to<lst>(e.op(0));
2877                                 ex arg = e.op(1);
2878
2879                                 // special cases if all parameters are either 0, 1 or -1
2880                                 bool allthesame = true;
2881                                 if (parameter.op(0) == 0) {
2882                                         for (int i=1; i<parameter.nops(); i++) {
2883                                                 if (parameter.op(i) != 0) {
2884                                                         allthesame = false;
2885                                                         break;
2886                                                 }
2887                                         }
2888                                         if (allthesame) {
2889                                                 map_trafo_H_mult unify;
2890                                                 return unify((pow(-H(lst(1),(1-arg)/(1+arg)).hold() - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
2891                                                        / factorial(parameter.nops())).expand());
2892                                         }
2893                                 } else if (parameter.op(0) == -1) {
2894                                         for (int i=1; i<parameter.nops(); i++) {
2895                                                 if (parameter.op(i) != -1) {
2896                                                         allthesame = false;
2897                                                         break;
2898                                                 }
2899                                         }
2900                                         if (allthesame) {
2901                                                 map_trafo_H_mult unify;
2902                                                 return unify((pow(log(2) - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
2903                                                        / factorial(parameter.nops())).expand());
2904                                         }
2905                                 } else {
2906                                         for (int i=1; i<parameter.nops(); i++) {
2907                                                 if (parameter.op(i) != 1) {
2908                                                         allthesame = false;
2909                                                         break;
2910                                                 }
2911                                         }
2912                                         if (allthesame) {
2913                                                 map_trafo_H_mult unify;
2914                                                 return unify((pow(-log(2) - H(lst(0),(1-arg)/(1+arg)).hold() + H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
2915                                                        / factorial(parameter.nops())).expand());
2916                                         }
2917                                 }
2918
2919                                 lst newparameter = parameter;
2920                                 newparameter.remove_first();
2921
2922                                 if (parameter.op(0) == 0) {
2923
2924                                         // leading zero
2925                                         ex res = convert_H_to_zeta(parameter);
2926                                         map_trafo_H_1mxt1px recursion;
2927                                         ex buffer = recursion(H(newparameter, arg).hold());
2928                                         if (is_a<add>(buffer)) {
2929                                                 for (int i=0; i<buffer.nops(); i++) {
2930                                                         res -= trafo_H_1mxt1px_prepend_one(buffer.op(i), arg) + trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
2931                                                 }
2932                                         } else {
2933                                                 res -= trafo_H_1mxt1px_prepend_one(buffer, arg) + trafo_H_1mxt1px_prepend_minusone(buffer, arg);
2934                                         }
2935                                         return res;
2936
2937                                 } else if (parameter.op(0) == -1) {
2938
2939                                         // leading negative one
2940                                         ex res = convert_H_to_zeta(parameter);
2941                                         map_trafo_H_1mxt1px recursion;
2942                                         ex buffer = recursion(H(newparameter, arg).hold());
2943                                         if (is_a<add>(buffer)) {
2944                                                 for (int i=0; i<buffer.nops(); i++) {
2945                                                         res -= trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
2946                                                 }
2947                                         } else {
2948                                   &n