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