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