Synced to 1.1
[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  *  The functions are: 
5  *    classical polylogarithm              Li(n,x)
6  *    multiple polylogarithm               Li(lst(n_1,...,n_k),lst(x_1,...,x_k)
7  *    nielsen's generalized polylogarithm  S(n,p,x)
8  *    harmonic polylogarithm               H(lst(m_1,...,m_k),x)
9  *    multiple zeta value                  mZeta(lst(m_1,...,m_k))
10  *
11  *  Some remarks:
12  *    - All formulae used can be looked up in the following publication:
13  *      Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
14  *      This document will be referenced as [Kol] throughout this source code.
15  *    - Classical polylogarithms (Li) and nielsen's generalized polylogarithms (S) can be numerically
16  *      evaluated in the whole complex plane. And of course, there is still room for speed optimizations ;-).
17  *    - The calculation of classical polylogarithms is speed up by using Euler-Maclaurin summation (EuMac).
18  *    - The remaining functions can only be numerically evaluated with arguments lying in the unit sphere
19  *      at the moment. Sorry. The evaluation especially for mZeta is very slow ... better not use it
20  *      right now.
21  *    - The functions have no series expansion. To do it, you have to convert these functions
22  *      into the appropriate objects from the nestedsums library, do the expansion and convert the
23  *      result back. 
24  *    - Numerical testing of this implementation has been performed by doing a comparison of results
25  *      between this software and the commercial M.......... 4.1.
26  *
27  */
28
29 /*
30  *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
31  *
32  *  This program is free software; you can redistribute it and/or modify
33  *  it under the terms of the GNU General Public License as published by
34  *  the Free Software Foundation; either version 2 of the License, or
35  *  (at your option) any later version.
36  *
37  *  This program is distributed in the hope that it will be useful,
38  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
39  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
40  *  GNU General Public License for more details.
41  *
42  *  You should have received a copy of the GNU General Public License
43  *  along with this program; if not, write to the Free Software
44  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
45  */
46
47 #include <stdexcept>
48 #include <vector>
49 #include <cln/cln.h>
50
51 #include "inifcns.h"
52 #include "lst.h"
53 #include "numeric.h"
54 #include "operators.h"
55 #include "relational.h"
56 #include "pseries.h"
57
58
59 namespace GiNaC {
60
61         
62 // lookup table for Euler-Maclaurin optimization
63 // see fill_Xn()
64 std::vector<std::vector<cln::cl_N> > Xn;
65 int xnsize = 0;
66
67
68 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
69 // see fill_Yn()
70 std::vector<std::vector<cln::cl_N> > Yn;
71 int ynsize = 0; // number of Yn[]
72 int ynlength = 100; // initial length of all Yn[i]
73
74
75 //////////////////////
76 // helper functions //
77 //////////////////////
78
79
80 // This function calculates the X_n. The X_n are needed for the Euler-Maclaurin summation (EMS) of
81 // classical polylogarithms.
82 // With EMS the polylogs can be calculated as follows:
83 //   Li_p (x)  =  \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with  u = -log(1-x)
84 //   X_0(n) = B_n (Bernoulli numbers)
85 //   X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
86 // The calculation of Xn depends on X0 and X{n-1}.
87 // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
88 // This results in a slightly more complicated algorithm for the X_n.
89 // The first index in Xn corresponds to the index of the polylog minus 2.
90 // The second index in Xn corresponds to the index from the EMS.
91 static void fill_Xn(int n)
92 {
93         // rule of thumb. needs to be improved. TODO
94         const int initsize = Digits * 3 / 2;
95
96         if (n>1) {
97                 // calculate X_2 and higher (corresponding to Li_4 and higher)
98                 std::vector<cln::cl_N> buf(initsize);
99                 std::vector<cln::cl_N>::iterator it = buf.begin();
100                 cln::cl_N result;
101                 *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
102                 it++;
103                 for (int i=2; i<=initsize; i++) {
104                         if (i&1) {
105                                 result = 0; // k == 0
106                         } else {
107                                 result = Xn[0][i/2-1]; // k == 0
108                         }
109                         for (int k=1; k<i-1; k++) {
110                                 if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
111                                         result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
112                                 }
113                         }
114                         result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
115                         result = result + Xn[n-1][i-1] / (i+1); // k == i
116                         
117                         *it = result;
118                         it++;
119                 }
120                 Xn.push_back(buf);
121         } else if (n==1) {
122                 // special case to handle the X_0 correct
123                 std::vector<cln::cl_N> buf(initsize);
124                 std::vector<cln::cl_N>::iterator it = buf.begin();
125                 cln::cl_N result;
126                 *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
127                 it++;
128                 *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
129                 it++;
130                 for (int i=3; i<=initsize; i++) {
131                         if (i & 1) {
132                                 result = -Xn[0][(i-3)/2]/2;
133                                 *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
134                                 it++;
135                         } else {
136                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
137                                 for (int k=1; k<i/2; k++) {
138                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
139                                 }
140                                 *it = result;
141                                 it++;
142                         }
143                 }
144                 Xn.push_back(buf);
145         } else {
146                 // calculate X_0
147                 std::vector<cln::cl_N> buf(initsize/2);
148                 std::vector<cln::cl_N>::iterator it = buf.begin();
149                 for (int i=1; i<=initsize/2; i++) {
150                         *it = bernoulli(i*2).to_cl_N();
151                         it++;
152                 }
153                 Xn.push_back(buf);
154         }
155
156         xnsize++;
157 }
158
159
160 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
161 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
162 // representing S_{n,p}(x).
163 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
164 // equivalent Z-sum.
165 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
166 // representing S_{n,p}(x).
167 // The calculation of Y_n uses the values from Y_{n-1}.
168 static void fill_Yn(int n, const cln::float_format_t& prec)
169 {
170         const int initsize = ynlength;
171         //const int initsize = initsize_Yn;
172         cln::cl_N one = cln::cl_float(1, prec);
173
174         if (n) {
175                 std::vector<cln::cl_N> buf(initsize);
176                 std::vector<cln::cl_N>::iterator it = buf.begin();
177                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
178                 *it = (*itprev) / cln::cl_N(n+1) * one;
179                 it++;
180                 itprev++;
181                 // sums with an index smaller than the depth are zero and need not to be calculated.
182                 // calculation starts with depth, which is n+2)
183                 for (int i=n+2; i<=initsize+n; i++) {
184                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
185                         it++;
186                         itprev++;
187                 }
188                 Yn.push_back(buf);
189         } else {
190                 std::vector<cln::cl_N> buf(initsize);
191                 std::vector<cln::cl_N>::iterator it = buf.begin();
192                 *it = 1 * one;
193                 it++;
194                 for (int i=2; i<=initsize; i++) {
195                         *it = *(it-1) + 1 / cln::cl_N(i) * one;
196                         it++;
197                 }
198                 Yn.push_back(buf);
199         }
200         ynsize++;
201 }
202
203
204 // make Yn longer ... 
205 static void make_Yn_longer(int newsize, const cln::float_format_t& prec)
206 {
207
208         cln::cl_N one = cln::cl_float(1, prec);
209
210         Yn[0].resize(newsize);
211         std::vector<cln::cl_N>::iterator it = Yn[0].begin();
212         it += ynlength;
213         for (int i=ynlength+1; i<=newsize; i++) {
214                 *it = *(it-1) + 1 / cln::cl_N(i) * one;
215                 it++;
216         }
217
218         for (int n=1; n<ynsize; n++) {
219                 Yn[n].resize(newsize);
220                 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
221                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
222                 it += ynlength;
223                 itprev += ynlength;
224                 for (int i=ynlength+n+1; i<=newsize+n; i++) {
225                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
226                         it++;
227                         itprev++;
228                 }
229         }
230         
231         ynlength = newsize;
232 }
233
234
235 // calculates Li(2,x) without EuMac
236 static cln::cl_N Li2_series(const cln::cl_N& x)
237 {
238         cln::cl_N res = x;
239         cln::cl_N resbuf;
240         cln::cl_N num = x;
241         cln::cl_I den = 1; // n^2 = 1
242         unsigned i = 3;
243         do {
244                 resbuf = res;
245                 num = num * x;
246                 den = den + i;  // n^2 = 4, 9, 16, ...
247                 i += 2;
248                 res = res + num / den;
249         } while (res != resbuf);
250         return res;
251 }
252
253
254 // calculates Li(2,x) with EuMac
255 static cln::cl_N Li2_series_EuMac(const cln::cl_N& x)
256 {
257         std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
258         cln::cl_N u = -cln::log(1-x);
259         cln::cl_N factor = u;
260         cln::cl_N res = u - u*u/4;
261         cln::cl_N resbuf;
262         unsigned i = 1;
263         do {
264                 resbuf = res;
265                 factor = factor * u*u / (2*i * (2*i+1));
266                 res = res + (*it) * factor;
267                 it++; // should we check it? or rely on initsize? ...
268                 i++;
269         } while (res != resbuf);
270         return res;
271 }
272
273
274 // calculates Li(n,x), n>2 without EuMac
275 static cln::cl_N Lin_series(int n, const cln::cl_N& x)
276 {
277         cln::cl_N factor = x;
278         cln::cl_N res = x;
279         cln::cl_N resbuf;
280         int i=2;
281         do {
282                 resbuf = res;
283                 factor = factor * x;
284                 res = res + factor / cln::expt(cln::cl_I(i),n);
285                 i++;
286         } while (res != resbuf);
287         return res;
288 }
289
290
291 // calculates Li(n,x), n>2 with EuMac
292 static cln::cl_N Lin_series_EuMac(int n, const cln::cl_N& x)
293 {
294         std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
295         cln::cl_N u = -cln::log(1-x);
296         cln::cl_N factor = u;
297         cln::cl_N res = u;
298         cln::cl_N resbuf;
299         unsigned i=2;
300         do {
301                 resbuf = res;
302                 factor = factor * u / i;
303                 res = res + (*it) * factor;
304                 it++; // should we check it? or rely on initsize? ...
305                 i++;
306         } while (res != resbuf);
307         return res;
308 }
309
310
311 // forward declaration needed by function Li_projection and C below
312 static numeric S_num(int n, int p, const numeric& x);
313
314
315 // helper function for classical polylog Li
316 static cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
317 {
318         // treat n=2 as special case
319         if (n == 2) {
320                 // check if precalculated X0 exists
321                 if (xnsize == 0) {
322                         fill_Xn(0);
323                 }
324
325                 if (cln::realpart(x) < 0.5) {
326                         // choose the faster algorithm
327                         // the switching point was empirically determined. the optimal point
328                         // depends on hardware, Digits, ... so an approx value is okay.
329                         // it solves also the problem with precision due to the u=-log(1-x) transformation
330                         if (cln::abs(cln::realpart(x)) < 0.25) {
331                                 
332                                 return Li2_series(x);
333                         } else {
334                                 return Li2_series_EuMac(x);
335                         }
336                 } else {
337                         // choose the faster algorithm
338                         if (cln::abs(cln::realpart(x)) > 0.75) {
339                                 return -Li2_series(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
340                         } else {
341                                 return -Li2_series_EuMac(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
342                         }
343                 }
344         } else {
345                 // check if precalculated Xn exist
346                 if (n > xnsize+1) {
347                         for (int i=xnsize; i<n-1; i++) {
348                                 fill_Xn(i);
349                         }
350                 }
351
352                 if (cln::realpart(x) < 0.5) {
353                         // choose the faster algorithm
354                         // with n>=12 the "normal" summation always wins against EuMac
355                         if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
356                                 return Lin_series(n, x);
357                         } else {
358                                 return Lin_series_EuMac(n, x);
359                         }
360                 } else {
361                         cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
362                         for (int j=0; j<n-1; j++) {
363                                 result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
364                                         * cln::expt(cln::log(x), j) / cln::factorial(j);
365                         }
366                         return result;
367                 }
368         }
369 }
370
371
372 // helper function for classical polylog Li
373 static numeric Li_num(int n, const numeric& x)
374 {
375         if (n == 1) {
376                 // just a log
377                 return -cln::log(1-x.to_cl_N());
378         }
379         if (x.is_zero()) {
380                 return 0;
381         }
382         if (x == 1) {
383                 // [Kol] (2.22)
384                 return cln::zeta(n);
385         }
386         else if (x == -1) {
387                 // [Kol] (2.22)
388                 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
389         }
390         
391         // what is the desired float format?
392         // first guess: default format
393         cln::float_format_t prec = cln::default_float_format;
394         const cln::cl_N value = x.to_cl_N();
395         // second guess: the argument's format
396         if (!x.real().is_rational())
397                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
398         else if (!x.imag().is_rational())
399                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
400         
401         // [Kol] (5.15)
402         if (cln::abs(value) > 1) {
403                 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
404                 // check if argument is complex. if it is real, the new polylog has to be conjugated.
405                 if (cln::zerop(cln::imagpart(value))) {
406                         if (n & 1) {
407                                 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
408                         }
409                         else {
410                                 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
411                         }
412                 }
413                 else {
414                         if (n & 1) {
415                                 result = result + Li_projection(n, cln::recip(value), prec);
416                         }
417                         else {
418                                 result = result - Li_projection(n, cln::recip(value), prec);
419                         }
420                 }
421                 cln::cl_N add;
422                 for (int j=0; j<n-1; j++) {
423                         add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
424                                         * Li_num(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
425                 }
426                 result = result - add;
427                 return result;
428         }
429         else {
430                 return Li_projection(n, value, prec);
431         }
432 }
433
434
435 // helper function for S(n,p,x)
436 static cln::cl_N numeric_nielsen(int n, int step)
437 {
438         if (step) {
439                 cln::cl_N res;
440                 for (int i=1; i<n; i++) {
441                         res = res + numeric_nielsen(i, step-1) / cln::cl_I(i);
442                 }
443                 return res;
444         }
445         else {
446                 return 1;
447         }
448 }
449
450
451 // helper function for S(n,p,x)
452 // [Kol] (7.2)
453 static cln::cl_N C(int n, int p)
454 {
455         cln::cl_N result;
456
457         for (int k=0; k<p; k++) {
458                 for (int j=0; j<=(n+k-1)/2; j++) {
459                         if (k == 0) {
460                                 if (n & 1) {
461                                         if (j & 1) {
462                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
463                                         }
464                                         else {
465                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
466                                         }
467                                 }
468                         }
469                         else {
470                                 if (k & 1) {
471                                         if (j & 1) {
472                                                 result = result + cln::factorial(n+k-1)
473                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
474                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
475                                         }
476                                         else {
477                                                 result = result - cln::factorial(n+k-1)
478                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
479                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
480                                         }
481                                 }
482                                 else {
483                                         if (j & 1) {
484                                                 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()
485                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
486                                         }
487                                         else {
488                                                 result = result + cln::factorial(n+k-1)
489                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
490                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
491                                         }
492                                 }
493                         }
494                 }
495         }
496         int np = n+p;
497         if ((np-1) & 1) {
498                 if (((np)/2+n) & 1) {
499                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
500                 }
501                 else {
502                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
503                 }
504         }
505
506         return result;
507 }
508
509
510 // helper function for S(n,p,x)
511 // [Kol] remark to (9.1)
512 static cln::cl_N a_k(int k)
513 {
514         cln::cl_N result;
515
516         if (k == 0) {
517                 return 1;
518         }
519
520         result = result;
521         for (int m=2; m<=k; m++) {
522                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
523         }
524
525         return -result / k;
526 }
527
528
529 // helper function for S(n,p,x)
530 // [Kol] remark to (9.1)
531 static cln::cl_N b_k(int k)
532 {
533         cln::cl_N result;
534
535         if (k == 0) {
536                 return 1;
537         }
538
539         result = result;
540         for (int m=2; m<=k; m++) {
541                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
542         }
543
544         return result / k;
545 }
546
547
548 // helper function for S(n,p,x)
549 static cln::cl_N S_series(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
550 {
551         if (p==1) {
552                 return Li_projection(n+1, x, prec);
553         }
554         
555         // check if precalculated values are sufficient
556         if (p > ynsize+1) {
557                 for (int i=ynsize; i<p-1; i++) {
558                         fill_Yn(i, prec);
559                 }
560         }
561
562         // should be done otherwise
563         cln::cl_N xf = x * cln::cl_float(1, prec);
564
565         cln::cl_N res;
566         cln::cl_N resbuf;
567         cln::cl_N factor = cln::expt(xf, p);
568         int i = p;
569         do {
570                 resbuf = res;
571                 if (i-p >= ynlength) {
572                         // make Yn longer
573                         make_Yn_longer(ynlength*2, prec);
574                 }
575                 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? ...
576                 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
577                 factor = factor * xf;
578                 i++;
579         } while (res != resbuf);
580         
581         return res;
582 }
583
584
585 // helper function for S(n,p,x)
586 static cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
587 {
588         // [Kol] (5.3)
589         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
590
591                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
592                         * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
593
594                 for (int s=0; s<n; s++) {
595                         cln::cl_N res2;
596                         for (int r=0; r<p; r++) {
597                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
598                                         * S_series(p-r,n-s,1-x,prec) / cln::factorial(r);
599                         }
600                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
601                 }
602
603                 return result;
604         }
605         
606         return S_series(n, p, x, prec);
607 }
608
609
610 // helper function for S(n,p,x)
611 static numeric S_num(int n, int p, const numeric& x)
612 {
613         if (x == 1) {
614                 if (n == 1) {
615                     // [Kol] (2.22) with (2.21)
616                         return cln::zeta(p+1);
617                 }
618
619                 if (p == 1) {
620                     // [Kol] (2.22)
621                         return cln::zeta(n+1);
622                 }
623
624                 // [Kol] (9.1)
625                 cln::cl_N result;
626                 for (int nu=0; nu<n; nu++) {
627                         for (int rho=0; rho<=p; rho++) {
628                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
629                                         * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
630                         }
631                 }
632                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
633
634                 return result;
635         }
636         else if (x == -1) {
637                 // [Kol] (2.22)
638                 if (p == 1) {
639                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
640                 }
641 //              throw std::runtime_error("don't know how to evaluate this function!");
642         }
643
644         // what is the desired float format?
645         // first guess: default format
646         cln::float_format_t prec = cln::default_float_format;
647         const cln::cl_N value = x.to_cl_N();
648         // second guess: the argument's format
649         if (!x.real().is_rational())
650                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
651         else if (!x.imag().is_rational())
652                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
653
654
655         // [Kol] (5.3)
656         if (cln::realpart(value) < -0.5) {
657
658                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
659                         * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
660
661                 for (int s=0; s<n; s++) {
662                         cln::cl_N res2;
663                         for (int r=0; r<p; r++) {
664                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
665                                         * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
666                         }
667                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
668                 }
669
670                 return result;
671                 
672         }
673         // [Kol] (5.12)
674         if (cln::abs(value) > 1) {
675                 
676                 cln::cl_N result;
677
678                 for (int s=0; s<p; s++) {
679                         for (int r=0; r<=s; r++) {
680                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
681                                         / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
682                                         * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
683                         }
684                 }
685                 result = result * cln::expt(cln::cl_I(-1),n);
686
687                 cln::cl_N res2;
688                 for (int r=0; r<n; r++) {
689                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
690                 }
691                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
692
693                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
694
695                 return result;
696         }
697         else {
698                 return S_projection(n, p, value, prec);
699         }
700 }
701
702
703 // helper function for multiple polylogarithm
704 static cln::cl_N numeric_zsum(int n, std::vector<cln::cl_N>& x, std::vector<cln::cl_N>& m)
705 {
706         cln::cl_N res;
707         if (x.empty()) {
708                 return 1;
709         }
710         for (int i=1; i<n; i++) {
711                 std::vector<cln::cl_N>::iterator be;
712                 std::vector<cln::cl_N>::iterator en;
713                 be = x.begin();
714                 be++;
715                 en = x.end();
716                 std::vector<cln::cl_N> xbuf(be, en);
717                 be = m.begin();
718                 be++;
719                 en = m.end();
720                 std::vector<cln::cl_N> mbuf(be, en);
721                 res = res + cln::expt(x[0],i) / cln::expt(i,m[0]) * numeric_zsum(i, xbuf, mbuf);
722         }
723         return res;
724 }
725
726
727 // helper function for harmonic polylogarithm
728 static cln::cl_N numeric_harmonic(int n, std::vector<cln::cl_N>& m)
729 {
730         cln::cl_N res;
731         if (m.empty()) {
732                 return 1;
733         }
734         for (int i=1; i<n; i++) {
735                 std::vector<cln::cl_N>::iterator be;
736                 std::vector<cln::cl_N>::iterator en;
737                 be = m.begin();
738                 be++;
739                 en = m.end();
740                 std::vector<cln::cl_N> mbuf(be, en);
741                 res = res + cln::recip(cln::expt(i,m[0])) * numeric_harmonic(i, mbuf);
742         }
743         return res;
744 }
745
746
747 /////////////////////////////
748 // end of helper functions //
749 /////////////////////////////
750
751
752 // Polylogarithm and multiple polylogarithm
753
754 static ex Li_eval(const ex& x1, const ex& x2)
755 {
756         if (x2.is_zero()) {
757                 return 0;
758         }
759         else {
760                 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational)))
761                         return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
762                 return Li(x1,x2).hold();
763         }
764 }
765
766 static ex Li_evalf(const ex& x1, const ex& x2)
767 {
768         // classical polylogs
769         if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
770                 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
771         }
772         // multiple polylogs
773         else if (is_a<lst>(x1) && is_a<lst>(x2)) {
774                 for (int i=0; i<x1.nops(); i++) {
775                         if (!is_a<numeric>(x1.op(i)))
776                                 return Li(x1,x2).hold();
777                         if (!is_a<numeric>(x2.op(i)))
778                                 return Li(x1,x2).hold();
779                         if (x2.op(i) >= 1)
780                                 return Li(x1,x2).hold();
781                 }
782
783                 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
784                 cln::cl_N x_1 = ex_to<numeric>(x2.op(x2.nops()-1)).to_cl_N();
785                 std::vector<cln::cl_N> x;
786                 std::vector<cln::cl_N> m;
787                 const int nops = ex_to<numeric>(x1.nops()).to_int();
788                 for (int i=nops-2; i>=0; i--) {
789                         m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
790                         x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
791                 }
792
793                 cln::cl_N res;
794                 cln::cl_N resbuf;
795                 for (int i=nops; true; i++) {
796                         resbuf = res;
797                         res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_zsum(i, x, m);
798                         if (cln::zerop(res-resbuf))
799                                 break;
800                 }
801
802                 return numeric(res);
803
804         }
805
806         return Li(x1,x2).hold();
807 }
808
809 static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
810 {
811         epvector seq;
812         seq.push_back(expair(Li(x1,x2), 0));
813         return pseries(rel,seq);
814 }
815
816 REGISTER_FUNCTION(Li, eval_func(Li_eval).evalf_func(Li_evalf).do_not_evalf_params().series_func(Li_series));
817
818
819 // Nielsen's generalized polylogarithm
820
821 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
822 {
823         if (x2 == 1) {
824                 return Li(x1+1,x3);
825         }
826         if (x3.info(info_flags::numeric) && (!x3.info(info_flags::crational)) && 
827                         x1.info(info_flags::posint) && x2.info(info_flags::posint)) {
828                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
829         }
830         return S(x1,x2,x3).hold();
831 }
832
833 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
834 {
835         if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
836                 if ((x3 == -1) && (x2 != 1)) {
837                         // no formula to evaluate this ... sorry
838 //                      return S(x1,x2,x3).hold();
839                 }
840                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
841         }
842         return S(x1,x2,x3).hold();
843 }
844
845 static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options)
846 {
847         epvector seq;
848         seq.push_back(expair(S(x1,x2,x3), 0));
849         return pseries(rel,seq);
850 }
851
852 REGISTER_FUNCTION(S, eval_func(S_eval).evalf_func(S_evalf).do_not_evalf_params().series_func(S_series));
853
854
855 // Harmonic polylogarithm
856
857 static ex H_eval(const ex& x1, const ex& x2)
858 {
859         if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational))) {
860                 return H(x1,x2).evalf();
861         }
862         return H(x1,x2).hold();
863 }
864
865 static ex H_evalf(const ex& x1, const ex& x2)
866 {
867         if (is_a<lst>(x1) && is_a<numeric>(x2)) {
868                 for (int i=0; i<x1.nops(); i++) {
869                         if (!is_a<numeric>(x1.op(i)))
870                                 return H(x1,x2).hold();
871                 }
872                 if (x2 >= 1) {
873                         return H(x1,x2).hold();
874                 }
875
876                 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
877                 cln::cl_N x_1 = ex_to<numeric>(x2).to_cl_N();
878                 std::vector<cln::cl_N> m;
879                 const int nops = ex_to<numeric>(x1.nops()).to_int();
880                 for (int i=nops-2; i>=0; i--) {
881                         m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
882                 }
883
884                 cln::cl_N res;
885                 cln::cl_N resbuf;
886                 for (int i=nops; true; i++) {
887                         resbuf = res;
888                         res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_harmonic(i, m);
889                         if (cln::zerop(res-resbuf))
890                                 break;
891                 }
892
893                 return numeric(res);
894
895         }
896
897         return H(x1,x2).hold();
898 }
899
900 static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
901 {
902         epvector seq;
903         seq.push_back(expair(H(x1,x2), 0));
904         return pseries(rel,seq);
905 }
906
907 REGISTER_FUNCTION(H, eval_func(H_eval).evalf_func(H_evalf).do_not_evalf_params().series_func(H_series));
908
909
910 // Multiple zeta value
911
912 static ex mZeta_eval(const ex& x1)
913 {
914         return mZeta(x1).hold();
915 }
916
917 static ex mZeta_evalf(const ex& x1)
918 {
919         if (is_a<lst>(x1)) {
920                 for (int i=0; i<x1.nops(); i++) {
921                         if (!x1.op(i).info(info_flags::posint))
922                                 return mZeta(x1).hold();
923                 }
924
925                 const int j = x1.nops();
926
927                 std::vector<int> r(j);
928                 for (int i=0; i<j; i++) {
929                         r[j-1-i] = ex_to<numeric>(x1.op(i)).to_int();
930                 }
931
932                 // check for divergence
933                 if (r[0] == 1) {
934                         return mZeta(x1).hold();
935                 }
936                 
937                 // buffer for subsums
938                 std::vector<cln::cl_N> t(j);
939                 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
940
941                 cln::cl_N t0buf;
942                 int q = 0;
943                 do {
944                         t0buf = t[0];
945                         q++;
946                         t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
947                         for (int k=j-2; k>=0; k--) {
948                                 t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
949                         }
950                 } while (t[0] != t0buf);
951
952                 return numeric(t[0]);
953         }
954
955         return mZeta(x1).hold();
956 }
957
958 static ex mZeta_series(const ex& x1, const relational& rel, int order, unsigned options)
959 {
960         epvector seq;
961         seq.push_back(expair(mZeta(x1), 0));
962         return pseries(rel,seq);
963 }
964
965 REGISTER_FUNCTION(mZeta, eval_func(mZeta_eval).evalf_func(mZeta_evalf).do_not_evalf_params().series_func(mZeta_series));
966
967
968 } // namespace GiNaC
969