synced to 1.2
[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 except for S(n,p,-1) when p is not unit (no formula yet
17  *      to tackle these points). And of course, there is still room for speed optimizations ;-).
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
57
58 namespace GiNaC {
59
60         
61 //////////////////////
62 // helper functions //
63 //////////////////////
64
65
66 // helper function for classical polylog Li
67 static cln::cl_N Li_series(int n, const cln::cl_N& x, const cln::float_format_t& prec)
68 {
69         // Note: argument must be in the unit circle
70         cln::cl_N aug, acc;
71         cln::cl_N num = cln::complex(cln::cl_float(1, prec), 0);
72         cln::cl_N den = 0;
73         int i = 1;
74         do {
75                 num = num * x;
76                 cln::cl_R ii = i;
77                 den = cln::expt(ii, n);
78                 i++;
79                 aug = num / den;
80                 acc = acc + aug;
81         } while (acc != acc+aug);
82         return acc;
83 }
84
85
86 // helper function for classical polylog Li
87 static cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
88 {
89         return Li_series(n, x, prec);
90 }
91
92
93 // helper function for classical polylog Li
94 static numeric Li_num(int n, const numeric& x)
95 {
96         if (n == 1) {
97                 // just a log
98                 return -cln::log(1-x.to_cl_N());
99         }
100         if (x.is_zero()) {
101                 return 0;
102         }
103         if (x == 1) {
104                 // [Kol] (2.22)
105                 return cln::zeta(n);
106         }
107         else if (x == -1) {
108                 // [Kol] (2.22)
109                 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
110         }
111         
112         // what is the desired float format?
113         // first guess: default format
114         cln::float_format_t prec = cln::default_float_format;
115         const cln::cl_N value = x.to_cl_N();
116         // second guess: the argument's format
117         if (!x.real().is_rational())
118                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
119         else if (!x.imag().is_rational())
120                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
121         
122         // [Kol] (5.15)
123         if (cln::abs(value) > 1) {
124                 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
125                 // check if argument is complex. if it is real, the new polylog has to be conjugated.
126                 if (cln::zerop(cln::imagpart(value))) {
127                         if (n & 1) {
128                                 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
129                         }
130                         else {
131                                 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
132                         }
133                 }
134                 else {
135                         if (n & 1) {
136                                 result = result + Li_projection(n, cln::recip(value), prec);
137                         }
138                         else {
139                                 result = result - Li_projection(n, cln::recip(value), prec);
140                         }
141                 }
142                 cln::cl_N add;
143                 for (int j=0; j<n-1; j++) {
144                         add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
145                                         * Li_num(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
146                 }
147                 result = result - add;
148                 return result;
149         }
150         else {
151                 return Li_projection(n, value, prec);
152         }
153 }
154
155
156 // helper function for S(n,p,x)
157 static cln::cl_N numeric_nielsen(int n, int step)
158 {
159         if (step) {
160                 cln::cl_N res;
161                 for (int i=1; i<n; i++) {
162                         res = res + numeric_nielsen(i, step-1) / cln::cl_I(i);
163                 }
164                 return res;
165         }
166         else {
167                 return 1;
168         }
169 }
170
171
172 // forward declaration needed by function C below
173 static numeric S_num(int n, int p, const numeric& x);
174
175         
176 // helper function for S(n,p,x)
177 // [Kol] (7.2)
178 static cln::cl_N C(int n, int p)
179 {
180         cln::cl_N result;
181
182         for (int k=0; k<p; k++) {
183                 for (int j=0; j<=(n+k-1)/2; j++) {
184                         if (k == 0) {
185                                 if (n & 1) {
186                                         if (j & 1) {
187                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
188                                         }
189                                         else {
190                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
191                                         }
192                                 }
193                         }
194                         else {
195                                 if (k & 1) {
196                                         if (j & 1) {
197                                                 result = result + cln::factorial(n+k-1)
198                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
199                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
200                                         }
201                                         else {
202                                                 result = result - cln::factorial(n+k-1)
203                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
204                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
205                                         }
206                                 }
207                                 else {
208                                         if (j & 1) {
209                                                 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()
210                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
211                                         }
212                                         else {
213                                                 result = result + cln::factorial(n+k-1)
214                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
215                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
216                                         }
217                                 }
218                         }
219                 }
220         }
221         int np = n+p;
222         if ((np-1) & 1) {
223                 if (((np)/2+n) & 1) {
224                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
225                 }
226                 else {
227                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
228                 }
229         }
230
231         return result;
232 }
233
234
235 // helper function for S(n,p,x)
236 // [Kol] remark to (9.1)
237 static cln::cl_N a_k(int k)
238 {
239         cln::cl_N result;
240
241         if (k == 0) {
242                 return 1;
243         }
244
245         result = result;
246         for (int m=2; m<=k; m++) {
247                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
248         }
249
250         return -result / k;
251 }
252
253
254 // helper function for S(n,p,x)
255 // [Kol] remark to (9.1)
256 static cln::cl_N b_k(int k)
257 {
258         cln::cl_N result;
259
260         if (k == 0) {
261                 return 1;
262         }
263
264         result = result;
265         for (int m=2; m<=k; m++) {
266                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
267         }
268
269         return result / k;
270 }
271
272
273 // helper function for S(n,p,x)
274 static cln::cl_N S_series(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
275 {
276         n++;
277         int i = p;
278         p--;
279         cln::cl_N aug, acc;
280         cln::cl_N num = cln::expt(x,p);
281         cln::cl_N converter = cln::complex(cln::cl_float(1, prec), 0);
282         cln::cl_N den = 0;
283         do {
284                 num = num * x;
285                 den = cln::expt(cln::cl_I(i), n);
286                 aug = num / den * numeric_nielsen(i, p);
287                 i++;
288                 acc = acc + aug;
289         } while (acc != acc+aug);
290
291         return acc;
292 }
293
294
295 // helper function for S(n,p,x)
296 static cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
297 {
298         // [Kol] (5.3)
299         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
300
301                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
302                         * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
303
304                 for (int s=0; s<n; s++) {
305                         cln::cl_N res2;
306                         for (int r=0; r<p; r++) {
307                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
308                                         * S_series(p-r,n-s,1-x,prec) / cln::factorial(r);
309                         }
310                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
311                 }
312
313                 return result;
314         }
315         
316         return S_series(n, p, x, prec);
317 }
318
319
320 // helper function for S(n,p,x)
321 static numeric S_num(int n, int p, const numeric& x)
322 {
323         if (x == 1) {
324                 if (n == 1) {
325                     // [Kol] (2.22) with (2.21)
326                         return cln::zeta(p+1);
327                 }
328
329                 if (p == 1) {
330                     // [Kol] (2.22)
331                         return cln::zeta(n+1);
332                 }
333
334                 // [Kol] (9.1)
335                 cln::cl_N result;
336                 for (int nu=0; nu<n; nu++) {
337                         for (int rho=0; rho<=p; rho++) {
338                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
339                                         * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
340                         }
341                 }
342                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
343
344                 return result;
345         }
346         else if (x == -1) {
347                 // [Kol] (2.22)
348                 if (p == 1) {
349                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
350                 }
351                 throw std::runtime_error("don't know how to evaluate this function!");
352         }
353
354         // what is the desired float format?
355         // first guess: default format
356         cln::float_format_t prec = cln::default_float_format;
357         const cln::cl_N value = x.to_cl_N();
358         // second guess: the argument's format
359         if (!x.real().is_rational())
360                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
361         else if (!x.imag().is_rational())
362                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
363
364
365         // [Kol] (5.3)
366         if (cln::realpart(value) < -0.5) {
367
368                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
369                         * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
370
371                 for (int s=0; s<n; s++) {
372                         cln::cl_N res2;
373                         for (int r=0; r<p; r++) {
374                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
375                                         * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
376                         }
377                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
378                 }
379
380                 return result;
381                 
382         }
383         // [Kol] (5.12)
384         else if (cln::abs(value) > 1) {
385                 
386                 cln::cl_N result;
387
388                 for (int s=0; s<p; s++) {
389                         for (int r=0; r<=s; r++) {
390                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
391                                         / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
392                                         * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
393                         }
394                 }
395                 result = result * cln::expt(cln::cl_I(-1),n);
396
397                 cln::cl_N res2;
398                 for (int r=0; r<n; r++) {
399                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
400                 }
401                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
402
403                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
404
405                 return result;
406         }
407         else {
408                 return S_projection(n, p, value, prec);
409         }
410 }
411
412
413 // helper function for multiple polylogarithm
414 static cln::cl_N numeric_zsum(int n, std::vector<cln::cl_N>& x, std::vector<cln::cl_N>& m)
415 {
416         cln::cl_N res;
417         if (x.empty()) {
418                 return 1;
419         }
420         for (int i=1; i<n; i++) {
421                 std::vector<cln::cl_N>::iterator be;
422                 std::vector<cln::cl_N>::iterator en;
423                 be = x.begin();
424                 be++;
425                 en = x.end();
426                 std::vector<cln::cl_N> xbuf(be, en);
427                 be = m.begin();
428                 be++;
429                 en = m.end();
430                 std::vector<cln::cl_N> mbuf(be, en);
431                 res = res + cln::expt(x[0],i) / cln::expt(i,m[0]) * numeric_zsum(i, xbuf, mbuf);
432         }
433         return res;
434 }
435
436
437 // helper function for harmonic polylogarithm
438 static cln::cl_N numeric_harmonic(int n, std::vector<cln::cl_N>& m)
439 {
440         cln::cl_N res;
441         if (m.empty()) {
442                 return 1;
443         }
444         for (int i=1; i<n; i++) {
445                 std::vector<cln::cl_N>::iterator be;
446                 std::vector<cln::cl_N>::iterator en;
447                 be = m.begin();
448                 be++;
449                 en = m.end();
450                 std::vector<cln::cl_N> mbuf(be, en);
451                 res = res + cln::recip(cln::expt(i,m[0])) * numeric_harmonic(i, mbuf);
452         }
453         return res;
454 }
455
456
457 /////////////////////////////
458 // end of helper functions //
459 /////////////////////////////
460
461
462 // Polylogarithm and multiple polylogarithm
463
464 static ex Li_eval(const ex& x1, const ex& x2)
465 {
466         if (x2.is_zero()) {
467                 return 0;
468         }
469         else {
470                 return Li(x1,x2).hold();
471         }
472 }
473
474 static ex Li_evalf(const ex& x1, const ex& x2)
475 {
476         // classical polylogs
477         if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
478                 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
479         }
480         // multiple polylogs
481         else if (is_a<lst>(x1) && is_a<lst>(x2)) {
482                 for (int i=0; i<x1.nops(); i++) {
483                         if (!is_a<numeric>(x1.op(i)))
484                                 return Li(x1,x2).hold();
485                         if (!is_a<numeric>(x2.op(i)))
486                                 return Li(x1,x2).hold();
487                         if (x2 >= 1)
488                                 return Li(x1,x2).hold();
489                 }
490
491                 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
492                 cln::cl_N x_1 = ex_to<numeric>(x2.op(x2.nops()-1)).to_cl_N();
493                 std::vector<cln::cl_N> x;
494                 std::vector<cln::cl_N> m;
495                 const int nops = ex_to<numeric>(x1.nops()).to_int();
496                 for (int i=nops-2; i>=0; i--) {
497                         m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
498                         x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
499                 }
500
501                 cln::cl_N res;
502                 cln::cl_N resbuf;
503                 for (int i=nops; true; i++) {
504                         resbuf = res;
505                         res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_zsum(i, x, m);
506                         if (cln::zerop(res-resbuf))
507                                 break;
508                 }
509
510                 return numeric(res);
511
512         }
513
514         return Li(x1,x2).hold();
515 }
516
517 REGISTER_FUNCTION(Li, eval_func(Li_eval).evalf_func(Li_evalf).do_not_evalf_params());
518
519
520 // Nielsen's generalized polylogarithm
521
522 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
523 {
524         if (x2 == 1) {
525                 return Li(x1+1,x3);
526         }
527         return S(x1,x2,x3).hold();
528 }
529
530 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
531 {
532         if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
533                 if ((x3 == -1) && (x2 != 1)) {
534                         // no formula to evaluate this ... sorry
535                         return S(x1,x2,x3).hold();
536                 }
537                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
538         }
539         return S(x1,x2,x3).hold();
540 }
541
542 REGISTER_FUNCTION(S, eval_func(S_eval).evalf_func(S_evalf).do_not_evalf_params());
543
544
545 // Harmonic polylogarithm
546
547 static ex H_eval(const ex& x1, const ex& x2)
548 {
549         return H(x1,x2).hold();
550 }
551
552 static ex H_evalf(const ex& x1, const ex& x2)
553 {
554         if (is_a<lst>(x1) && is_a<numeric>(x2)) {
555                 for (int i=0; i<x1.nops(); i++) {
556                         if (!is_a<numeric>(x1.op(i)))
557                                 return H(x1,x2).hold();
558                 }
559
560                 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
561                 cln::cl_N x_1 = ex_to<numeric>(x2).to_cl_N();
562                 std::vector<cln::cl_N> m;
563                 const int nops = ex_to<numeric>(x1.nops()).to_int();
564                 for (int i=nops-2; i>=0; i--) {
565                         m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
566                 }
567
568                 cln::cl_N res;
569                 cln::cl_N resbuf;
570                 for (int i=nops; true; i++) {
571                         resbuf = res;
572                         res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_harmonic(i, m);
573                         if (cln::zerop(res-resbuf))
574                                 break;
575                 }
576
577                 return numeric(res);
578
579         }
580
581         return H(x1,x2).hold();
582 }
583
584 REGISTER_FUNCTION(H, eval_func(H_eval).evalf_func(H_evalf).do_not_evalf_params());
585
586
587 // Multiple zeta value
588
589 static ex mZeta_eval(const ex& x1)
590 {
591         return mZeta(x1).hold();
592 }
593
594 static ex mZeta_evalf(const ex& x1)
595 {
596         if (is_a<lst>(x1)) {
597                 for (int i=0; i<x1.nops(); i++) {
598                         if (!is_a<numeric>(x1.op(i)))
599                                 return mZeta(x1).hold();
600                 }
601
602                 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
603                 std::vector<cln::cl_N> m;
604                 const int nops = ex_to<numeric>(x1.nops()).to_int();
605                 for (int i=nops-2; i>=0; i--) {
606                         m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
607                 }
608
609                 cln::float_format_t prec = cln::default_float_format;
610                 cln::cl_N res = cln::complex(cln::cl_float(0, prec), 0);
611                 cln::cl_N resbuf;
612                 for (int i=nops; true; i++) {
613                         // to infinity and beyond ... timewise
614                         resbuf = res;
615                         res = res + cln::recip(cln::expt(i,m_1)) * numeric_harmonic(i, m);
616                         if (cln::zerop(res-resbuf))
617                                 break;
618                 }
619
620                 return numeric(res);
621
622         }
623
624         return mZeta(x1).hold();
625 }
626
627 REGISTER_FUNCTION(mZeta, eval_func(mZeta_eval).evalf_func(mZeta_evalf).do_not_evalf_params());
628
629
630 } // namespace GiNaC
631