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