ginac/inifcns_nstdsums.cpp

   1 /** @file inifcns_nstdsums.cpp
   2  *
   3  *  Implementation of some special functions that have a representation as nested sums.
   4  *
   5  *  The functions are:
   6  *    classical polylogarithm              Li(n,x)
   7  *    multiple polylogarithm               Li(lst(n_1,...,n_k),lst(x_1,...,x_k))
   8  *    nielsen's generalized polylogarithm  S(n,p,x)
   9  *    harmonic polylogarithm               H(n,x) or H(lst(n_1,...,n_k),x)
  10  *    multiple zeta value                  zeta(n) or zeta(lst(n_1,...,n_k))
  11  *
  12  *  Some remarks:
  13  *
  14  *    - All formulae used can be looked up in the following publications:
  15  *      [Kol] Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
  16  *      [Cra] Fast Evaluation of Multiple Zeta Sums, R.E.Crandall, Math.Comp. 67 (1998), pp. 1163-1172.
  17  *      [ReV] Harmonic Polylogarithms, E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
  18  *
  19  *    - The order of parameters and arguments of H, Li and zeta is defined according to their order in the
  20  *      nested sums representation.
  21  *
  22  *    - Except for the multiple polylogarithm all functions can be nummerically evaluated with arguments in
  23  *      the whole complex plane. Multiple polylogarithms evaluate only if each argument x_i is smaller than
  24  *      one. The parameters for every function (n, p or n_i) must be positive integers.
  25  *
  26  *    - The calculation of classical polylogarithms is speed up by using Bernoulli numbers and
  27  *      look-up tables. S uses look-up tables as well. The zeta function applies the algorithm in
  28  *      [Cra] for speed up.
  29  *
  30  *    - The functions have no series expansion as nested sums. To do it, you have to convert these functions
  31  *      into the appropriate objects from the nestedsums library, do the expansion and convert the
  32  *      result back.
  33  *
  34  *    - Numerical testing of this implementation has been performed by doing a comparison of results
  35  *      between this software and the commercial M.......... 4.1. Multiple zeta values have been checked
  36  *      by means of evaluations into simple zeta values. Harmonic polylogarithms have been checked by
  37  *      comparison to S(n,p,x) for corresponding parameter combinations and by continuity checks
  38  *      around |x|=1 along with comparisons to corresponding zeta functions.
  39  *
  40  */
  41
  42 /*
  43  *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
  44  *
  45  *  This program is free software; you can redistribute it and/or modify
  46  *  it under the terms of the GNU General Public License as published by
  47  *  the Free Software Foundation; either version 2 of the License, or
  48  *  (at your option) any later version.
  49  *
  50  *  This program is distributed in the hope that it will be useful,
  51  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  52  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  53  *  GNU General Public License for more details.
  54  *
  55  *  You should have received a copy of the GNU General Public License
  56  *  along with this program; if not, write to the Free Software
  57  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  58  */
  59
  60 #include <stdexcept>
  61 #include <vector>
  62 #include <cln/cln.h>
  63
  64 #include "inifcns.h"
  65
  66 #include "add.h"
  67 #include "constant.h"
  68 #include "lst.h"
  69 #include "mul.h"
  70 #include "numeric.h"
  71 #include "operators.h"
  72 #include "power.h"
  73 #include "pseries.h"
  74 #include "relational.h"
  75 #include "symbol.h"
  76 #include "utils.h"
  77 #include "wildcard.h"
  78
  79
  80 namespace GiNaC {
  81
  82
  83 //////////////////////////////////////////////////////////////////////
  84 //
  85 // Classical polylogarithm  Li
  86 //
  87 // helper functions
  88 //
  89 //////////////////////////////////////////////////////////////////////
  90
  91
  92 // anonymous namespace for helper functions
  93 namespace {
  94
  95
  96 // lookup table for factors built from Bernoulli numbers
  97 // see fill_Xn()
  98 std::vector<std::vector<cln::cl_N> > Xn;
  99 int xnsize = 0;
 100
 101
 102 // This function calculates the X_n. The X_n are needed for speed up of classical polylogarithms.
 103 // With these numbers the polylogs can be calculated as follows:
 104 //   Li_p (x)  =  \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with  u = -log(1-x)
 105 //   X_0(n) = B_n (Bernoulli numbers)
 106 //   X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
 107 // The calculation of Xn depends on X0 and X{n-1}.
 108 // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
 109 // This results in a slightly more complicated algorithm for the X_n.
 110 // The first index in Xn corresponds to the index of the polylog minus 2.
 111 // The second index in Xn corresponds to the index from the actual sum.
 112 void fill_Xn(int n)
 113 {
 114         // rule of thumb. needs to be improved. TODO
 115         const int initsize = Digits * 3 / 2;
 116
 117         if (n>1) {
 118                 // calculate X_2 and higher (corresponding to Li_4 and higher)
 119                 std::vector<cln::cl_N> buf(initsize);
 120                 std::vector<cln::cl_N>::iterator it = buf.begin();
 121                 cln::cl_N result;
 122                 *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
 123                 it++;
 124                 for (int i=2; i<=initsize; i++) {
 125                         if (i&1) {
 126                                 result = 0; // k == 0
 127                         } else {
 128                                 result = Xn[0][i/2-1]; // k == 0
 129                         }
 130                         for (int k=1; k<i-1; k++) {
 131                                 if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
 132                                         result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
 133                                 }
 134                         }
 135                         result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
 136                         result = result + Xn[n-1][i-1] / (i+1); // k == i
 137
 138                         *it = result;
 139                         it++;
 140                 }
 141                 Xn.push_back(buf);
 142         } else if (n==1) {
 143                 // special case to handle the X_0 correct
 144                 std::vector<cln::cl_N> buf(initsize);
 145                 std::vector<cln::cl_N>::iterator it = buf.begin();
 146                 cln::cl_N result;
 147                 *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
 148                 it++;
 149                 *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
 150                 it++;
 151                 for (int i=3; i<=initsize; i++) {
 152                         if (i & 1) {
 153                                 result = -Xn[0][(i-3)/2]/2;
 154                                 *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
 155                                 it++;
 156                         } else {
 157                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
 158                                 for (int k=1; k<i/2; k++) {
 159                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
 160                                 }
 161                                 *it = result;
 162                                 it++;
 163                         }
 164                 }
 165                 Xn.push_back(buf);
 166         } else {
 167                 // calculate X_0
 168                 std::vector<cln::cl_N> buf(initsize/2);
 169                 std::vector<cln::cl_N>::iterator it = buf.begin();
 170                 for (int i=1; i<=initsize/2; i++) {
 171                         *it = bernoulli(i*2).to_cl_N();
 172                         it++;
 173                 }
 174                 Xn.push_back(buf);
 175         }
 176
 177         xnsize++;
 178 }
 179
 180
 181 // calculates Li(2,x) without Xn
 182 cln::cl_N Li2_do_sum(const cln::cl_N& x)
 183 {
 184         cln::cl_N res = x;
 185         cln::cl_N resbuf;
 186         cln::cl_N num = x;
 187         cln::cl_I den = 1; // n^2 = 1
 188         unsigned i = 3;
 189         do {
 190                 resbuf = res;
 191                 num = num * x;
 192                 den = den + i;  // n^2 = 4, 9, 16, ...
 193                 i += 2;
 194                 res = res + num / den;
 195         } while (res != resbuf);
 196         return res;
 197 }
 198
 199
 200 // calculates Li(2,x) with Xn
 201 cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x)
 202 {
 203         std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
 204         cln::cl_N u = -cln::log(1-x);
 205         cln::cl_N factor = u;
 206         cln::cl_N res = u - u*u/4;
 207         cln::cl_N resbuf;
 208         unsigned i = 1;
 209         do {
 210                 resbuf = res;
 211                 factor = factor * u*u / (2*i * (2*i+1));
 212                 res = res + (*it) * factor;
 213                 it++; // should we check it? or rely on initsize? ...
 214                 i++;
 215         } while (res != resbuf);
 216         return res;
 217 }
 218
 219
 220 // calculates Li(n,x), n>2 without Xn
 221 cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
 222 {
 223         cln::cl_N factor = x;
 224         cln::cl_N res = x;
 225         cln::cl_N resbuf;
 226         int i=2;
 227         do {
 228                 resbuf = res;
 229                 factor = factor * x;
 230                 res = res + factor / cln::expt(cln::cl_I(i),n);
 231                 i++;
 232         } while (res != resbuf);
 233         return res;
 234 }
 235
 236
 237 // calculates Li(n,x), n>2 with Xn
 238 cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
 239 {
 240         std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
 241         cln::cl_N u = -cln::log(1-x);
 242         cln::cl_N factor = u;
 243         cln::cl_N res = u;
 244         cln::cl_N resbuf;
 245         unsigned i=2;
 246         do {
 247                 resbuf = res;
 248                 factor = factor * u / i;
 249                 res = res + (*it) * factor;
 250                 it++; // should we check it? or rely on initsize? ...
 251                 i++;
 252         } while (res != resbuf);
 253         return res;
 254 }
 255
 256
 257 // forward declaration needed by function Li_projection and C below
 258 numeric S_num(int n, int p, const numeric& x);
 259
 260
 261 // helper function for classical polylog Li
 262 cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
 263 {
 264         // treat n=2 as special case
 265         if (n == 2) {
 266                 // check if precalculated X0 exists
 267                 if (xnsize == 0) {
 268                         fill_Xn(0);
 269                 }
 270
 271                 if (cln::realpart(x) < 0.5) {
 272                         // choose the faster algorithm
 273                         // the switching point was empirically determined. the optimal point
 274                         // depends on hardware, Digits, ... so an approx value is okay.
 275                         // it solves also the problem with precision due to the u=-log(1-x) transformation
 276                         if (cln::abs(cln::realpart(x)) < 0.25) {
 277
 278                                 return Li2_do_sum(x);
 279                         } else {
 280                                 return Li2_do_sum_Xn(x);
 281                         }
 282                 } else {
 283                         // choose the faster algorithm
 284                         if (cln::abs(cln::realpart(x)) > 0.75) {
 285                                 return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
 286                         } else {
 287                                 return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
 288                         }
 289                 }
 290         } else {
 291                 // check if precalculated Xn exist
 292                 if (n > xnsize+1) {
 293                         for (int i=xnsize; i<n-1; i++) {
 294                                 fill_Xn(i);
 295                         }
 296                 }
 297
 298                 if (cln::realpart(x) < 0.5) {
 299                         // choose the faster algorithm
 300                         // with n>=12 the "normal" summation always wins against the method with Xn
 301                         if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
 302                                 return Lin_do_sum(n, x);
 303                         } else {
 304                                 return Lin_do_sum_Xn(n, x);
 305                         }
 306                 } else {
 307                         cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
 308                         for (int j=0; j<n-1; j++) {
 309                                 result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
 310                                         * cln::expt(cln::log(x), j) / cln::factorial(j);
 311                         }
 312                         return result;
 313                 }
 314         }
 315 }
 316
 317
 318 // helper function for classical polylog Li
 319 numeric Li_num(int n, const numeric& x)
 320 {
 321         if (n == 1) {
 322                 // just a log
 323                 return -cln::log(1-x.to_cl_N());
 324         }
 325         if (x.is_zero()) {
 326                 return 0;
 327         }
 328         if (x == 1) {
 329                 // [Kol] (2.22)
 330                 return cln::zeta(n);
 331         }
 332         else if (x == -1) {
 333                 // [Kol] (2.22)
 334                 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
 335         }
 336
 337         // what is the desired float format?
 338         // first guess: default format
 339         cln::float_format_t prec = cln::default_float_format;
 340         const cln::cl_N value = x.to_cl_N();
 341         // second guess: the argument's format
 342         if (!x.real().is_rational())
 343                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
 344         else if (!x.imag().is_rational())
 345                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 346
 347         // [Kol] (5.15)
 348         if (cln::abs(value) > 1) {
 349                 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
 350                 // check if argument is complex. if it is real, the new polylog has to be conjugated.
 351                 if (cln::zerop(cln::imagpart(value))) {
 352                         if (n & 1) {
 353                                 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
 354                         }
 355                         else {
 356                                 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
 357                         }
 358                 }
 359                 else {
 360                         if (n & 1) {
 361                                 result = result + Li_projection(n, cln::recip(value), prec);
 362                         }
 363                         else {
 364                                 result = result - Li_projection(n, cln::recip(value), prec);
 365                         }
 366                 }
 367                 cln::cl_N add;
 368                 for (int j=0; j<n-1; j++) {
 369                         add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
 370                                         * Li_num(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
 371                 }
 372                 result = result - add;
 373                 return result;
 374         }
 375         else {
 376                 return Li_projection(n, value, prec);
 377         }
 378 }
 379
 380
 381 } // end of anonymous namespace
 382
 383
 384 //////////////////////////////////////////////////////////////////////
 385 //
 386 // Multiple polylogarithm  Li
 387 //
 388 // helper function
 389 //
 390 //////////////////////////////////////////////////////////////////////
 391
 392
 393 // anonymous namespace for helper function
 394 namespace {
 395
 396
 397 cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
 398 {
 399         const int j = s.size();
 400
 401         std::vector<cln::cl_N> t(j);
 402         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
 403
 404         cln::cl_N t0buf;
 405         int q = 0;
 406         do {
 407                 t0buf = t[0];
 408                 q++;
 409                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
 410                 for (int k=j-2; k>=0; k--) {
 411                         t[k] = t[k] + t[k+1] * cln::expt(x[k], q+j-1-k) / cln::expt(cln::cl_I(q+j-1-k), s[k]);
 412                 }
 413         } while (t[0] != t0buf);
 414
 415         return t[0];
 416 }
 417
 418
 419 } // end of anonymous namespace
 420
 421
 422 //////////////////////////////////////////////////////////////////////
 423 //
 424 // Classical polylogarithm and multiple polylogarithm  Li
 425 //
 426 // GiNaC function
 427 //
 428 //////////////////////////////////////////////////////////////////////
 429
 430
 431 static ex Li_eval(const ex& x1, const ex& x2)
 432 {
 433         if (x2.is_zero()) {
 434                 return _ex0;
 435         }
 436         else {
 437                 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational)))
 438                         return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
 439                 if (is_a<lst>(x2)) {
 440                         for (int i=0; i<x2.nops(); i++) {
 441                                 if (!is_a<numeric>(x2.op(i))) {
 442                                         return Li(x1,x2).hold();
 443                                 }
 444                         }
 445                         return Li(x1,x2).evalf();
 446                 }
 447                 return Li(x1,x2).hold();
 448         }
 449 }
 450
 451
 452 static ex Li_evalf(const ex& x1, const ex& x2)
 453 {
 454         // classical polylogs
 455         if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
 456                 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
 457         }
 458         // multiple polylogs
 459         else if (is_a<lst>(x1) && is_a<lst>(x2)) {
 460                 for (int i=0; i<x1.nops(); i++) {
 461                         if (!x1.op(i).info(info_flags::posint)) {
 462                                 return Li(x1,x2).hold();
 463                         }
 464                         if (!is_a<numeric>(x2.op(i))) {
 465                                 return Li(x1,x2).hold();
 466                         }
 467                         if (x2.op(i) >= 1) {
 468                                 return Li(x1,x2).hold();
 469                         }
 470                 }
 471
 472                 std::vector<int> m;
 473                 std::vector<cln::cl_N> x;
 474                 for (int i=0; i<ex_to<numeric>(x1.nops()).to_int(); i++) {
 475                         m.push_back(ex_to<numeric>(x1.op(i)).to_int());
 476                         x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
 477                 }
 478
 479                 return numeric(multipleLi_do_sum(m, x));
 480         }
 481
 482         return Li(x1,x2).hold();
 483 }
 484
 485
 486 static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
 487 {
 488         epvector seq;
 489         seq.push_back(expair(Li(x1,x2), 0));
 490         return pseries(rel,seq);
 491 }
 492
 493
 494 static ex Li_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
 495 {
 496         GINAC_ASSERT(deriv_param < 2);
 497         if (deriv_param == 0) {
 498                 return _ex0;
 499         }
 500         if (x1 > 0) {
 501                 return Li(x1-1, x2) / x2;
 502         } else {
 503                 return 1/(1-x2);
 504         }
 505 }
 506
 507
 508 REGISTER_FUNCTION(Li,
 509                 eval_func(Li_eval).
 510                 evalf_func(Li_evalf).
 511                 do_not_evalf_params().
 512                 series_func(Li_series).
 513                 derivative_func(Li_deriv));
 514
 515
 516 //////////////////////////////////////////////////////////////////////
 517 //
 518 // Nielsen's generalized polylogarithm  S
 519 //
 520 // helper functions
 521 //
 522 //////////////////////////////////////////////////////////////////////
 523
 524
 525 // anonymous namespace for helper functions
 526 namespace {
 527
 528
 529 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
 530 // see fill_Yn()
 531 std::vector<std::vector<cln::cl_N> > Yn;
 532 int ynsize = 0; // number of Yn[]
 533 int ynlength = 100; // initial length of all Yn[i]
 534
 535
 536 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
 537 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
 538 // representing S_{n,p}(x).
 539 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
 540 // equivalent Z-sum.
 541 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
 542 // representing S_{n,p}(x).
 543 // The calculation of Y_n uses the values from Y_{n-1}.
 544 void fill_Yn(int n, const cln::float_format_t& prec)
 545 {
 546         const int initsize = ynlength;
 547         //const int initsize = initsize_Yn;
 548         cln::cl_N one = cln::cl_float(1, prec);
 549
 550         if (n) {
 551                 std::vector<cln::cl_N> buf(initsize);
 552                 std::vector<cln::cl_N>::iterator it = buf.begin();
 553                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
 554                 *it = (*itprev) / cln::cl_N(n+1) * one;
 555                 it++;
 556                 itprev++;
 557                 // sums with an index smaller than the depth are zero and need not to be calculated.
 558                 // calculation starts with depth, which is n+2)
 559                 for (int i=n+2; i<=initsize+n; i++) {
 560                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
 561                         it++;
 562                         itprev++;
 563                 }
 564                 Yn.push_back(buf);
 565         } else {
 566                 std::vector<cln::cl_N> buf(initsize);
 567                 std::vector<cln::cl_N>::iterator it = buf.begin();
 568                 *it = 1 * one;
 569                 it++;
 570                 for (int i=2; i<=initsize; i++) {
 571                         *it = *(it-1) + 1 / cln::cl_N(i) * one;
 572                         it++;
 573                 }
 574                 Yn.push_back(buf);
 575         }
 576         ynsize++;
 577 }
 578
 579
 580 // make Yn longer ...
 581 void make_Yn_longer(int newsize, const cln::float_format_t& prec)
 582 {
 583
 584         cln::cl_N one = cln::cl_float(1, prec);
 585
 586         Yn[0].resize(newsize);
 587         std::vector<cln::cl_N>::iterator it = Yn[0].begin();
 588         it += ynlength;
 589         for (int i=ynlength+1; i<=newsize; i++) {
 590                 *it = *(it-1) + 1 / cln::cl_N(i) * one;
 591                 it++;
 592         }
 593
 594         for (int n=1; n<ynsize; n++) {
 595                 Yn[n].resize(newsize);
 596                 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
 597                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
 598                 it += ynlength;
 599                 itprev += ynlength;
 600                 for (int i=ynlength+n+1; i<=newsize+n; i++) {
 601                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
 602                         it++;
 603                         itprev++;
 604                 }
 605         }
 606
 607         ynlength = newsize;
 608 }
 609
 610
 611 // helper function for S(n,p,x)
 612 // [Kol] (7.2)
 613 cln::cl_N C(int n, int p)
 614 {
 615         cln::cl_N result;
 616
 617         for (int k=0; k<p; k++) {
 618                 for (int j=0; j<=(n+k-1)/2; j++) {
 619                         if (k == 0) {
 620                                 if (n & 1) {
 621                                         if (j & 1) {
 622                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 623                                         }
 624                                         else {
 625                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 626                                         }
 627                                 }
 628                         }
 629                         else {
 630                                 if (k & 1) {
 631                                         if (j & 1) {
 632                                                 result = result + cln::factorial(n+k-1)
 633                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 634                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 635                                         }
 636                                         else {
 637                                                 result = result - cln::factorial(n+k-1)
 638                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 639                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 640                                         }
 641                                 }
 642                                 else {
 643                                         if (j & 1) {
 644                                                 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()
 645                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 646                                         }
 647                                         else {
 648                                                 result = result + cln::factorial(n+k-1)
 649                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 650                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 651                                         }
 652                                 }
 653                         }
 654                 }
 655         }
 656         int np = n+p;
 657         if ((np-1) & 1) {
 658                 if (((np)/2+n) & 1) {
 659                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 660                 }
 661                 else {
 662                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 663                 }
 664         }
 665
 666         return result;
 667 }
 668
 669
 670 // helper function for S(n,p,x)
 671 // [Kol] remark to (9.1)
 672 cln::cl_N a_k(int k)
 673 {
 674         cln::cl_N result;
 675
 676         if (k == 0) {
 677                 return 1;
 678         }
 679
 680         result = result;
 681         for (int m=2; m<=k; m++) {
 682                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
 683         }
 684
 685         return -result / k;
 686 }
 687
 688
 689 // helper function for S(n,p,x)
 690 // [Kol] remark to (9.1)
 691 cln::cl_N b_k(int k)
 692 {
 693         cln::cl_N result;
 694
 695         if (k == 0) {
 696                 return 1;
 697         }
 698
 699         result = result;
 700         for (int m=2; m<=k; m++) {
 701                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
 702         }
 703
 704         return result / k;
 705 }
 706
 707
 708 // helper function for S(n,p,x)
 709 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 710 {
 711         if (p==1) {
 712                 return Li_projection(n+1, x, prec);
 713         }
 714
 715         // check if precalculated values are sufficient
 716         if (p > ynsize+1) {
 717                 for (int i=ynsize; i<p-1; i++) {
 718                         fill_Yn(i, prec);
 719                 }
 720         }
 721
 722         // should be done otherwise
 723         cln::cl_N xf = x * cln::cl_float(1, prec);
 724
 725         cln::cl_N res;
 726         cln::cl_N resbuf;
 727         cln::cl_N factor = cln::expt(xf, p);
 728         int i = p;
 729         do {
 730                 resbuf = res;
 731                 if (i-p >= ynlength) {
 732                         // make Yn longer
 733                         make_Yn_longer(ynlength*2, prec);
 734                 }
 735                 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? ...
 736                 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
 737                 factor = factor * xf;
 738                 i++;
 739         } while (res != resbuf);
 740
 741         return res;
 742 }
 743
 744
 745 // helper function for S(n,p,x)
 746 cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 747 {
 748         // [Kol] (5.3)
 749         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
 750
 751                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
 752                         * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
 753
 754                 for (int s=0; s<n; s++) {
 755                         cln::cl_N res2;
 756                         for (int r=0; r<p; r++) {
 757                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
 758                                         * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
 759                         }
 760                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
 761                 }
 762
 763                 return result;
 764         }
 765
 766         return S_do_sum(n, p, x, prec);
 767 }
 768
 769
 770 // helper function for S(n,p,x)
 771 numeric S_num(int n, int p, const numeric& x)
 772 {
 773         if (x == 1) {
 774                 if (n == 1) {
 775                     // [Kol] (2.22) with (2.21)
 776                         return cln::zeta(p+1);
 777                 }
 778
 779                 if (p == 1) {
 780                     // [Kol] (2.22)
 781                         return cln::zeta(n+1);
 782                 }
 783
 784                 // [Kol] (9.1)
 785                 cln::cl_N result;
 786                 for (int nu=0; nu<n; nu++) {
 787                         for (int rho=0; rho<=p; rho++) {
 788                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
 789                                         * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
 790                         }
 791                 }
 792                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
 793
 794                 return result;
 795         }
 796         else if (x == -1) {
 797                 // [Kol] (2.22)
 798                 if (p == 1) {
 799                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
 800                 }
 801 //              throw std::runtime_error("don't know how to evaluate this function!");
 802         }
 803
 804         // what is the desired float format?
 805         // first guess: default format
 806         cln::float_format_t prec = cln::default_float_format;
 807         const cln::cl_N value = x.to_cl_N();
 808         // second guess: the argument's format
 809         if (!x.real().is_rational())
 810                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
 811         else if (!x.imag().is_rational())
 812                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 813
 814
 815         // [Kol] (5.3)
 816         if (cln::realpart(value) < -0.5) {
 817
 818                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
 819                         * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
 820
 821                 for (int s=0; s<n; s++) {
 822                         cln::cl_N res2;
 823                         for (int r=0; r<p; r++) {
 824                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
 825                                         * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
 826                         }
 827                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
 828                 }
 829
 830                 return result;
 831
 832         }
 833         // [Kol] (5.12)
 834         if (cln::abs(value) > 1) {
 835
 836                 cln::cl_N result;
 837
 838                 for (int s=0; s<p; s++) {
 839                         for (int r=0; r<=s; r++) {
 840                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
 841                                         / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
 842                                         * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
 843                         }
 844                 }
 845                 result = result * cln::expt(cln::cl_I(-1),n);
 846
 847                 cln::cl_N res2;
 848                 for (int r=0; r<n; r++) {
 849                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
 850                 }
 851                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
 852
 853                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
 854
 855                 return result;
 856         }
 857         else {
 858                 return S_projection(n, p, value, prec);
 859         }
 860 }
 861
 862
 863 } // end of anonymous namespace
 864
 865
 866 //////////////////////////////////////////////////////////////////////
 867 //
 868 // Nielsen's generalized polylogarithm  S
 869 //
 870 // GiNaC function
 871 //
 872 //////////////////////////////////////////////////////////////////////
 873
 874
 875 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
 876 {
 877         if (x2 == 1) {
 878                 return Li(x1+1,x3);
 879         }
 880         if (x3.info(info_flags::numeric) && (!x3.info(info_flags::crational)) &&
 881                         x1.info(info_flags::posint) && x2.info(info_flags::posint)) {
 882                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
 883         }
 884         return S(x1,x2,x3).hold();
 885 }
 886
 887
 888 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
 889 {
 890         if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
 891                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
 892         }
 893         return S(x1,x2,x3).hold();
 894 }
 895
 896
 897 static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options)
 898 {
 899         epvector seq;
 900         seq.push_back(expair(S(x1,x2,x3), 0));
 901         return pseries(rel,seq);
 902 }
 903
 904
 905 static ex S_deriv(const ex& x1, const ex& x2, const ex& x3, unsigned deriv_param)
 906 {
 907         GINAC_ASSERT(deriv_param < 3);
 908         if (deriv_param < 2) {
 909                 return _ex0;
 910         }
 911         if (x1 > 0) {
 912                 return S(x1-1, x2, x3) / x3;
 913         } else {
 914                 return S(x1, x2-1, x3) / (1-x3);
 915         }
 916 }
 917
 918
 919 REGISTER_FUNCTION(S,
 920                 eval_func(S_eval).
 921                 evalf_func(S_evalf).
 922                 do_not_evalf_params().
 923                 series_func(S_series).
 924                 derivative_func(S_deriv));
 925
 926
 927 //////////////////////////////////////////////////////////////////////
 928 //
 929 // Harmonic polylogarithm  H
 930 //
 931 // helper function
 932 //
 933 //////////////////////////////////////////////////////////////////////
 934
 935
 936 // anonymous namespace for helper functions
 937 namespace {
 938
 939
 940 // forward declaration
 941 ex convert_from_RV(const lst& parameterlst, const ex& arg);
 942
 943
 944 // multiplies an one-dimensional H with another H
 945 // [ReV] (18)
 946 ex trafo_H_mult(const ex& h1, const ex& h2)
 947 {
 948         ex res;
 949         ex hshort;
 950         lst hlong;
 951         ex h1nops = h1.op(0).nops();
 952         ex h2nops = h2.op(0).nops();
 953         if (h1nops > 1) {
 954                 hshort = h2.op(0).op(0);
 955                 hlong = ex_to<lst>(h1.op(0));
 956         } else {
 957                 hshort = h1.op(0).op(0);
 958                 if (h2nops > 1) {
 959                         hlong = ex_to<lst>(h2.op(0));
 960                 } else {
 961                         hlong = h2.op(0).op(0);
 962                 }
 963         }
 964         for (int i=0; i<=hlong.nops(); i++) {
 965                 lst newparameter;
 966                 int j=0;
 967                 for (; j<i; j++) {
 968                         newparameter.append(hlong[j]);
 969                 }
 970                 newparameter.append(hshort);
 971                 for (; j<hlong.nops(); j++) {
 972                         newparameter.append(hlong[j]);
 973                 }
 974                 res += H(newparameter, h1.op(1)).hold();
 975         }
 976         return res;
 977 }
 978
 979
 980 // applies trafo_H_mult recursively on expressions
 981 struct map_trafo_H_mult : public map_function
 982 {
 983         ex operator()(const ex& e)
 984         {
 985                 if (is_a<add>(e)) {
 986                         return e.map(*this);
 987                 }
 988
 989                 if (is_a<mul>(e)) {
 990
 991                         ex result = 1;
 992                         ex firstH;
 993                         lst Hlst;
 994                         for (int pos=0; pos<e.nops(); pos++) {
 995                                 if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
 996                                         std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
 997                                         if (name == "H") {
 998                                                 for (ex i=0; i<e.op(pos).op(1); i++) {
 999                                                         Hlst.append(e.op(pos).op(0));
1000                                                 }
1001                                                 continue;
1002                                         }
1003                                 } else if (is_a<function>(e.op(pos))) {
1004                                         std::string name = ex_to<function>(e.op(pos)).get_name();
1005                                         if (name == "H") {
1006                                                 if (e.op(pos).op(0).nops() > 1) {
1007                                                         firstH = e.op(pos);
1008                                                 } else {
1009                                                         Hlst.append(e.op(pos));
1010                                                 }
1011                                                 continue;
1012                                         }
1013                                 }
1014                                 result *= e.op(pos);
1015                         }
1016                         if (firstH == 0) {
1017                                 if (Hlst.nops() > 0) {
1018                                         firstH = Hlst[Hlst.nops()-1];
1019                                         Hlst.remove_last();
1020                                 } else {
1021                                         return e;
1022                                 }
1023                         }
1024
1025                         if (Hlst.nops() > 0) {
1026                                 ex buffer = trafo_H_mult(firstH, Hlst.op(0));
1027                                 result *= buffer;
1028                                 for (int i=1; i<Hlst.nops(); i++) {
1029                                         result *= Hlst.op(i);
1030                                 }
1031                                 result = result.expand();
1032                                 map_trafo_H_mult recursion;
1033                                 return recursion(result);
1034                         } else {
1035                                 return e;
1036                         }
1037
1038                 }
1039                 return e;
1040         }
1041 };
1042
1043
1044 // do integration [ReV] (49)
1045 // put parameter 1 in front of existing parameters
1046 ex trafo_H_prepend_one(const ex& e, const ex& arg)
1047 {
1048         ex h;
1049         std::string name;
1050         if (is_a<function>(e)) {
1051                 name = ex_to<function>(e).get_name();
1052         }
1053         if (name == "H") {
1054                 h = e;
1055         } else {
1056                 for (int i=0; i<e.nops(); i++) {
1057                         if (is_a<function>(e.op(i))) {
1058                                 std::string name = ex_to<function>(e.op(i)).get_name();
1059                                 if (name == "H") {
1060                                         h = e.op(i);
1061                                 }
1062                         }
1063                 }
1064         }
1065         if (h != 0) {
1066                 lst newparameter = ex_to<lst>(h.op(0));
1067                 newparameter.prepend(1);
1068                 return e.subs(h == H(newparameter, h.op(1)).hold());
1069         } else {
1070                 return e * H(lst(1),1-arg).hold();
1071         }
1072 }
1073
1074
1075 // do integration [ReV] (55)
1076 // put parameter 0 in front of existing parameters
1077 ex trafo_H_prepend_zero(const ex& e, const ex& arg)
1078 {
1079         ex h;
1080         std::string name;
1081         if (is_a<function>(e)) {
1082                 name = ex_to<function>(e).get_name();
1083         }
1084         if (name == "H") {
1085                 h = e;
1086         } else {
1087                 for (int i=0; i<e.nops(); i++) {
1088                         if (is_a<function>(e.op(i))) {
1089                                 std::string name = ex_to<function>(e.op(i)).get_name();
1090                                 if (name == "H") {
1091                                         h = e.op(i);
1092                                 }
1093                         }
1094                 }
1095         }
1096         if (h != 0) {
1097                 lst newparameter = ex_to<lst>(h.op(0));
1098                 newparameter.prepend(0);
1099                 ex addzeta = convert_from_RV(newparameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1100                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
1101         } else {
1102                 return e * (-H(lst(0),1/arg).hold());
1103         }
1104 }
1105
1106
1107 // do x -> 1-x transformation
1108 struct map_trafo_H_1mx : public map_function
1109 {
1110         ex operator()(const ex& e)
1111         {
1112                 if (is_a<add>(e) || is_a<mul>(e)) {
1113                         return e.map(*this);
1114                 }
1115
1116                 if (is_a<function>(e)) {
1117                         std::string name = ex_to<function>(e).get_name();
1118                         if (name == "H") {
1119
1120                                 lst parameter = ex_to<lst>(e.op(0));
1121                                 ex arg = e.op(1);
1122
1123                                 // if all parameters are either zero or one return the transformed function
1124                                 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1125                                         lst newparameter;
1126                                         for (int i=parameter.nops(); i>0; i--) {
1127                                                 newparameter.append(0);
1128                                         }
1129                                         return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1130                                 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1131                                         lst newparameter;
1132                                         for (int i=parameter.nops(); i>0; i--) {
1133                                                 newparameter.append(1);
1134                                         }
1135                                         return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1136                                 }
1137
1138                                 lst newparameter = parameter;
1139                                 newparameter.remove_first();
1140
1141                                 if (parameter.op(0) == 0) {
1142
1143                                         // leading zero
1144                                         ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1145                                         map_trafo_H_1mx recursion;
1146                                         ex buffer = recursion(H(newparameter, arg).hold());
1147                                         if (is_a<add>(buffer)) {
1148                                                 for (int i=0; i<buffer.nops(); i++) {
1149                                                         res -= trafo_H_prepend_one(buffer.op(i), arg);
1150                                                 }
1151                                         } else {
1152                                                 res -= trafo_H_prepend_one(buffer, arg);
1153                                         }
1154                                         return res;
1155
1156                                 } else {
1157
1158                                         // leading one
1159                                         map_trafo_H_1mx recursion;
1160                                         map_trafo_H_mult unify;
1161                                         ex res;
1162                                         int firstzero = 0;
1163                                         while (parameter.op(firstzero) == 1) {
1164                                                 firstzero++;
1165                                         }
1166                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1167                                                 lst newparameter;
1168                                                 int j=0;
1169                                                 for (; j<=i; j++) {
1170                                                         newparameter.append(parameter[j+1]);
1171                                                 }
1172                                                 newparameter.append(1);
1173                                                 for (; j<parameter.nops()-1; j++) {
1174                                                         newparameter.append(parameter[j+1]);
1175                                                 }
1176                                                 res -= H(newparameter, arg).hold();
1177                                         }
1178                                         return (unify((-H(lst(0), 1-arg).hold() * recursion(H(newparameter, arg).hold())).expand()) +
1179                                                         recursion(res)) / firstzero;
1180
1181                                 }
1182
1183                         }
1184                 }
1185                 return e;
1186         }
1187 };
1188
1189
1190 // do x -> 1/x transformation
1191 struct map_trafo_H_1overx : public map_function
1192 {
1193         ex operator()(const ex& e)
1194         {
1195                 if (is_a<add>(e) || is_a<mul>(e)) {
1196                         return e.map(*this);
1197                 }
1198
1199                 if (is_a<function>(e)) {
1200                         std::string name = ex_to<function>(e).get_name();
1201                         if (name == "H") {
1202
1203                                 lst parameter = ex_to<lst>(e.op(0));
1204                                 ex arg = e.op(1);
1205
1206                                 // if all parameters are either zero or one return the transformed function
1207                                 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1208                                         map_trafo_H_mult unify;
1209                                         return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() - I*Pi, parameter.nops()) /
1210                                                                 factorial(parameter.nops())).expand());
1211                                 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1212                                         return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
1213                                 }
1214
1215                                 lst newparameter = parameter;
1216                                 newparameter.remove_first();
1217
1218                                 if (parameter.op(0) == 0) {
1219
1220                                         // leading zero
1221                                         ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1222                                         map_trafo_H_1overx recursion;
1223                                         ex buffer = recursion(H(newparameter, arg).hold());
1224                                         if (is_a<add>(buffer)) {
1225                                                 for (int i=0; i<buffer.nops(); i++) {
1226                                                         res += trafo_H_prepend_zero(buffer.op(i), arg);
1227                                                 }
1228                                         } else {
1229                                                 res += trafo_H_prepend_zero(buffer, arg);
1230                                         }
1231                                         return res;
1232
1233                                 } else {
1234
1235                                         // leading one
1236                                         map_trafo_H_1overx recursion;
1237                                         map_trafo_H_mult unify;
1238                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
1239                                         int firstzero = 0;
1240                                         while (parameter.op(firstzero) == 1) {
1241                                                 firstzero++;
1242                                         }
1243                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1244                                                 lst newparameter;
1245                                                 int j=0;
1246                                                 for (; j<=i; j++) {
1247                                                         newparameter.append(parameter[j+1]);
1248                                                 }
1249                                                 newparameter.append(1);
1250                                                 for (; j<parameter.nops()-1; j++) {
1251                                                         newparameter.append(parameter[j+1]);
1252                                                 }
1253                                                 res -= H(newparameter, arg).hold();
1254                                         }
1255                                         res = recursion(res).expand() / firstzero;
1256                                         return unify(res);
1257
1258                                 }
1259
1260                         }
1261                 }
1262                 return e;
1263         }
1264 };
1265
1266
1267 // remove trailing zeros from H-parameters
1268 struct map_trafo_H_reduce_trailing_zeros : public map_function
1269 {
1270         ex operator()(const ex& e)
1271         {
1272                 if (is_a<add>(e) || is_a<mul>(e)) {
1273                         return e.map(*this);
1274                 }
1275                 if (is_a<function>(e)) {
1276                         std::string name = ex_to<function>(e).get_name();
1277                         if (name == "H") {
1278                                 lst parameter;
1279                                 if (is_a<lst>(e.op(0))) {
1280                                                 parameter = ex_to<lst>(e.op(0));
1281                                 } else {
1282                                         parameter = lst(e.op(0));
1283                                 }
1284                                 ex arg = e.op(1);
1285                                 if (parameter.op(parameter.nops()-1) == 0) {
1286
1287                                         //
1288                                         if (parameter.nops() == 1) {
1289                                                 return log(arg);
1290                                         }
1291
1292                                         //
1293                                         lst::const_iterator it = parameter.begin();
1294                                         while ((it != parameter.end()) && (*it == 0)) {
1295                                                 it++;
1296                                         }
1297                                         if (it == parameter.end()) {
1298                                                 return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
1299                                         }
1300
1301                                         //
1302                                         parameter.remove_last();
1303                                         int lastentry = parameter.nops();
1304                                         while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
1305                                                 lastentry--;
1306                                         }
1307
1308                                         //
1309                                         ex result = log(arg) * H(parameter,arg).hold();
1310                                         for (ex i=0; i<lastentry; i++) {
1311                                                 parameter[i]++;
1312                                                 result -= (parameter[i]-1) * H(parameter, arg).hold();
1313                                                 parameter[i]--;
1314                                         }
1315
1316                                         if (lastentry < parameter.nops()) {
1317                                                 result = result / (parameter.nops()-lastentry+1);
1318                                                 return result.map(*this);
1319                                         } else {
1320                                                 return result;
1321                                         }
1322                                 }
1323                         }
1324                 }
1325                 return e;
1326         }
1327 };
1328
1329
1330 // recursively call convert_from_RV on expression
1331 struct map_trafo_H_convert : public map_function
1332 {
1333         ex operator()(const ex& e)
1334         {
1335                 if (is_a<add>(e) || is_a<mul>(e) || is_a<power>(e)) {
1336                         return e.map(*this);
1337                 }
1338                 if (is_a<function>(e)) {
1339                         std::string name = ex_to<function>(e).get_name();
1340                         if (name == "H") {
1341                                 lst parameter = ex_to<lst>(e.op(0));
1342                                 ex arg = e.op(1);
1343                                 return convert_from_RV(parameter, arg);
1344                         }
1345                 }
1346                 return e;
1347         }
1348 };
1349
1350
1351 // translate notation from nested sums to Remiddi/Vermaseren
1352 lst convert_to_RV(const lst& o)
1353 {
1354         lst res;
1355         for (lst::const_iterator it = o.begin(); it != o.end(); it++) {
1356                 for (ex i=0; i<(*it)-1; i++) {
1357                         res.append(0);
1358                 }
1359                 res.append(1);
1360         }
1361         return res;
1362 }
1363
1364
1365 // translate notation from Remiddi/Vermaseren to nested sums
1366 ex convert_from_RV(const lst& parameterlst, const ex& arg)
1367 {
1368         lst newparameterlst;
1369
1370         lst::const_iterator it = parameterlst.begin();
1371         int count = 1;
1372         while (it != parameterlst.end()) {
1373                 if (*it == 0) {
1374                         count++;
1375                 } else {
1376                         newparameterlst.append((*it>0) ? count : -count);
1377                         count = 1;
1378                 }
1379                 it++;
1380         }
1381         for (int i=1; i<count; i++) {
1382                 newparameterlst.append(0);
1383         }
1384
1385         map_trafo_H_reduce_trailing_zeros filter;
1386         return filter(H(newparameterlst, arg).hold());
1387 }
1388
1389
1390 // do the actual summation.
1391 cln::cl_N H_do_sum(const std::vector<int>& s, const cln::cl_N& x)
1392 {
1393         const int j = s.size();
1394
1395         std::vector<cln::cl_N> t(j);
1396
1397         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1398         cln::cl_N factor = cln::expt(x, j) * one;
1399         cln::cl_N t0buf;
1400         int q = 0;
1401         do {
1402                 t0buf = t[0];
1403                 q++;
1404                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1405                 for (int k=j-2; k>=1; k--) {
1406                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1407                 }
1408                 t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), s[0]);
1409                 factor = factor * x;
1410         } while (t[0] != t0buf);
1411
1412         return t[0];
1413 }
1414
1415
1416 } // end of anonymous namespace
1417
1418
1419 //////////////////////////////////////////////////////////////////////
1420 //
1421 // Harmonic polylogarithm  H
1422 //
1423 // GiNaC function
1424 //
1425 //////////////////////////////////////////////////////////////////////
1426
1427
1428 static ex H_eval(const ex& x1, const ex& x2)
1429 {
1430         if (x2 == 0) {
1431                 return 0;
1432         }
1433         if (x2 == 1) {
1434                 return zeta(x1);
1435         }
1436         if (x1.nops() == 1) {
1437                 return Li(x1.op(0), x2);
1438         }
1439         if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational))) {
1440                 return H(x1,x2).evalf();
1441         }
1442         return H(x1,x2).hold();
1443 }
1444
1445
1446 static ex H_evalf(const ex& x1, const ex& x2)
1447 {
1448         if (is_a<lst>(x1) && is_a<numeric>(x2)) {
1449                 for (int i=0; i<x1.nops(); i++) {
1450                         if (!x1.op(i).info(info_flags::posint)) {
1451                                 return H(x1,x2).hold();
1452                         }
1453                 }
1454                 if (x1.nops() < 1) {
1455                         return _ex1;
1456                 }
1457                 if (x1.nops() == 1) {
1458                         return Li(x1.op(0), x2).evalf();
1459                 }
1460                 cln::cl_N x = ex_to<numeric>(x2).to_cl_N();
1461                 if (x == 1) {
1462                         return zeta(x1).evalf();
1463                 }
1464
1465                 // choose trafo
1466                 if (cln::abs(x) > 1) {
1467                         symbol xtemp("xtemp");
1468                         map_trafo_H_1overx trafo;
1469                         ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1470                         map_trafo_H_convert converter;
1471                         res = converter(res);
1472                         return res.subs(xtemp==x2).evalf();
1473                 }
1474
1475                 // since the x->1-x transformation produces a lot of terms, it is only
1476                 // efficient for argument near one.
1477                 if (cln::realpart(x) > 0.95) {
1478                         symbol xtemp("xtemp");
1479                         map_trafo_H_1mx trafo;
1480                         ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1481                         map_trafo_H_convert converter;
1482                         res = converter(res);
1483                         return res.subs(xtemp==x2).evalf();
1484                 }
1485
1486                 // no trafo -> do summation
1487                 int count = x1.nops();
1488                 std::vector<int> r(count);
1489                 for (int i=0; i<count; i++) {
1490                         r[i] = ex_to<numeric>(x1.op(i)).to_int();
1491                 }
1492
1493                 return numeric(H_do_sum(r,x));
1494         }
1495
1496         return H(x1,x2).hold();
1497 }
1498
1499
1500 static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
1501 {
1502         epvector seq;
1503         seq.push_back(expair(H(x1,x2), 0));
1504         return pseries(rel,seq);
1505 }
1506
1507
1508 static ex H_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
1509 {
1510         GINAC_ASSERT(deriv_param < 2);
1511         if (deriv_param == 0) {
1512                 return _ex0;
1513         }
1514         if (is_a<lst>(x1)) {
1515                 lst newparameter = ex_to<lst>(x1);
1516                 if (x1.op(0) == 1) {
1517                         newparameter.remove_first();
1518                         return 1/(1-x2) * H(newparameter, x2);
1519                 } else {
1520                         newparameter[0]--;
1521                         return H(newparameter, x2).hold() / x2;
1522                 }
1523         } else {
1524                 if (x1 == 1) {
1525                         return 1/(1-x2);
1526                 } else {
1527                         return H(x1-1, x2).hold() / x2;
1528                 }
1529         }
1530 }
1531
1532
1533 REGISTER_FUNCTION(H,
1534                 eval_func(H_eval).
1535                 evalf_func(H_evalf).
1536                 do_not_evalf_params().
1537                 series_func(H_series).
1538                 derivative_func(H_deriv));
1539
1540
1541 //////////////////////////////////////////////////////////////////////
1542 //
1543 // Multiple zeta values  zeta
1544 //
1545 // helper functions
1546 //
1547 //////////////////////////////////////////////////////////////////////
1548
1549
1550 // anonymous namespace for helper functions
1551 namespace {
1552
1553
1554 // parameters and data for [Cra] algorithm
1555 const cln::cl_N lambda = cln::cl_N("319/320");
1556 int L1;
1557 int L2;
1558 std::vector<std::vector<cln::cl_N> > f_kj;
1559 std::vector<cln::cl_N> crB;
1560 std::vector<std::vector<cln::cl_N> > crG;
1561 std::vector<cln::cl_N> crX;
1562
1563
1564 void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
1565 {
1566         const int size = a.size();
1567         for (int n=0; n<size; n++) {
1568                 c[n] = 0;
1569                 for (int m=0; m<=n; m++) {
1570                         c[n] = c[n] + a[m]*b[n-m];
1571                 }
1572         }
1573 }
1574
1575
1576 // [Cra] section 4
1577 void initcX(const std::vector<int>& s)
1578 {
1579         const int k = s.size();
1580
1581         crX.clear();
1582         crG.clear();
1583         crB.clear();
1584
1585         for (int i=0; i<=L2; i++) {
1586                 crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
1587         }
1588
1589         int Sm = 0;
1590         int Smp1 = 0;
1591         for (int m=0; m<k-1; m++) {
1592                 std::vector<cln::cl_N> crGbuf;
1593                 Sm = Sm + s[m];
1594                 Smp1 = Sm + s[m+1];
1595                 for (int i=0; i<=L2; i++) {
1596                         crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
1597                 }
1598                 crG.push_back(crGbuf);
1599         }
1600
1601         crX = crB;
1602
1603         for (int m=0; m<k-1; m++) {
1604                 std::vector<cln::cl_N> Xbuf;
1605                 for (int i=0; i<=L2; i++) {
1606                         Xbuf.push_back(crX[i] * crG[m][i]);
1607                 }
1608                 halfcyclic_convolute(Xbuf, crB, crX);
1609         }
1610 }
1611
1612
1613 // [Cra] section 4
1614 cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
1615 {
1616         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1617         cln::cl_N factor = cln::expt(lambda, Sqk);
1618         cln::cl_N res = factor / Sqk * crX[0] * one;
1619         cln::cl_N resbuf;
1620         int N = 0;
1621         do {
1622                 resbuf = res;
1623                 factor = factor * lambda;
1624                 N++;
1625                 res = res + crX[N] * factor / (N+Sqk);
1626         } while ((res != resbuf) || cln::zerop(crX[N]));
1627         return res;
1628 }
1629
1630
1631 // [Cra] section 4
1632 void calc_f(int maxr)
1633 {
1634         f_kj.clear();
1635         f_kj.resize(L1);
1636
1637         cln::cl_N t0, t1, t2, t3, t4;
1638         int i, j, k;
1639         std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
1640         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1641
1642         t0 = cln::exp(-lambda);
1643         t2 = 1;
1644         for (k=1; k<=L1; k++) {
1645                 t1 = k * lambda;
1646                 t2 = t0 * t2;
1647                 for (j=1; j<=maxr; j++) {
1648                         t3 = 1;
1649                         t4 = 1;
1650                         for (i=2; i<=j; i++) {
1651                                 t4 = t4 * (j-i+1);
1652                                 t3 = t1 * t3 + t4;
1653                         }
1654                         (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one);
1655                 }
1656                 it++;
1657         }
1658 }
1659
1660
1661 // [Cra] (3.1)
1662 cln::cl_N crandall_Z(const std::vector<int>& s)
1663 {
1664         const int j = s.size();
1665
1666         if (j == 1) {
1667                 cln::cl_N t0;
1668                 cln::cl_N t0buf;
1669                 int q = 0;
1670                 do {
1671                         t0buf = t0;
1672                         q++;
1673                         t0 = t0 + f_kj[q+j-2][s[0]-1];
1674                 } while (t0 != t0buf);
1675
1676                 return t0 / cln::factorial(s[0]-1);
1677         }
1678
1679         std::vector<cln::cl_N> t(j);
1680
1681         cln::cl_N t0buf;
1682         int q = 0;
1683         do {
1684                 t0buf = t[0];
1685                 q++;
1686                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1687                 for (int k=j-2; k>=1; k--) {
1688                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1689                 }
1690                 t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
1691         } while (t[0] != t0buf);
1692
1693         return t[0] / cln::factorial(s[0]-1);
1694 }
1695
1696
1697 // [Cra] (2.4)
1698 cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
1699 {
1700         std::vector<int> r = s;
1701         const int j = r.size();
1702
1703         // decide on maximal size of f_kj for crandall_Z
1704         if (Digits < 50) {
1705                 L1 = 150;
1706         } else {
1707                 L1 = Digits * 3 + j*2;
1708         }
1709
1710         // decide on maximal size of crX for crandall_Y
1711         if (Digits < 38) {
1712                 L2 = 63;
1713         } else if (Digits < 86) {
1714                 L2 = 127;
1715         } else if (Digits < 192) {
1716                 L2 = 255;
1717         } else if (Digits < 394) {
1718                 L2 = 511;
1719         } else if (Digits < 808) {
1720                 L2 = 1023;
1721         } else {
1722                 L2 = 2047;
1723         }
1724
1725         cln::cl_N res;
1726
1727         int maxr = 0;
1728         int S = 0;
1729         for (int i=0; i<j; i++) {
1730                 S += r[i];
1731                 if (r[i] > maxr) {
1732                         maxr = r[i];
1733                 }
1734         }
1735
1736         calc_f(maxr);
1737
1738         const cln::cl_N r0factorial = cln::factorial(r[0]-1);
1739
1740         std::vector<int> rz;
1741         int skp1buf;
1742         int Srun = S;
1743         for (int k=r.size()-1; k>0; k--) {
1744
1745                 rz.insert(rz.begin(), r.back());
1746                 skp1buf = rz.front();
1747                 Srun -= skp1buf;
1748                 r.pop_back();
1749
1750                 initcX(r);
1751
1752                 for (int q=0; q<skp1buf; q++) {
1753
1754                         cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
1755                         cln::cl_N pp2 = crandall_Z(rz);
1756
1757                         rz.front()--;
1758
1759                         if (q & 1) {
1760                                 res = res - pp1 * pp2 / cln::factorial(q);
1761                         } else {
1762                                 res = res + pp1 * pp2 / cln::factorial(q);
1763                         }
1764                 }
1765                 rz.front() = skp1buf;
1766         }
1767         rz.insert(rz.begin(), r.back());
1768
1769         initcX(rz);
1770
1771         res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
1772
1773         return res;
1774 }
1775
1776
1777 cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
1778 {
1779         const int j = r.size();
1780
1781         // buffer for subsums
1782         std::vector<cln::cl_N> t(j);
1783         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1784
1785         cln::cl_N t0buf;
1786         int q = 0;
1787         do {
1788                 t0buf = t[0];
1789                 q++;
1790                 t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
1791                 for (int k=j-2; k>=0; k--) {
1792                         t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
1793                 }
1794         } while (t[0] != t0buf);
1795
1796         return t[0];
1797 }
1798
1799
1800 } // end of anonymous namespace
1801
1802
1803 //////////////////////////////////////////////////////////////////////
1804 //
1805 // Multiple zeta values  zeta
1806 //
1807 // GiNaC function
1808 //
1809 //////////////////////////////////////////////////////////////////////
1810
1811
1812 static ex zeta1_evalf(const ex& x)
1813 {
1814         if (is_exactly_a<lst>(x) && (x.nops()>1)) {
1815
1816                 // multiple zeta value
1817                 const int count = x.nops();
1818                 const lst& xlst = ex_to<lst>(x);
1819                 std::vector<int> r(count);
1820
1821                 // check parameters and convert them
1822                 lst::const_iterator it1 = xlst.begin();
1823                 std::vector<int>::iterator it2 = r.begin();
1824                 do {
1825                         if (!(*it1).info(info_flags::posint)) {
1826                                 return zeta(x).hold();
1827                         }
1828                         *it2 = ex_to<numeric>(*it1).to_int();
1829                         it1++;
1830                         it2++;
1831                 } while (it2 != r.end());
1832
1833                 // check for divergence
1834                 if (r[0] == 1) {
1835                         return zeta(x).hold();
1836                 }
1837
1838                 // decide on summation algorithm
1839                 // this is still a bit clumsy
1840                 int limit = (Digits>17) ? 10 : 6;
1841                 if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) {
1842                         return numeric(zeta_do_sum_Crandall(r));
1843                 } else {
1844                         return numeric(zeta_do_sum_simple(r));
1845                 }
1846         }
1847
1848         // single zeta value
1849         if (is_exactly_a<numeric>(x) && (x != 1)) {
1850                 try {
1851                         return zeta(ex_to<numeric>(x));
1852                 } catch (const dunno &e) { }
1853         }
1854
1855         return zeta(x).hold();
1856 }
1857
1858
1859 static ex zeta1_eval(const ex& x)
1860 {
1861         if (is_exactly_a<lst>(x)) {
1862                 if (x.nops() == 1) {
1863                         return zeta(x.op(0));
1864                 }
1865                 return zeta(x).hold();
1866         }
1867
1868         if (x.info(info_flags::numeric)) {
1869                 const numeric& y = ex_to<numeric>(x);
1870                 // trap integer arguments:
1871                 if (y.is_integer()) {
1872                         if (y.is_zero()) {
1873                                 return _ex_1_2;
1874                         }
1875                         if (y.is_equal(_num1)) {
1876                                 return zeta(x).hold();
1877                         }
1878                         if (y.info(info_flags::posint)) {
1879                                 if (y.info(info_flags::odd)) {
1880                                         return zeta(x).hold();
1881                                 } else {
1882                                         return abs(bernoulli(y)) * pow(Pi, y) * pow(_num2, y-_num1) / factorial(y);
1883                                 }
1884                         } else {
1885                                 if (y.info(info_flags::odd)) {
1886                                         return -bernoulli(_num1-y) / (_num1-y);
1887                                 } else {
1888                                         return _ex0;
1889                                 }
1890                         }
1891                 }
1892                 // zeta(float)
1893                 if (y.info(info_flags::numeric) && !y.info(info_flags::crational))
1894                         return zeta1_evalf(x);
1895         }
1896         return zeta(x).hold();
1897 }
1898
1899
1900 static ex zeta1_deriv(const ex& x, unsigned deriv_param)
1901 {
1902         GINAC_ASSERT(deriv_param==0);
1903
1904         if (is_exactly_a<lst>(x)) {
1905                 return _ex0;
1906         } else {
1907                 return zeta(_ex1, x);
1908         }
1909 }
1910
1911
1912 unsigned zeta1_SERIAL::serial =
1913                         function::register_new(function_options("zeta").
1914                                                 eval_func(zeta1_eval).
1915                                                 evalf_func(zeta1_evalf).
1916                                                 do_not_evalf_params().
1917                                                 derivative_func(zeta1_deriv).
1918                                                 latex_name("\\zeta").
1919                                                 overloaded(2));
1920
1921
1922 //////////////////////////////////////////////////////////////////////
1923 //
1924 // Multiple zeta values  mZeta
1925 //
1926 // The use of mZeta is deprecated! This function will be removed
1927 // from GiNaC source soon. Use zeta instead!!
1928 //
1929 // GiNaC function
1930 //
1931 //////////////////////////////////////////////////////////////////////
1932
1933
1934 static ex mZeta_eval(const ex& x1)
1935 {
1936         return mZeta(x1).hold();
1937 }
1938
1939
1940 static ex mZeta_evalf(const ex& x1)
1941 {
1942         if (is_a<lst>(x1)) {
1943                 for (int i=0; i<x1.nops(); i++) {
1944                         if (!x1.op(i).info(info_flags::posint))
1945                                 return mZeta(x1).hold();
1946                 }
1947
1948                 const int j = x1.nops();
1949
1950                 std::vector<int> r(j);
1951                 for (int i=0; i<j; i++) {
1952                         r[j-1-i] = ex_to<numeric>(x1.op(i)).to_int();
1953                 }
1954
1955                 // check for divergence
1956                 if (r[0] == 1) {
1957                         return mZeta(x1).hold();
1958                 }
1959
1960                 // if only one argument, use cln::zeta
1961                 if (j == 1) {
1962                         return numeric(cln::zeta(r[0]));
1963                 }
1964
1965                 // decide on summation algorithm
1966                 // this is still a bit clumsy
1967                 int limit = (Digits>17) ? 10 : 6;
1968                 if (r[0]<limit || (j>3 && r[1]<limit/2)) {
1969                         return numeric(zeta_do_sum_Crandall(r));
1970                 } else {
1971                         return numeric(zeta_do_sum_simple(r));
1972                 }
1973         } else if (x1.info(info_flags::posint) && (x1 != 1)) {
1974                 return numeric(cln::zeta(ex_to<numeric>(x1).to_int()));
1975         }
1976
1977         return mZeta(x1).hold();
1978 }
1979
1980
1981 static ex mZeta_deriv(const ex& x, unsigned deriv_param)
1982 {
1983         return 0;
1984 }
1985
1986
1987 REGISTER_FUNCTION(mZeta,
1988                 eval_func(mZeta_eval).
1989                 evalf_func(mZeta_evalf).
1990                 do_not_evalf_params().
1991                 derivative_func(mZeta_deriv));
1992
1993
1994 } // namespace GiNaC
1995