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