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(m_1,...,m_k),lst(x_1,...,x_k))
   8  *                                         G(lst(a_1,...,a_k),y) or G(lst(a_1,...,a_k),lst(s_1,...,s_k),y)
   9  *    Nielsen's generalized polylogarithm  S(n,p,x)
  10  *    harmonic polylogarithm               H(m,x) or H(lst(m_1,...,m_k),x)
  11  *    multiple zeta value                  zeta(m) or zeta(lst(m_1,...,m_k))
  12  *    alternating Euler sum                zeta(m,s) or zeta(lst(m_1,...,m_k),lst(s_1,...,s_k))
  13  *
  14  *  Some remarks:
  15  *
  16  *    - All formulae used can be looked up in the following publications:
  17  *      [Kol] Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
  18  *      [Cra] Fast Evaluation of Multiple Zeta Sums, R.E.Crandall, Math.Comp. 67 (1998), pp. 1163-1172.
  19  *      [ReV] Harmonic Polylogarithms, E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
  20  *      [BBB] Special Values of Multiple Polylogarithms, J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001), pp. 907-941
  21  *      [VSW] Numerical evaluation of multiple polylogarithms, J.Vollinga, S.Weinzierl, hep-ph/0410259
  22  *
  23  *    - The order of parameters and arguments of Li and zeta is defined according to the nested sums
  24  *      representation. The parameters for H are understood as in [ReV]. They can be in expanded --- only
  25  *      0, 1 and -1 --- or in compactified --- a string with zeros in front of 1 or -1 is written as a single
  26  *      number --- notation.
  27  *
  28  *    - All functions can be nummerically evaluated with arguments in the whole complex plane. The parameters
  29  *      for Li, zeta and S must be positive integers. If you want to have an alternating Euler sum, you have
  30  *      to give the signs of the parameters as a second argument s to zeta(m,s) containing 1 and -1.
  31  *
  32  *    - The calculation of classical polylogarithms is speeded up by using Bernoulli numbers and
  33  *      look-up tables. S uses look-up tables as well. The zeta function applies the algorithms in
  34  *      [Cra] and [BBB] for speed up. Multiple polylogarithms use Hoelder convolution [BBB].
  35  *
  36  *    - The functions have no means to do a series expansion into nested sums. To do this, you have to convert
  37  *      these functions into the appropriate objects from the nestedsums library, do the expansion and convert
  38  *      the result back.
  39  *
  40  *    - Numerical testing of this implementation has been performed by doing a comparison of results
  41  *      between this software and the commercial M.......... 4.1. Multiple zeta values have been checked
  42  *      by means of evaluations into simple zeta values. Harmonic polylogarithms have been checked by
  43  *      comparison to S(n,p,x) for corresponding parameter combinations and by continuity checks
  44  *      around |x|=1 along with comparisons to corresponding zeta functions. Multiple polylogarithms were
  45  *      checked against H and zeta and by means of shuffle and quasi-shuffle relations.
  46  *
  47  */
  48
  49 /*
  50  *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
  51  *
  52  *  This program is free software; you can redistribute it and/or modify
  53  *  it under the terms of the GNU General Public License as published by
  54  *  the Free Software Foundation; either version 2 of the License, or
  55  *  (at your option) any later version.
  56  *
  57  *  This program is distributed in the hope that it will be useful,
  58  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  59  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  60  *  GNU General Public License for more details.
  61  *
  62  *  You should have received a copy of the GNU General Public License
  63  *  along with this program; if not, write to the Free Software
  64  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  65  */
  66
  67 #include <sstream>
  68 #include <stdexcept>
  69 #include <vector>
  70 #include <cln/cln.h>
  71
  72 #include "inifcns.h"
  73
  74 #include "add.h"
  75 #include "constant.h"
  76 #include "lst.h"
  77 #include "mul.h"
  78 #include "numeric.h"
  79 #include "operators.h"
  80 #include "power.h"
  81 #include "pseries.h"
  82 #include "relational.h"
  83 #include "symbol.h"
  84 #include "utils.h"
  85 #include "wildcard.h"
  86
  87
  88 namespace GiNaC {
  89
  90
  91 //////////////////////////////////////////////////////////////////////
  92 //
  93 // Classical polylogarithm  Li(n,x)
  94 //
  95 // helper functions
  96 //
  97 //////////////////////////////////////////////////////////////////////
  98
  99
 100 // anonymous namespace for helper functions
 101 namespace {
 102
 103
 104 // lookup table for factors built from Bernoulli numbers
 105 // see fill_Xn()
 106 std::vector<std::vector<cln::cl_N> > Xn;
 107 // initial size of Xn that should suffice for 32bit machines (must be even)
 108 const int xninitsizestep = 26;
 109 int xninitsize = xninitsizestep;
 110 int xnsize = 0;
 111
 112
 113 // This function calculates the X_n. The X_n are needed for speed up of classical polylogarithms.
 114 // With these numbers the polylogs can be calculated as follows:
 115 //   Li_p (x)  =  \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with  u = -log(1-x)
 116 //   X_0(n) = B_n (Bernoulli numbers)
 117 //   X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
 118 // The calculation of Xn depends on X0 and X{n-1}.
 119 // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
 120 // This results in a slightly more complicated algorithm for the X_n.
 121 // The first index in Xn corresponds to the index of the polylog minus 2.
 122 // The second index in Xn corresponds to the index from the actual sum.
 123 void fill_Xn(int n)
 124 {
 125         if (n>1) {
 126                 // calculate X_2 and higher (corresponding to Li_4 and higher)
 127                 std::vector<cln::cl_N> buf(xninitsize);
 128                 std::vector<cln::cl_N>::iterator it = buf.begin();
 129                 cln::cl_N result;
 130                 *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
 131                 it++;
 132                 for (int i=2; i<=xninitsize; i++) {
 133                         if (i&1) {
 134                                 result = 0; // k == 0
 135                         } else {
 136                                 result = Xn[0][i/2-1]; // k == 0
 137                         }
 138                         for (int k=1; k<i-1; k++) {
 139                                 if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
 140                                         result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
 141                                 }
 142                         }
 143                         result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
 144                         result = result + Xn[n-1][i-1] / (i+1); // k == i
 145
 146                         *it = result;
 147                         it++;
 148                 }
 149                 Xn.push_back(buf);
 150         } else if (n==1) {
 151                 // special case to handle the X_0 correct
 152                 std::vector<cln::cl_N> buf(xninitsize);
 153                 std::vector<cln::cl_N>::iterator it = buf.begin();
 154                 cln::cl_N result;
 155                 *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
 156                 it++;
 157                 *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
 158                 it++;
 159                 for (int i=3; i<=xninitsize; i++) {
 160                         if (i & 1) {
 161                                 result = -Xn[0][(i-3)/2]/2;
 162                                 *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
 163                                 it++;
 164                         } else {
 165                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
 166                                 for (int k=1; k<i/2; k++) {
 167                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
 168                                 }
 169                                 *it = result;
 170                                 it++;
 171                         }
 172                 }
 173                 Xn.push_back(buf);
 174         } else {
 175                 // calculate X_0
 176                 std::vector<cln::cl_N> buf(xninitsize/2);
 177                 std::vector<cln::cl_N>::iterator it = buf.begin();
 178                 for (int i=1; i<=xninitsize/2; i++) {
 179                         *it = bernoulli(i*2).to_cl_N();
 180                         it++;
 181                 }
 182                 Xn.push_back(buf);
 183         }
 184
 185         xnsize++;
 186 }
 187
 188 // doubles the number of entries in each Xn[]
 189 void double_Xn()
 190 {
 191         const int pos0 = xninitsize / 2;
 192         // X_0
 193         for (int i=1; i<=xninitsizestep/2; ++i) {
 194                 Xn[0].push_back(bernoulli((i+pos0)*2).to_cl_N());
 195         }
 196         if (Xn.size() > 1) {
 197                 int xend = xninitsize + xninitsizestep;
 198                 cln::cl_N result;
 199                 // X_1
 200                 for (int i=xninitsize+1; i<=xend; ++i) {
 201                         if (i & 1) {
 202                                 result = -Xn[0][(i-3)/2]/2;
 203                                 Xn[1].push_back((cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result);
 204                         } else {
 205                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
 206                                 for (int k=1; k<i/2; k++) {
 207                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
 208                                 }
 209                                 Xn[1].push_back(result);
 210                         }
 211                 }
 212                 // X_n
 213                 for (int n=2; n<Xn.size(); ++n) {
 214                         for (int i=xninitsize+1; i<=xend; ++i) {
 215                                 if (i & 1) {
 216                                         result = 0; // k == 0
 217                                 } else {
 218                                         result = Xn[0][i/2-1]; // k == 0
 219                                 }
 220                                 for (int k=1; k<i-1; ++k) {
 221                                         if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
 222                                                 result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
 223                                         }
 224                                 }
 225                                 result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
 226                                 result = result + Xn[n-1][i-1] / (i+1); // k == i
 227                                 Xn[n].push_back(result);
 228                         }
 229                 }
 230         }
 231         xninitsize += xninitsizestep;
 232 }
 233
 234
 235 // calculates Li(2,x) without Xn
 236 cln::cl_N Li2_do_sum(const cln::cl_N& x)
 237 {
 238         cln::cl_N res = x;
 239         cln::cl_N resbuf;
 240         cln::cl_N num = x * cln::cl_float(1, cln::float_format(Digits));
 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 Xn
 255 cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x)
 256 {
 257         std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
 258         std::vector<cln::cl_N>::const_iterator xend = Xn[0].end();
 259         cln::cl_N u = -cln::log(1-x);
 260         cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
 261         cln::cl_N uu = cln::square(u);
 262         cln::cl_N res = u - uu/4;
 263         cln::cl_N resbuf;
 264         unsigned i = 1;
 265         do {
 266                 resbuf = res;
 267                 factor = factor * uu / (2*i * (2*i+1));
 268                 res = res + (*it) * factor;
 269                 i++;
 270                 if (++it == xend) {
 271                         double_Xn();
 272                         it = Xn[0].begin() + (i-1);
 273                         xend = Xn[0].end();
 274                 }
 275         } while (res != resbuf);
 276         return res;
 277 }
 278
 279
 280 // calculates Li(n,x), n>2 without Xn
 281 cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
 282 {
 283         cln::cl_N factor = x * cln::cl_float(1, cln::float_format(Digits));
 284         cln::cl_N res = x;
 285         cln::cl_N resbuf;
 286         int i=2;
 287         do {
 288                 resbuf = res;
 289                 factor = factor * x;
 290                 res = res + factor / cln::expt(cln::cl_I(i),n);
 291                 i++;
 292         } while (res != resbuf);
 293         return res;
 294 }
 295
 296
 297 // calculates Li(n,x), n>2 with Xn
 298 cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
 299 {
 300         std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
 301         std::vector<cln::cl_N>::const_iterator xend = Xn[n-2].end();
 302         cln::cl_N u = -cln::log(1-x);
 303         cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
 304         cln::cl_N res = u;
 305         cln::cl_N resbuf;
 306         unsigned i=2;
 307         do {
 308                 resbuf = res;
 309                 factor = factor * u / i;
 310                 res = res + (*it) * factor;
 311                 i++;
 312                 if (++it == xend) {
 313                         double_Xn();
 314                         it = Xn[n-2].begin() + (i-2);
 315                         xend = Xn[n-2].end();
 316                 }
 317         } while (res != resbuf);
 318         return res;
 319 }
 320
 321
 322 // forward declaration needed by function Li_projection and C below
 323 const cln::cl_N S_num(int n, int p, const cln::cl_N& x);
 324
 325
 326 // helper function for classical polylog Li
 327 cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
 328 {
 329         // treat n=2 as special case
 330         if (n == 2) {
 331                 // check if precalculated X0 exists
 332                 if (xnsize == 0) {
 333                         fill_Xn(0);
 334                 }
 335
 336                 if (cln::realpart(x) < 0.5) {
 337                         // choose the faster algorithm
 338                         // the switching point was empirically determined. the optimal point
 339                         // depends on hardware, Digits, ... so an approx value is okay.
 340                         // it solves also the problem with precision due to the u=-log(1-x) transformation
 341                         if (cln::abs(cln::realpart(x)) < 0.25) {
 342
 343                                 return Li2_do_sum(x);
 344                         } else {
 345                                 return Li2_do_sum_Xn(x);
 346                         }
 347                 } else {
 348                         // choose the faster algorithm
 349                         if (cln::abs(cln::realpart(x)) > 0.75) {
 350                                 return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
 351                         } else {
 352                                 return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
 353                         }
 354                 }
 355         } else {
 356                 // check if precalculated Xn exist
 357                 if (n > xnsize+1) {
 358                         for (int i=xnsize; i<n-1; i++) {
 359                                 fill_Xn(i);
 360                         }
 361                 }
 362
 363                 if (cln::realpart(x) < 0.5) {
 364                         // choose the faster algorithm
 365                         // with n>=12 the "normal" summation always wins against the method with Xn
 366                         if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
 367                                 return Lin_do_sum(n, x);
 368                         } else {
 369                                 return Lin_do_sum_Xn(n, x);
 370                         }
 371                 } else {
 372                         cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
 373                         for (int j=0; j<n-1; j++) {
 374                                 result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
 375                                                   * cln::expt(cln::log(x), j) / cln::factorial(j);
 376                         }
 377                         return result;
 378                 }
 379         }
 380 }
 381
 382 // helper function for classical polylog Li
 383 const cln::cl_N Lin_numeric(const int n, const cln::cl_N& x)
 384 {
 385         if (n == 1) {
 386                 // just a log
 387                 return -cln::log(1-x);
 388         }
 389         if (zerop(x)) {
 390                 return 0;
 391         }
 392         if (x == 1) {
 393                 // [Kol] (2.22)
 394                 return cln::zeta(n);
 395         }
 396         else if (x == -1) {
 397                 // [Kol] (2.22)
 398                 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
 399         }
 400         if (cln::abs(realpart(x)) < 0.4 && cln::abs(cln::abs(x)-1) < 0.01) {
 401                 cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
 402                 for (int j=0; j<n-1; j++) {
 403                         result = result + (S_num(n-j-1, 1, 1) - S_num(1, n-j-1, 1-x))
 404                                 * cln::expt(cln::log(x), j) / cln::factorial(j);
 405                 }
 406                 return result;
 407         }
 408
 409         // what is the desired float format?
 410         // first guess: default format
 411         cln::float_format_t prec = cln::default_float_format;
 412         const cln::cl_N value = x;
 413         // second guess: the argument's format
 414         if (!instanceof(realpart(x), cln::cl_RA_ring))
 415                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
 416         else if (!instanceof(imagpart(x), cln::cl_RA_ring))
 417                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 418
 419         // [Kol] (5.15)
 420         if (cln::abs(value) > 1) {
 421                 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
 422                 // check if argument is complex. if it is real, the new polylog has to be conjugated.
 423                 if (cln::zerop(cln::imagpart(value))) {
 424                         if (n & 1) {
 425                                 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
 426                         }
 427                         else {
 428                                 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
 429                         }
 430                 }
 431                 else {
 432                         if (n & 1) {
 433                                 result = result + Li_projection(n, cln::recip(value), prec);
 434                         }
 435                         else {
 436                                 result = result - Li_projection(n, cln::recip(value), prec);
 437                         }
 438                 }
 439                 cln::cl_N add;
 440                 for (int j=0; j<n-1; j++) {
 441                         add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
 442                                     * Lin_numeric(n-j,1) * cln::expt(cln::log(-value),j) / cln::factorial(j);
 443                 }
 444                 result = result - add;
 445                 return result;
 446         }
 447         else {
 448                 return Li_projection(n, value, prec);
 449         }
 450 }
 451
 452
 453 } // end of anonymous namespace
 454
 455
 456 //////////////////////////////////////////////////////////////////////
 457 //
 458 // Multiple polylogarithm  Li(n,x)
 459 //
 460 // helper function
 461 //
 462 //////////////////////////////////////////////////////////////////////
 463
 464
 465 // anonymous namespace for helper function
 466 namespace {
 467
 468
 469 // performs the actual series summation for multiple polylogarithms
 470 cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
 471 {
 472         // ensure all x <> 0.
 473         for (std::vector<cln::cl_N>::const_iterator it = x.begin(); it != x.end(); ++it) {
 474                 if ( *it == 0 ) return cln::cl_float(0, cln::float_format(Digits));
 475         }
 476
 477         const int j = s.size();
 478         bool flag_accidental_zero = false;
 479
 480         std::vector<cln::cl_N> t(j);
 481         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
 482
 483         cln::cl_N t0buf;
 484         int q = 0;
 485         do {
 486                 t0buf = t[0];
 487                 q++;
 488                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
 489                 for (int k=j-2; k>=0; k--) {
 490                         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]);
 491                 }
 492                 q++;
 493                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
 494                 for (int k=j-2; k>=0; k--) {
 495                         flag_accidental_zero = cln::zerop(t[k+1]);
 496                         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]);
 497                 }
 498         } while ( (t[0] != t0buf) || cln::zerop(t[0]) || flag_accidental_zero );
 499
 500         return t[0];
 501 }
 502
 503
 504 // converts parameter types and calls multipleLi_do_sum (convenience function for G_numeric)
 505 cln::cl_N mLi_do_summation(const lst& m, const lst& x)
 506 {
 507         std::vector<int> m_int;
 508         std::vector<cln::cl_N> x_cln;
 509         for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
 510                 m_int.push_back(ex_to<numeric>(*itm).to_int());
 511                 x_cln.push_back(ex_to<numeric>(*itx).to_cl_N());
 512         }
 513         return multipleLi_do_sum(m_int, x_cln);
 514 }
 515
 516
 517 // forward declaration for Li_eval()
 518 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf);
 519
 520
 521 // holding dummy-symbols for the G/Li transformations
 522 std::vector<ex> gsyms;
 523
 524
 525 // type used by the transformation functions for G
 526 typedef std::vector<int> Gparameter;
 527
 528
 529 // G_eval1-function for G transformations
 530 ex G_eval1(int a, int scale)
 531 {
 532         if (a != 0) {
 533                 const ex& scs = gsyms[std::abs(scale)];
 534                 const ex& as = gsyms[std::abs(a)];
 535                 if (as != scs) {
 536                         return -log(1 - scs/as);
 537                 } else {
 538                         return -zeta(1);
 539                 }
 540         } else {
 541                 return log(gsyms[std::abs(scale)]);
 542         }
 543 }
 544
 545
 546 // G_eval-function for G transformations
 547 ex G_eval(const Gparameter& a, int scale)
 548 {
 549         // check for properties of G
 550         ex sc = gsyms[std::abs(scale)];
 551         lst newa;
 552         bool all_zero = true;
 553         bool all_ones = true;
 554         int count_ones = 0;
 555         for (Gparameter::const_iterator it = a.begin(); it != a.end(); ++it) {
 556                 if (*it != 0) {
 557                         const ex sym = gsyms[std::abs(*it)];
 558                         newa.append(sym);
 559                         all_zero = false;
 560                         if (sym != sc) {
 561                                 all_ones = false;
 562                         }
 563                         if (all_ones) {
 564                                 ++count_ones;
 565                         }
 566                 } else {
 567                         all_ones = false;
 568                 }
 569         }
 570
 571         // care about divergent G: shuffle to separate divergencies that will be canceled
 572         // later on in the transformation
 573         if (newa.nops() > 1 && newa.op(0) == sc && !all_ones && a.front()!=0) {
 574                 // do shuffle
 575                 Gparameter short_a;
 576                 Gparameter::const_iterator it = a.begin();
 577                 ++it;
 578                 for (; it != a.end(); ++it) {
 579                         short_a.push_back(*it);
 580                 }
 581                 ex result = G_eval1(a.front(), scale) * G_eval(short_a, scale);
 582                 it = short_a.begin();
 583                 for (int i=1; i<count_ones; ++i) {
 584                         ++it;
 585                 }
 586                 for (; it != short_a.end(); ++it) {
 587
 588                         Gparameter newa;
 589                         Gparameter::const_iterator it2 = short_a.begin();
 590                         for (--it2; it2 != it;) {
 591                                 ++it2;
 592                                 newa.push_back(*it2);
 593                         }
 594                         newa.push_back(a[0]);
 595                         ++it2;
 596                         for (; it2 != short_a.end(); ++it2) {
 597                                 newa.push_back(*it2);
 598                         }
 599                         result -= G_eval(newa, scale);
 600                 }
 601                 return result / count_ones;
 602         }
 603
 604         // G({1,...,1};y) -> G({1};y)^k / k!
 605         if (all_ones && a.size() > 1) {
 606                 return pow(G_eval1(a.front(),scale), count_ones) / factorial(count_ones);
 607         }
 608
 609         // G({0,...,0};y) -> log(y)^k / k!
 610         if (all_zero) {
 611                 return pow(log(gsyms[std::abs(scale)]), a.size()) / factorial(a.size());
 612         }
 613
 614         // no special cases anymore -> convert it into Li
 615         lst m;
 616         lst x;
 617         ex argbuf = gsyms[std::abs(scale)];
 618         ex mval = _ex1;
 619         for (Gparameter::const_iterator it=a.begin(); it!=a.end(); ++it) {
 620                 if (*it != 0) {
 621                         const ex& sym = gsyms[std::abs(*it)];
 622                         x.append(argbuf / sym);
 623                         m.append(mval);
 624                         mval = _ex1;
 625                         argbuf = sym;
 626                 } else {
 627                         ++mval;
 628                 }
 629         }
 630         return pow(-1, x.nops()) * Li(m, x);
 631 }
 632
 633
 634 // converts data for G: pending_integrals -> a
 635 Gparameter convert_pending_integrals_G(const Gparameter& pending_integrals)
 636 {
 637         GINAC_ASSERT(pending_integrals.size() != 1);
 638
 639         if (pending_integrals.size() > 0) {
 640                 // get rid of the first element, which would stand for the new upper limit
 641                 Gparameter new_a(pending_integrals.begin()+1, pending_integrals.end());
 642                 return new_a;
 643         } else {
 644                 // just return empty parameter list
 645                 Gparameter new_a;
 646                 return new_a;
 647         }
 648 }
 649
 650
 651 // check the parameters a and scale for G and return information about convergence, depth, etc.
 652 // convergent     : true if G(a,scale) is convergent
 653 // depth          : depth of G(a,scale)
 654 // trailing_zeros : number of trailing zeros of a
 655 // min_it         : iterator of a pointing on the smallest element in a
 656 Gparameter::const_iterator check_parameter_G(const Gparameter& a, int scale,
 657                 bool& convergent, int& depth, int& trailing_zeros, Gparameter::const_iterator& min_it)
 658 {
 659         convergent = true;
 660         depth = 0;
 661         trailing_zeros = 0;
 662         min_it = a.end();
 663         Gparameter::const_iterator lastnonzero = a.end();
 664         for (Gparameter::const_iterator it = a.begin(); it != a.end(); ++it) {
 665                 if (std::abs(*it) > 0) {
 666                         ++depth;
 667                         trailing_zeros = 0;
 668                         lastnonzero = it;
 669                         if (std::abs(*it) < scale) {
 670                                 convergent = false;
 671                                 if ((min_it == a.end()) || (std::abs(*it) < std::abs(*min_it))) {
 672                                         min_it = it;
 673                                 }
 674                         }
 675                 } else {
 676                         ++trailing_zeros;
 677                 }
 678         }
 679         return ++lastnonzero;
 680 }
 681
 682
 683 // add scale to pending_integrals if pending_integrals is empty
 684 Gparameter prepare_pending_integrals(const Gparameter& pending_integrals, int scale)
 685 {
 686         GINAC_ASSERT(pending_integrals.size() != 1);
 687
 688         if (pending_integrals.size() > 0) {
 689                 return pending_integrals;
 690         } else {
 691                 Gparameter new_pending_integrals;
 692                 new_pending_integrals.push_back(scale);
 693                 return new_pending_integrals;
 694         }
 695 }
 696
 697
 698 // handles trailing zeroes for an otherwise convergent integral
 699 ex trailing_zeros_G(const Gparameter& a, int scale)
 700 {
 701         bool convergent;
 702         int depth, trailing_zeros;
 703         Gparameter::const_iterator last, dummyit;
 704         last = check_parameter_G(a, scale, convergent, depth, trailing_zeros, dummyit);
 705
 706         GINAC_ASSERT(convergent);
 707
 708         if ((trailing_zeros > 0) && (depth > 0)) {
 709                 ex result;
 710                 Gparameter new_a(a.begin(), a.end()-1);
 711                 result += G_eval1(0, scale) * trailing_zeros_G(new_a, scale);
 712                 for (Gparameter::const_iterator it = a.begin(); it != last; ++it) {
 713                         Gparameter new_a(a.begin(), it);
 714                         new_a.push_back(0);
 715                         new_a.insert(new_a.end(), it, a.end()-1);
 716                         result -= trailing_zeros_G(new_a, scale);
 717                 }
 718
 719                 return result / trailing_zeros;
 720         } else {
 721                 return G_eval(a, scale);
 722         }
 723 }
 724
 725
 726 // G transformation [VSW] (57),(58)
 727 ex depth_one_trafo_G(const Gparameter& pending_integrals, const Gparameter& a, int scale)
 728 {
 729         // pendint = ( y1, b1, ..., br )
 730         //       a = ( 0, ..., 0, amin )
 731         //   scale = y2
 732         //
 733         // int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(0, ..., 0, sr; y2)
 734         // where sr replaces amin
 735
 736         GINAC_ASSERT(a.back() != 0);
 737         GINAC_ASSERT(a.size() > 0);
 738
 739         ex result;
 740         Gparameter new_pending_integrals = prepare_pending_integrals(pending_integrals, std::abs(a.back()));
 741         const int psize = pending_integrals.size();
 742
 743         // length == 1
 744         // G(sr_{+-}; y2 ) = G(y2_{-+}; sr) - G(0; sr) + ln(-y2_{-+})
 745
 746         if (a.size() == 1) {
 747
 748           // ln(-y2_{-+})
 749           result += log(gsyms[ex_to<numeric>(scale).to_int()]);
 750                 if (a.back() > 0) {
 751                         new_pending_integrals.push_back(-scale);
 752                         result += I*Pi;
 753                 } else {
 754                         new_pending_integrals.push_back(scale);
 755                         result -= I*Pi;
 756                 }
 757                 if (psize) {
 758                         result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front());
 759                 }
 760
 761                 // G(y2_{-+}; sr)
 762                 result += trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front());
 763
 764                 // G(0; sr)
 765                 new_pending_integrals.back() = 0;
 766                 result -= trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals), new_pending_integrals.front());
 767
 768                 return result;
 769         }
 770
 771         // length > 1
 772         // G_m(sr_{+-}; y2) = -zeta_m + int_0^y2 dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
 773         //                            - int_0^sr dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
 774
 775         //term zeta_m
 776         result -= zeta(a.size());
 777         if (psize) {
 778                 result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front());
 779         }
 780
 781         // term int_0^sr dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
 782         //    = int_0^sr dt/t G_{m-1}( t_{+-}; y2 )
 783         Gparameter new_a(a.begin()+1, a.end());
 784         new_pending_integrals.push_back(0);
 785         result -= depth_one_trafo_G(new_pending_integrals, new_a, scale);
 786
 787         // term int_0^y2 dt/t G_{m-1}( (1/y2)_{+-}; 1/t )
 788         //    = int_0^y2 dt/t G_{m-1}( t_{+-}; y2 )
 789         Gparameter new_pending_integrals_2;
 790         new_pending_integrals_2.push_back(scale);
 791         new_pending_integrals_2.push_back(0);
 792         if (psize) {
 793                 result += trailing_zeros_G(convert_pending_integrals_G(pending_integrals), pending_integrals.front())
 794                           * depth_one_trafo_G(new_pending_integrals_2, new_a, scale);
 795         } else {
 796                 result += depth_one_trafo_G(new_pending_integrals_2, new_a, scale);
 797         }
 798
 799         return result;
 800 }
 801
 802
 803 // forward declaration
 804 ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2,
 805              const Gparameter& pendint, const Gparameter& a_old, int scale);
 806
 807
 808 // G transformation [VSW]
 809 ex G_transform(const Gparameter& pendint, const Gparameter& a, int scale)
 810 {
 811         // main recursion routine
 812         //
 813         // pendint = ( y1, b1, ..., br )
 814         //       a = ( a1, ..., amin, ..., aw )
 815         //   scale = y2
 816         //
 817         // int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(a1,...,sr,...,aw,y2)
 818         // where sr replaces amin
 819
 820         // find smallest alpha, determine depth and trailing zeros, and check for convergence
 821         bool convergent;
 822         int depth, trailing_zeros;
 823         Gparameter::const_iterator min_it;
 824         Gparameter::const_iterator firstzero =
 825                 check_parameter_G(a, scale, convergent, depth, trailing_zeros, min_it);
 826         int min_it_pos = min_it - a.begin();
 827
 828         // special case: all a's are zero
 829         if (depth == 0) {
 830                 ex result;
 831
 832                 if (a.size() == 0) {
 833                   result = 1;
 834                 } else {
 835                   result = G_eval(a, scale);
 836                 }
 837                 if (pendint.size() > 0) {
 838                   result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front());
 839                 }
 840                 return result;
 841         }
 842
 843         // handle trailing zeros
 844         if (trailing_zeros > 0) {
 845                 ex result;
 846                 Gparameter new_a(a.begin(), a.end()-1);
 847                 result += G_eval1(0, scale) * G_transform(pendint, new_a, scale);
 848                 for (Gparameter::const_iterator it = a.begin(); it != firstzero; ++it) {
 849                         Gparameter new_a(a.begin(), it);
 850                         new_a.push_back(0);
 851                         new_a.insert(new_a.end(), it, a.end()-1);
 852                         result -= G_transform(pendint, new_a, scale);
 853                 }
 854                 return result / trailing_zeros;
 855         }
 856
 857         // convergence case
 858         if (convergent) {
 859                 if (pendint.size() > 0) {
 860                         return G_eval(convert_pending_integrals_G(pendint), pendint.front()) * G_eval(a, scale);
 861                 } else {
 862                         return G_eval(a, scale);
 863                 }
 864         }
 865
 866         // call basic transformation for depth equal one
 867         if (depth == 1) {
 868                 return depth_one_trafo_G(pendint, a, scale);
 869         }
 870
 871         // do recursion
 872         // int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(a1,...,sr,...,aw,y2)
 873         //  =  int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) G(a1,...,0,...,aw,y2)
 874         //   + int_0^y1 ds1/(s1-b1) ... int dsr/(sr-br) int_0^{sr} ds_{r+1} d/ds_{r+1} G(a1,...,s_{r+1},...,aw,y2)
 875
 876         // smallest element in last place
 877         if (min_it + 1 == a.end()) {
 878                 do { --min_it; } while (*min_it == 0);
 879                 Gparameter empty;
 880                 Gparameter a1(a.begin(),min_it+1);
 881                 Gparameter a2(min_it+1,a.end());
 882
 883                 ex result = G_transform(pendint,a2,scale)*G_transform(empty,a1,scale);
 884
 885                 result -= shuffle_G(empty,a1,a2,pendint,a,scale);
 886                 return result;
 887         }
 888
 889         Gparameter empty;
 890         Gparameter::iterator changeit;
 891
 892         // first term G(a_1,..,0,...,a_w;a_0)
 893         Gparameter new_pendint = prepare_pending_integrals(pendint, a[min_it_pos]);
 894         Gparameter new_a = a;
 895         new_a[min_it_pos] = 0;
 896         ex result = G_transform(empty, new_a, scale);
 897         if (pendint.size() > 0) {
 898                 result *= trailing_zeros_G(convert_pending_integrals_G(pendint), pendint.front());
 899         }
 900
 901         // other terms
 902         changeit = new_a.begin() + min_it_pos;
 903         changeit = new_a.erase(changeit);
 904         if (changeit != new_a.begin()) {
 905                 // smallest in the middle
 906                 new_pendint.push_back(*changeit);
 907                 result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
 908                         * G_transform(empty, new_a, scale);
 909                 int buffer = *changeit;
 910                 *changeit = *min_it;
 911                 result += G_transform(new_pendint, new_a, scale);
 912                 *changeit = buffer;
 913                 new_pendint.pop_back();
 914                 --changeit;
 915                 new_pendint.push_back(*changeit);
 916                 result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
 917                         * G_transform(empty, new_a, scale);
 918                 *changeit = *min_it;
 919                 result -= G_transform(new_pendint, new_a, scale);
 920         } else {
 921                 // smallest at the front
 922                 new_pendint.push_back(scale);
 923                 result += trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
 924                         * G_transform(empty, new_a, scale);
 925                 new_pendint.back() =  *changeit;
 926                 result -= trailing_zeros_G(convert_pending_integrals_G(new_pendint), new_pendint.front())
 927                         * G_transform(empty, new_a, scale);
 928                 *changeit = *min_it;
 929                 result += G_transform(new_pendint, new_a, scale);
 930         }
 931         return result;
 932 }
 933
 934
 935 // shuffles the two parameter list a1 and a2 and calls G_transform for every term except
 936 // for the one that is equal to a_old
 937 ex shuffle_G(const Gparameter & a0, const Gparameter & a1, const Gparameter & a2,
 938              const Gparameter& pendint, const Gparameter& a_old, int scale)
 939 {
 940         if (a1.size()==0 && a2.size()==0) {
 941                 // veto the one configuration we don't want
 942                 if ( a0 == a_old ) return 0;
 943
 944                 return G_transform(pendint,a0,scale);
 945         }
 946
 947         if (a2.size()==0) {
 948                 Gparameter empty;
 949                 Gparameter aa0 = a0;
 950                 aa0.insert(aa0.end(),a1.begin(),a1.end());
 951                 return shuffle_G(aa0,empty,empty,pendint,a_old,scale);
 952         }
 953
 954         if (a1.size()==0) {
 955                 Gparameter empty;
 956                 Gparameter aa0 = a0;
 957                 aa0.insert(aa0.end(),a2.begin(),a2.end());
 958                 return shuffle_G(aa0,empty,empty,pendint,a_old,scale);
 959         }
 960
 961         Gparameter a1_removed(a1.begin()+1,a1.end());
 962         Gparameter a2_removed(a2.begin()+1,a2.end());
 963
 964         Gparameter a01 = a0;
 965         Gparameter a02 = a0;
 966
 967         a01.push_back( a1[0] );
 968         a02.push_back( a2[0] );
 969
 970         return shuffle_G(a01,a1_removed,a2,pendint,a_old,scale)
 971              + shuffle_G(a02,a1,a2_removed,pendint,a_old,scale);
 972 }
 973
 974
 975 // handles the transformations and the numerical evaluation of G
 976 // the parameter x, s and y must only contain numerics
 977 ex G_numeric(const lst& x, const lst& s, const ex& y)
 978 {
 979         // check for convergence and necessary accelerations
 980         bool need_trafo = false;
 981         bool need_hoelder = false;
 982         int depth = 0;
 983         for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
 984                 if (!(*it).is_zero()) {
 985                         ++depth;
 986                         if (abs(*it) - y < -pow(10,-Digits+1)) {
 987                                 need_trafo = true;
 988                         }
 989                         if (abs((abs(*it) - y)/y) < 0.01) {
 990                                 need_hoelder = true;
 991                         }
 992                 }
 993         }
 994         if (x.op(x.nops()-1).is_zero()) {
 995                 need_trafo = true;
 996         }
 997         if (depth == 1 && x.nops() == 2 && !need_trafo) {
 998                 return -Li(x.nops(), y / x.op(x.nops()-1)).evalf();
 999         }
1000
1001         // do acceleration transformation (hoelder convolution [BBB])
1002         if (need_hoelder) {
1003
1004                 ex result;
1005                 const int size = x.nops();
1006                 lst newx;
1007                 for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
1008                         newx.append(*it / y);
1009                 }
1010
1011                 for (int r=0; r<=size; ++r) {
1012                         ex buffer = pow(-1, r);
1013                         ex p = 2;
1014                         bool adjustp;
1015                         do {
1016                                 adjustp = false;
1017                                 for (lst::const_iterator it = newx.begin(); it != newx.end(); ++it) {
1018                                         if (*it == 1/p) {
1019                                                 p += (3-p)/2;
1020                                                 adjustp = true;
1021                                                 continue;
1022                                         }
1023                                 }
1024                         } while (adjustp);
1025                         ex q = p / (p-1);
1026                         lst qlstx;
1027                         lst qlsts;
1028                         for (int j=r; j>=1; --j) {
1029                                 qlstx.append(1-newx.op(j-1));
1030                                 if (newx.op(j-1).info(info_flags::real) && newx.op(j-1) > 1 && newx.op(j-1) <= 2) {
1031                                         qlsts.append( s.op(j-1));
1032                                 } else {
1033                                         qlsts.append( -s.op(j-1));
1034                                 }
1035                         }
1036                         if (qlstx.nops() > 0) {
1037                                 buffer *= G_numeric(qlstx, qlsts, 1/q);
1038                         }
1039                         lst plstx;
1040                         lst plsts;
1041                         for (int j=r+1; j<=size; ++j) {
1042                                 plstx.append(newx.op(j-1));
1043                                 plsts.append(s.op(j-1));
1044                         }
1045                         if (plstx.nops() > 0) {
1046                                 buffer *= G_numeric(plstx, plsts, 1/p);
1047                         }
1048                         result += buffer;
1049                 }
1050                 return result;
1051         }
1052
1053         // convergence transformation
1054         if (need_trafo) {
1055
1056                 // sort (|x|<->position) to determine indices
1057                 std::multimap<ex,int> sortmap;
1058                 int size = 0;
1059                 for (int i=0; i<x.nops(); ++i) {
1060                         if (!x[i].is_zero()) {
1061                                 sortmap.insert(std::pair<ex,int>(abs(x[i]), i));
1062                                 ++size;
1063                         }
1064                 }
1065                 // include upper limit (scale)
1066                 sortmap.insert(std::pair<ex,int>(abs(y), x.nops()));
1067
1068                 // generate missing dummy-symbols
1069                 int i = 1;
1070                 gsyms.clear();
1071                 gsyms.push_back(symbol("GSYMS_ERROR"));
1072                 ex lastentry;
1073                 for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
1074                         if (it != sortmap.begin()) {
1075                                 if (it->second < x.nops()) {
1076                                         if (x[it->second] == lastentry) {
1077                                                 gsyms.push_back(gsyms.back());
1078                                                 continue;
1079                                         }
1080                                 } else {
1081                                         if (y == lastentry) {
1082                                                 gsyms.push_back(gsyms.back());
1083                                                 continue;
1084                                         }
1085                                 }
1086                         }
1087                         std::ostringstream os;
1088                         os << "a" << i;
1089                         gsyms.push_back(symbol(os.str()));
1090                         ++i;
1091                         if (it->second < x.nops()) {
1092                                 lastentry = x[it->second];
1093                         } else {
1094                                 lastentry = y;
1095                         }
1096                 }
1097
1098                 // fill position data according to sorted indices and prepare substitution list
1099                 Gparameter a(x.nops());
1100                 lst subslst;
1101                 int pos = 1;
1102                 int scale;
1103                 for (std::multimap<ex,int>::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
1104                         if (it->second < x.nops()) {
1105                                 if (s[it->second] > 0) {
1106                                         a[it->second] = pos;
1107                                 } else {
1108                                         a[it->second] = -pos;
1109                                 }
1110                                 subslst.append(gsyms[pos] == x[it->second]);
1111                         } else {
1112                                 scale = pos;
1113                                 subslst.append(gsyms[pos] == y);
1114                         }
1115                         ++pos;
1116                 }
1117
1118                 // do transformation
1119                 Gparameter pendint;
1120                 ex result = G_transform(pendint, a, scale);
1121                 // replace dummy symbols with their values
1122                 result = result.eval().expand();
1123                 result = result.subs(subslst).evalf();
1124
1125                 return result;
1126         }
1127
1128         // do summation
1129         lst newx;
1130         lst m;
1131         int mcount = 1;
1132         ex sign = 1;
1133         ex factor = y;
1134         for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
1135                 if ((*it).is_zero()) {
1136                         ++mcount;
1137                 } else {
1138                         newx.append(factor / (*it));
1139                         factor = *it;
1140                         m.append(mcount);
1141                         mcount = 1;
1142                         sign = -sign;
1143                 }
1144         }
1145
1146         return sign * numeric(mLi_do_summation(m, newx));
1147 }
1148
1149
1150 ex mLi_numeric(const lst& m, const lst& x)
1151 {
1152         // let G_numeric do the transformation
1153         lst newx;
1154         lst s;
1155         ex factor = 1;
1156         for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
1157                 for (int i = 1; i < *itm; ++i) {
1158                         newx.append(0);
1159                         s.append(1);
1160                 }
1161                 newx.append(factor / *itx);
1162                 factor /= *itx;
1163                 s.append(1);
1164         }
1165         return pow(-1, m.nops()) * G_numeric(newx, s, _ex1);
1166 }
1167
1168
1169 } // end of anonymous namespace
1170
1171
1172 //////////////////////////////////////////////////////////////////////
1173 //
1174 // Generalized multiple polylogarithm  G(x, y) and G(x, s, y)
1175 //
1176 // GiNaC function
1177 //
1178 //////////////////////////////////////////////////////////////////////
1179
1180
1181 static ex G2_evalf(const ex& x_, const ex& y)
1182 {
1183         if (!y.info(info_flags::positive)) {
1184                 return G(x_, y).hold();
1185         }
1186         lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
1187         if (x.nops() == 0) {
1188                 return _ex1;
1189         }
1190         if (x.op(0) == y) {
1191                 return G(x_, y).hold();
1192         }
1193         lst s;
1194         bool all_zero = true;
1195         for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
1196                 if (!(*it).info(info_flags::numeric)) {
1197                         return G(x_, y).hold();
1198                 }
1199                 if (*it != _ex0) {
1200                         all_zero = false;
1201                 }
1202                 s.append(+1);
1203         }
1204         if (all_zero) {
1205                 return pow(log(y), x.nops()) / factorial(x.nops());
1206         }
1207         return G_numeric(x, s, y);
1208 }
1209
1210
1211 static ex G2_eval(const ex& x_, const ex& y)
1212 {
1213         //TODO eval to MZV or H or S or Lin
1214
1215         if (!y.info(info_flags::positive)) {
1216                 return G(x_, y).hold();
1217         }
1218         lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
1219         if (x.nops() == 0) {
1220                 return _ex1;
1221         }
1222         if (x.op(0) == y) {
1223                 return G(x_, y).hold();
1224         }
1225         lst s;
1226         bool all_zero = true;
1227         bool crational = true;
1228         for (lst::const_iterator it = x.begin(); it != x.end(); ++it) {
1229                 if (!(*it).info(info_flags::numeric)) {
1230                         return G(x_, y).hold();
1231                 }
1232                 if (!(*it).info(info_flags::crational)) {
1233                         crational = false;
1234                 }
1235                 if (*it != _ex0) {
1236                         all_zero = false;
1237                 }
1238                 s.append(+1);
1239         }
1240         if (all_zero) {
1241                 return pow(log(y), x.nops()) / factorial(x.nops());
1242         }
1243         if (!y.info(info_flags::crational)) {
1244                 crational = false;
1245         }
1246         if (crational) {
1247                 return G(x_, y).hold();
1248         }
1249         return G_numeric(x, s, y);
1250 }
1251
1252
1253 unsigned G2_SERIAL::serial = function::register_new(function_options("G", 2).
1254                                 evalf_func(G2_evalf).
1255                                 eval_func(G2_eval).
1256                                 do_not_evalf_params().
1257                                 overloaded(2));
1258 //TODO
1259 //                                derivative_func(G2_deriv).
1260 //                                print_func<print_latex>(G2_print_latex).
1261
1262
1263 static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
1264 {
1265         if (!y.info(info_flags::positive)) {
1266                 return G(x_, s_, y).hold();
1267         }
1268         lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
1269         lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst(s_);
1270         if (x.nops() != s.nops()) {
1271                 return G(x_, s_, y).hold();
1272         }
1273         if (x.nops() == 0) {
1274                 return _ex1;
1275         }
1276         if (x.op(0) == y) {
1277                 return G(x_, s_, y).hold();
1278         }
1279         lst sn;
1280         bool all_zero = true;
1281         for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
1282                 if (!(*itx).info(info_flags::numeric)) {
1283                         return G(x_, y).hold();
1284                 }
1285                 if (!(*its).info(info_flags::real)) {
1286                         return G(x_, y).hold();
1287                 }
1288                 if (*itx != _ex0) {
1289                         all_zero = false;
1290                 }
1291                 if (*its >= 0) {
1292                         sn.append(+1);
1293                 } else {
1294                         sn.append(-1);
1295                 }
1296         }
1297         if (all_zero) {
1298                 return pow(log(y), x.nops()) / factorial(x.nops());
1299         }
1300         return G_numeric(x, sn, y);
1301 }
1302
1303
1304 static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
1305 {
1306         //TODO eval to MZV or H or S or Lin
1307
1308         if (!y.info(info_flags::positive)) {
1309                 return G(x_, s_, y).hold();
1310         }
1311         lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst(x_);
1312         lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst(s_);
1313         if (x.nops() != s.nops()) {
1314                 return G(x_, s_, y).hold();
1315         }
1316         if (x.nops() == 0) {
1317                 return _ex1;
1318         }
1319         if (x.op(0) == y) {
1320                 return G(x_, s_, y).hold();
1321         }
1322         lst sn;
1323         bool all_zero = true;
1324         bool crational = true;
1325         for (lst::const_iterator itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
1326                 if (!(*itx).info(info_flags::numeric)) {
1327                         return G(x_, s_, y).hold();
1328                 }
1329                 if (!(*its).info(info_flags::real)) {
1330                         return G(x_, s_, y).hold();
1331                 }
1332                 if (!(*itx).info(info_flags::crational)) {
1333                         crational = false;
1334                 }
1335                 if (*itx != _ex0) {
1336                         all_zero = false;
1337                 }
1338                 if (*its >= 0) {
1339                         sn.append(+1);
1340                 } else {
1341                         sn.append(-1);
1342                 }
1343         }
1344         if (all_zero) {
1345                 return pow(log(y), x.nops()) / factorial(x.nops());
1346         }
1347         if (!y.info(info_flags::crational)) {
1348                 crational = false;
1349         }
1350         if (crational) {
1351                 return G(x_, s_, y).hold();
1352         }
1353         return G_numeric(x, sn, y);
1354 }
1355
1356
1357 unsigned G3_SERIAL::serial = function::register_new(function_options("G", 3).
1358                                 evalf_func(G3_evalf).
1359                                 eval_func(G3_eval).
1360                                 do_not_evalf_params().
1361                                 overloaded(2));
1362 //TODO
1363 //                                derivative_func(G3_deriv).
1364 //                                print_func<print_latex>(G3_print_latex).
1365
1366
1367 //////////////////////////////////////////////////////////////////////
1368 //
1369 // Classical polylogarithm and multiple polylogarithm  Li(m,x)
1370 //
1371 // GiNaC function
1372 //
1373 //////////////////////////////////////////////////////////////////////
1374
1375
1376 static ex Li_evalf(const ex& m_, const ex& x_)
1377 {
1378         // classical polylogs
1379         if (m_.info(info_flags::posint)) {
1380                 if (x_.info(info_flags::numeric)) {
1381                         int m__ = ex_to<numeric>(m_).to_int();
1382                         const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
1383                         const cln::cl_N result = Lin_numeric(m__, x__);
1384                         return numeric(result);
1385                 } else {
1386                         // try to numerically evaluate second argument
1387                         ex x_val = x_.evalf();
1388                         if (x_val.info(info_flags::numeric)) {
1389                                 int m__ = ex_to<numeric>(m_).to_int();
1390                                 const cln::cl_N x__ = ex_to<numeric>(x_val).to_cl_N();
1391                                 const cln::cl_N result = Lin_numeric(m__, x__);
1392                                 return numeric(result);
1393                         }
1394                 }
1395         }
1396         // multiple polylogs
1397         if (is_a<lst>(m_) && is_a<lst>(x_)) {
1398
1399                 const lst& m = ex_to<lst>(m_);
1400                 const lst& x = ex_to<lst>(x_);
1401                 if (m.nops() != x.nops()) {
1402                         return Li(m_,x_).hold();
1403                 }
1404                 if (x.nops() == 0) {
1405                         return _ex1;
1406                 }
1407                 if ((m.op(0) == _ex1) && (x.op(0) == _ex1)) {
1408                         return Li(m_,x_).hold();
1409                 }
1410
1411                 for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
1412                         if (!(*itm).info(info_flags::posint)) {
1413                                 return Li(m_, x_).hold();
1414                         }
1415                         if (!(*itx).info(info_flags::numeric)) {
1416                                 return Li(m_, x_).hold();
1417                         }
1418                         if (*itx == _ex0) {
1419                                 return _ex0;
1420                         }
1421                 }
1422
1423                 return mLi_numeric(m, x);
1424         }
1425
1426         return Li(m_,x_).hold();
1427 }
1428
1429
1430 static ex Li_eval(const ex& m_, const ex& x_)
1431 {
1432         if (is_a<lst>(m_)) {
1433                 if (is_a<lst>(x_)) {
1434                         // multiple polylogs
1435                         const lst& m = ex_to<lst>(m_);
1436                         const lst& x = ex_to<lst>(x_);
1437                         if (m.nops() != x.nops()) {
1438                                 return Li(m_,x_).hold();
1439                         }
1440                         if (x.nops() == 0) {
1441                                 return _ex1;
1442                         }
1443                         bool is_H = true;
1444                         bool is_zeta = true;
1445                         bool do_evalf = true;
1446                         bool crational = true;
1447                         for (lst::const_iterator itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
1448                                 if (!(*itm).info(info_flags::posint)) {
1449                                         return Li(m_,x_).hold();
1450                                 }
1451                                 if ((*itx != _ex1) && (*itx != _ex_1)) {
1452                                         if (itx != x.begin()) {
1453                                                 is_H = false;
1454                                         }
1455                                         is_zeta = false;
1456                                 }
1457                                 if (*itx == _ex0) {
1458                                         return _ex0;
1459                                 }
1460                                 if (!(*itx).info(info_flags::numeric)) {
1461                                         do_evalf = false;
1462                                 }
1463                                 if (!(*itx).info(info_flags::crational)) {
1464                                         crational = false;
1465                                 }
1466                         }
1467                         if (is_zeta) {
1468                                 return zeta(m_,x_);
1469                         }
1470                         if (is_H) {
1471                                 ex prefactor;
1472                                 lst newm = convert_parameter_Li_to_H(m, x, prefactor);
1473                                 return prefactor * H(newm, x[0]);
1474                         }
1475                         if (do_evalf && !crational) {
1476                                 return mLi_numeric(m,x);
1477                         }
1478                 }
1479                 return Li(m_, x_).hold();
1480         } else if (is_a<lst>(x_)) {
1481                 return Li(m_, x_).hold();
1482         }
1483
1484         // classical polylogs
1485         if (x_ == _ex0) {
1486                 return _ex0;
1487         }
1488         if (x_ == _ex1) {
1489                 return zeta(m_);
1490         }
1491         if (x_ == _ex_1) {
1492                 return (pow(2,1-m_)-1) * zeta(m_);
1493         }
1494         if (m_ == _ex1) {
1495                 return -log(1-x_);
1496         }
1497         if (m_ == _ex2) {
1498                 if (x_.is_equal(I)) {
1499                         return power(Pi,_ex2)/_ex_48 + Catalan*I;
1500                 }
1501                 if (x_.is_equal(-I)) {
1502                         return power(Pi,_ex2)/_ex_48 - Catalan*I;
1503                 }
1504         }
1505         if (m_.info(info_flags::posint) && x_.info(info_flags::numeric) && !x_.info(info_flags::crational)) {
1506                 int m__ = ex_to<numeric>(m_).to_int();
1507                 const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
1508                 const cln::cl_N result = Lin_numeric(m__, x__);
1509                 return numeric(result);
1510         }
1511
1512         return Li(m_, x_).hold();
1513 }
1514
1515
1516 static ex Li_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
1517 {
1518         if (is_a<lst>(m) || is_a<lst>(x)) {
1519                 // multiple polylog
1520                 epvector seq;
1521                 seq.push_back(expair(Li(m, x), 0));
1522                 return pseries(rel, seq);
1523         }
1524
1525         // classical polylog
1526         const ex x_pt = x.subs(rel, subs_options::no_pattern);
1527         if (m.info(info_flags::numeric) && x_pt.info(info_flags::numeric)) {
1528                 // First special case: x==0 (derivatives have poles)
1529                 if (x_pt.is_zero()) {
1530                         const symbol s;
1531                         ex ser;
1532                         // manually construct the primitive expansion
1533                         for (int i=1; i<order; ++i)
1534                                 ser += pow(s,i) / pow(numeric(i), m);
1535                         // substitute the argument's series expansion
1536                         ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
1537                         // maybe that was terminating, so add a proper order term
1538                         epvector nseq;
1539                         nseq.push_back(expair(Order(_ex1), order));
1540                         ser += pseries(rel, nseq);
1541                         // reexpanding it will collapse the series again
1542                         return ser.series(rel, order);
1543                 }
1544                 // TODO special cases: x==1 (branch point) and x real, >=1 (branch cut)
1545                 throw std::runtime_error("Li_series: don't know how to do the series expansion at this point!");
1546         }
1547         // all other cases should be safe, by now:
1548         throw do_taylor();  // caught by function::series()
1549 }
1550
1551
1552 static ex Li_deriv(const ex& m_, const ex& x_, unsigned deriv_param)
1553 {
1554         GINAC_ASSERT(deriv_param < 2);
1555         if (deriv_param == 0) {
1556                 return _ex0;
1557         }
1558         if (m_.nops() > 1) {
1559                 throw std::runtime_error("don't know how to derivate multiple polylogarithm!");
1560         }
1561         ex m;
1562         if (is_a<lst>(m_)) {
1563                 m = m_.op(0);
1564         } else {
1565                 m = m_;
1566         }
1567         ex x;
1568         if (is_a<lst>(x_)) {
1569                 x = x_.op(0);
1570         } else {
1571                 x = x_;
1572         }
1573         if (m > 0) {
1574                 return Li(m-1, x) / x;
1575         } else {
1576                 return 1/(1-x);
1577         }
1578 }
1579
1580
1581 static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c)
1582 {
1583         lst m;
1584         if (is_a<lst>(m_)) {
1585                 m = ex_to<lst>(m_);
1586         } else {
1587                 m = lst(m_);
1588         }
1589         lst x;
1590         if (is_a<lst>(x_)) {
1591                 x = ex_to<lst>(x_);
1592         } else {
1593                 x = lst(x_);
1594         }
1595         c.s << "\\mbox{Li}_{";
1596         lst::const_iterator itm = m.begin();
1597         (*itm).print(c);
1598         itm++;
1599         for (; itm != m.end(); itm++) {
1600                 c.s << ",";
1601                 (*itm).print(c);
1602         }
1603         c.s << "}(";
1604         lst::const_iterator itx = x.begin();
1605         (*itx).print(c);
1606         itx++;
1607         for (; itx != x.end(); itx++) {
1608                 c.s << ",";
1609                 (*itx).print(c);
1610         }
1611         c.s << ")";
1612 }
1613
1614
1615 REGISTER_FUNCTION(Li,
1616                   evalf_func(Li_evalf).
1617                   eval_func(Li_eval).
1618                   series_func(Li_series).
1619                   derivative_func(Li_deriv).
1620                   print_func<print_latex>(Li_print_latex).
1621                   do_not_evalf_params());
1622
1623
1624 //////////////////////////////////////////////////////////////////////
1625 //
1626 // Nielsen's generalized polylogarithm  S(n,p,x)
1627 //
1628 // helper functions
1629 //
1630 //////////////////////////////////////////////////////////////////////
1631
1632
1633 // anonymous namespace for helper functions
1634 namespace {
1635
1636
1637 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
1638 // see fill_Yn()
1639 std::vector<std::vector<cln::cl_N> > Yn;
1640 int ynsize = 0; // number of Yn[]
1641 int ynlength = 100; // initial length of all Yn[i]
1642
1643
1644 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
1645 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
1646 // representing S_{n,p}(x).
1647 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
1648 // equivalent Z-sum.
1649 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
1650 // representing S_{n,p}(x).
1651 // The calculation of Y_n uses the values from Y_{n-1}.
1652 void fill_Yn(int n, const cln::float_format_t& prec)
1653 {
1654         const int initsize = ynlength;
1655         //const int initsize = initsize_Yn;
1656         cln::cl_N one = cln::cl_float(1, prec);
1657
1658         if (n) {
1659                 std::vector<cln::cl_N> buf(initsize);
1660                 std::vector<cln::cl_N>::iterator it = buf.begin();
1661                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
1662                 *it = (*itprev) / cln::cl_N(n+1) * one;
1663                 it++;
1664                 itprev++;
1665                 // sums with an index smaller than the depth are zero and need not to be calculated.
1666                 // calculation starts with depth, which is n+2)
1667                 for (int i=n+2; i<=initsize+n; i++) {
1668                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
1669                         it++;
1670                         itprev++;
1671                 }
1672                 Yn.push_back(buf);
1673         } else {
1674                 std::vector<cln::cl_N> buf(initsize);
1675                 std::vector<cln::cl_N>::iterator it = buf.begin();
1676                 *it = 1 * one;
1677                 it++;
1678                 for (int i=2; i<=initsize; i++) {
1679                         *it = *(it-1) + 1 / cln::cl_N(i) * one;
1680                         it++;
1681                 }
1682                 Yn.push_back(buf);
1683         }
1684         ynsize++;
1685 }
1686
1687
1688 // make Yn longer ...
1689 void make_Yn_longer(int newsize, const cln::float_format_t& prec)
1690 {
1691
1692         cln::cl_N one = cln::cl_float(1, prec);
1693
1694         Yn[0].resize(newsize);
1695         std::vector<cln::cl_N>::iterator it = Yn[0].begin();
1696         it += ynlength;
1697         for (int i=ynlength+1; i<=newsize; i++) {
1698                 *it = *(it-1) + 1 / cln::cl_N(i) * one;
1699                 it++;
1700         }
1701
1702         for (int n=1; n<ynsize; n++) {
1703                 Yn[n].resize(newsize);
1704                 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
1705                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
1706                 it += ynlength;
1707                 itprev += ynlength;
1708                 for (int i=ynlength+n+1; i<=newsize+n; i++) {
1709                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
1710                         it++;
1711                         itprev++;
1712                 }
1713         }
1714
1715         ynlength = newsize;
1716 }
1717
1718
1719 // helper function for S(n,p,x)
1720 // [Kol] (7.2)
1721 cln::cl_N C(int n, int p)
1722 {
1723         cln::cl_N result;
1724
1725         for (int k=0; k<p; k++) {
1726                 for (int j=0; j<=(n+k-1)/2; j++) {
1727                         if (k == 0) {
1728                                 if (n & 1) {
1729                                         if (j & 1) {
1730                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
1731                                         }
1732                                         else {
1733                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
1734                                         }
1735                                 }
1736                         }
1737                         else {
1738                                 if (k & 1) {
1739                                         if (j & 1) {
1740                                                 result = result + cln::factorial(n+k-1)
1741                                                                   * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1742                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1743                                         }
1744                                         else {
1745                                                 result = result - cln::factorial(n+k-1)
1746                                                                   * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1747                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1748                                         }
1749                                 }
1750                                 else {
1751                                         if (j & 1) {
1752                                                 result = result - cln::factorial(n+k-1) * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1753                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1754                                         }
1755                                         else {
1756                                                 result = result + cln::factorial(n+k-1)
1757                                                                   * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1758                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1759                                         }
1760                                 }
1761                         }
1762                 }
1763         }
1764         int np = n+p;
1765         if ((np-1) & 1) {
1766                 if (((np)/2+n) & 1) {
1767                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
1768                 }
1769                 else {
1770                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
1771                 }
1772         }
1773
1774         return result;
1775 }
1776
1777
1778 // helper function for S(n,p,x)
1779 // [Kol] remark to (9.1)
1780 cln::cl_N a_k(int k)
1781 {
1782         cln::cl_N result;
1783
1784         if (k == 0) {
1785                 return 1;
1786         }
1787
1788         result = result;
1789         for (int m=2; m<=k; m++) {
1790                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
1791         }
1792
1793         return -result / k;
1794 }
1795
1796
1797 // helper function for S(n,p,x)
1798 // [Kol] remark to (9.1)
1799 cln::cl_N b_k(int k)
1800 {
1801         cln::cl_N result;
1802
1803         if (k == 0) {
1804                 return 1;
1805         }
1806
1807         result = result;
1808         for (int m=2; m<=k; m++) {
1809                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
1810         }
1811
1812         return result / k;
1813 }
1814
1815
1816 // helper function for S(n,p,x)
1817 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
1818 {
1819         static cln::float_format_t oldprec = cln::default_float_format;
1820
1821         if (p==1) {
1822                 return Li_projection(n+1, x, prec);
1823         }
1824
1825         // precision has changed, we need to clear lookup table Yn
1826         if ( oldprec != prec ) {
1827                 Yn.clear();
1828                 ynsize = 0;
1829                 ynlength = 100;
1830                 oldprec = prec;
1831         }
1832
1833         // check if precalculated values are sufficient
1834         if (p > ynsize+1) {
1835                 for (int i=ynsize; i<p-1; i++) {
1836                         fill_Yn(i, prec);
1837                 }
1838         }
1839
1840         // should be done otherwise
1841         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1842         cln::cl_N xf = x * one;
1843         //cln::cl_N xf = x * cln::cl_float(1, prec);
1844
1845         cln::cl_N res;
1846         cln::cl_N resbuf;
1847         cln::cl_N factor = cln::expt(xf, p);
1848         int i = p;
1849         do {
1850                 resbuf = res;
1851                 if (i-p >= ynlength) {
1852                         // make Yn longer
1853                         make_Yn_longer(ynlength*2, prec);
1854                 }
1855                 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? ...
1856                 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
1857                 factor = factor * xf;
1858                 i++;
1859         } while (res != resbuf);
1860
1861         return res;
1862 }
1863
1864
1865 // helper function for S(n,p,x)
1866 cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
1867 {
1868         // [Kol] (5.3)
1869         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
1870
1871                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
1872                                    * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
1873
1874                 for (int s=0; s<n; s++) {
1875                         cln::cl_N res2;
1876                         for (int r=0; r<p; r++) {
1877                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
1878                                               * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
1879                         }
1880                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
1881                 }
1882
1883                 return result;
1884         }
1885
1886         return S_do_sum(n, p, x, prec);
1887 }
1888
1889
1890 // helper function for S(n,p,x)
1891 const cln::cl_N S_num(int n, int p, const cln::cl_N& x)
1892 {
1893         if (x == 1) {
1894                 if (n == 1) {
1895                     // [Kol] (2.22) with (2.21)
1896                         return cln::zeta(p+1);
1897                 }
1898
1899                 if (p == 1) {
1900                     // [Kol] (2.22)
1901                         return cln::zeta(n+1);
1902                 }
1903
1904                 // [Kol] (9.1)
1905                 cln::cl_N result;
1906                 for (int nu=0; nu<n; nu++) {
1907                         for (int rho=0; rho<=p; rho++) {
1908                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
1909                                                   * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
1910                         }
1911                 }
1912                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
1913
1914                 return result;
1915         }
1916         else if (x == -1) {
1917                 // [Kol] (2.22)
1918                 if (p == 1) {
1919                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
1920                 }
1921 //              throw std::runtime_error("don't know how to evaluate this function!");
1922         }
1923
1924         // what is the desired float format?
1925         // first guess: default format
1926         cln::float_format_t prec = cln::default_float_format;
1927         const cln::cl_N value = x;
1928         // second guess: the argument's format
1929         if (!instanceof(realpart(value), cln::cl_RA_ring))
1930                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
1931         else if (!instanceof(imagpart(value), cln::cl_RA_ring))
1932                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
1933
1934         // [Kol] (5.3)
1935         if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95))) {
1936
1937                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
1938                                    * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
1939
1940                 for (int s=0; s<n; s++) {
1941                         cln::cl_N res2;
1942                         for (int r=0; r<p; r++) {
1943                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
1944                                               * S_num(p-r,n-s,1-value) / cln::factorial(r);
1945                         }
1946                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
1947                 }
1948
1949                 return result;
1950
1951         }
1952         // [Kol] (5.12)
1953         if (cln::abs(value) > 1) {
1954
1955                 cln::cl_N result;
1956
1957                 for (int s=0; s<p; s++) {
1958                         for (int r=0; r<=s; r++) {
1959                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
1960                                                   / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
1961                                                   * S_num(n+s-r,p-s,cln::recip(value));
1962                         }
1963                 }
1964                 result = result * cln::expt(cln::cl_I(-1),n);
1965
1966                 cln::cl_N res2;
1967                 for (int r=0; r<n; r++) {
1968                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
1969                 }
1970                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
1971
1972                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
1973
1974                 return result;
1975         }
1976         else {
1977                 return S_projection(n, p, value, prec);
1978         }
1979 }
1980
1981
1982 } // end of anonymous namespace
1983
1984
1985 //////////////////////////////////////////////////////////////////////
1986 //
1987 // Nielsen's generalized polylogarithm  S(n,p,x)
1988 //
1989 // GiNaC function
1990 //
1991 //////////////////////////////////////////////////////////////////////
1992
1993
1994 static ex S_evalf(const ex& n, const ex& p, const ex& x)
1995 {
1996         if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
1997                 const int n_ = ex_to<numeric>(n).to_int();
1998                 const int p_ = ex_to<numeric>(p).to_int();
1999                 if (is_a<numeric>(x)) {
2000                         const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
2001                         const cln::cl_N result = S_num(n_, p_, x_);
2002                         return numeric(result);
2003                 } else {
2004                         ex x_val = x.evalf();
2005                         if (is_a<numeric>(x_val)) {
2006                                 const cln::cl_N x_val_ = ex_to<numeric>(x_val).to_cl_N();
2007                                 const cln::cl_N result = S_num(n_, p_, x_val_);
2008                                 return numeric(result);
2009                         }
2010                 }
2011         }
2012         return S(n, p, x).hold();
2013 }
2014
2015
2016 static ex S_eval(const ex& n, const ex& p, const ex& x)
2017 {
2018         if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
2019                 if (x == 0) {
2020                         return _ex0;
2021                 }
2022                 if (x == 1) {
2023                         lst m(n+1);
2024                         for (int i=ex_to<numeric>(p).to_int()-1; i>0; i--) {
2025                                 m.append(1);
2026                         }
2027                         return zeta(m);
2028                 }
2029                 if (p == 1) {
2030                         return Li(n+1, x);
2031                 }
2032                 if (x.info(info_flags::numeric) && (!x.info(info_flags::crational))) {
2033                         int n_ = ex_to<numeric>(n).to_int();
2034                         int p_ = ex_to<numeric>(p).to_int();
2035                         const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
2036                         const cln::cl_N result = S_num(n_, p_, x_);
2037                         return numeric(result);
2038                 }
2039         }
2040         if (n.is_zero()) {
2041                 // [Kol] (5.3)
2042                 return pow(-log(1-x), p) / factorial(p);
2043         }
2044         return S(n, p, x).hold();
2045 }
2046
2047
2048 static ex S_series(const ex& n, const ex& p, const ex& x, const relational& rel, int order, unsigned options)
2049 {
2050         if (p == _ex1) {
2051                 return Li(n+1, x).series(rel, order, options);
2052         }
2053
2054         const ex x_pt = x.subs(rel, subs_options::no_pattern);
2055         if (n.info(info_flags::posint) && p.info(info_flags::posint) && x_pt.info(info_flags::numeric)) {
2056                 // First special case: x==0 (derivatives have poles)
2057                 if (x_pt.is_zero()) {
2058                         const symbol s;
2059                         ex ser;
2060                         // manually construct the primitive expansion
2061                         // subsum = Euler-Zagier-Sum is needed
2062                         // dirty hack (slow ...) calculation of subsum:
2063                         std::vector<ex> presubsum, subsum;
2064                         subsum.push_back(0);
2065                         for (int i=1; i<order-1; ++i) {
2066                                 subsum.push_back(subsum[i-1] + numeric(1, i));
2067                         }
2068                         for (int depth=2; depth<p; ++depth) {
2069                                 presubsum = subsum;
2070                                 for (int i=1; i<order-1; ++i) {
2071                                         subsum[i] = subsum[i-1] + numeric(1, i) * presubsum[i-1];
2072                                 }
2073                         }
2074
2075                         for (int i=1; i<order; ++i) {
2076                                 ser += pow(s,i) / pow(numeric(i), n+1) * subsum[i-1];
2077                         }
2078                         // substitute the argument's series expansion
2079                         ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
2080                         // maybe that was terminating, so add a proper order term
2081                         epvector nseq;
2082                         nseq.push_back(expair(Order(_ex1), order));
2083                         ser += pseries(rel, nseq);
2084                         // reexpanding it will collapse the series again
2085                         return ser.series(rel, order);
2086                 }
2087                 // TODO special cases: x==1 (branch point) and x real, >=1 (branch cut)
2088                 throw std::runtime_error("S_series: don't know how to do the series expansion at this point!");
2089         }
2090         // all other cases should be safe, by now:
2091         throw do_taylor();  // caught by function::series()
2092 }
2093
2094
2095 static ex S_deriv(const ex& n, const ex& p, const ex& x, unsigned deriv_param)
2096 {
2097         GINAC_ASSERT(deriv_param < 3);
2098         if (deriv_param < 2) {
2099                 return _ex0;
2100         }
2101         if (n > 0) {
2102                 return S(n-1, p, x) / x;
2103         } else {
2104                 return S(n, p-1, x) / (1-x);
2105         }
2106 }
2107
2108
2109 static void S_print_latex(const ex& n, const ex& p, const ex& x, const print_context& c)
2110 {
2111         c.s << "\\mbox{S}_{";
2112         n.print(c);
2113         c.s << ",";
2114         p.print(c);
2115         c.s << "}(";
2116         x.print(c);
2117         c.s << ")";
2118 }
2119
2120
2121 REGISTER_FUNCTION(S,
2122                   evalf_func(S_evalf).
2123                   eval_func(S_eval).
2124                   series_func(S_series).
2125                   derivative_func(S_deriv).
2126                   print_func<print_latex>(S_print_latex).
2127                   do_not_evalf_params());
2128
2129
2130 //////////////////////////////////////////////////////////////////////
2131 //
2132 // Harmonic polylogarithm  H(m,x)
2133 //
2134 // helper functions
2135 //
2136 //////////////////////////////////////////////////////////////////////
2137
2138
2139 // anonymous namespace for helper functions
2140 namespace {
2141
2142
2143 // regulates the pole (used by 1/x-transformation)
2144 symbol H_polesign("IMSIGN");
2145
2146
2147 // convert parameters from H to Li representation
2148 // parameters are expected to be in expanded form, i.e. only 0, 1 and -1
2149 // returns true if some parameters are negative
2150 bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
2151 {
2152         // expand parameter list
2153         lst mexp;
2154         for (lst::const_iterator it = l.begin(); it != l.end(); it++) {
2155                 if (*it > 1) {
2156                         for (ex count=*it-1; count > 0; count--) {
2157                                 mexp.append(0);
2158                         }
2159                         mexp.append(1);
2160                 } else if (*it < -1) {
2161                         for (ex count=*it+1; count < 0; count++) {
2162                                 mexp.append(0);
2163                         }
2164                         mexp.append(-1);
2165                 } else {
2166                         mexp.append(*it);
2167                 }
2168         }
2169
2170         ex signum = 1;
2171         pf = 1;
2172         bool has_negative_parameters = false;
2173         ex acc = 1;
2174         for (lst::const_iterator it = mexp.begin(); it != mexp.end(); it++) {
2175                 if (*it == 0) {
2176                         acc++;
2177                         continue;
2178                 }
2179                 if (*it > 0) {
2180                         m.append((*it+acc-1) * signum);
2181                 } else {
2182                         m.append((*it-acc+1) * signum);
2183                 }
2184                 acc = 1;
2185                 signum = *it;
2186                 pf *= *it;
2187                 if (pf < 0) {
2188                         has_negative_parameters = true;
2189                 }
2190         }
2191         if (has_negative_parameters) {
2192                 for (int i=0; i<m.nops(); i++) {
2193                         if (m.op(i) < 0) {
2194                                 m.let_op(i) = -m.op(i);
2195                                 s.append(-1);
2196                         } else {
2197                                 s.append(1);
2198                         }
2199                 }
2200         }
2201
2202         return has_negative_parameters;
2203 }
2204
2205
2206 // recursivly transforms H to corresponding multiple polylogarithms
2207 struct map_trafo_H_convert_to_Li : public map_function
2208 {
2209         ex operator()(const ex& e)
2210         {
2211                 if (is_a<add>(e) || is_a<mul>(e)) {
2212                         return e.map(*this);
2213                 }
2214                 if (is_a<function>(e)) {
2215                         std::string name = ex_to<function>(e).get_name();
2216                         if (name == "H") {
2217                                 lst parameter;
2218                                 if (is_a<lst>(e.op(0))) {
2219                                                 parameter = ex_to<lst>(e.op(0));
2220                                 } else {
2221                                         parameter = lst(e.op(0));
2222                                 }
2223                                 ex arg = e.op(1);
2224
2225                                 lst m;
2226                                 lst s;
2227                                 ex pf;
2228                                 if (convert_parameter_H_to_Li(parameter, m, s, pf)) {
2229                                         s.let_op(0) = s.op(0) * arg;
2230                                         return pf * Li(m, s).hold();
2231                                 } else {
2232                                         for (int i=0; i<m.nops(); i++) {
2233                                                 s.append(1);
2234                                         }
2235                                         s.let_op(0) = s.op(0) * arg;
2236                                         return Li(m, s).hold();
2237                                 }
2238                         }
2239                 }
2240                 return e;
2241         }
2242 };
2243
2244
2245 // recursivly transforms H to corresponding zetas
2246 struct map_trafo_H_convert_to_zeta : public map_function
2247 {
2248         ex operator()(const ex& e)
2249         {
2250                 if (is_a<add>(e) || is_a<mul>(e)) {
2251                         return e.map(*this);
2252                 }
2253                 if (is_a<function>(e)) {
2254                         std::string name = ex_to<function>(e).get_name();
2255                         if (name == "H") {
2256                                 lst parameter;
2257                                 if (is_a<lst>(e.op(0))) {
2258                                                 parameter = ex_to<lst>(e.op(0));
2259                                 } else {
2260                                         parameter = lst(e.op(0));
2261                                 }
2262
2263                                 lst m;
2264                                 lst s;
2265                                 ex pf;
2266                                 if (convert_parameter_H_to_Li(parameter, m, s, pf)) {
2267                                         return pf * zeta(m, s);
2268                                 } else {
2269                                         return zeta(m);
2270                                 }
2271                         }
2272                 }
2273                 return e;
2274         }
2275 };
2276
2277
2278 // remove trailing zeros from H-parameters
2279 struct map_trafo_H_reduce_trailing_zeros : public map_function
2280 {
2281         ex operator()(const ex& e)
2282         {
2283                 if (is_a<add>(e) || is_a<mul>(e)) {
2284                         return e.map(*this);
2285                 }
2286                 if (is_a<function>(e)) {
2287                         std::string name = ex_to<function>(e).get_name();
2288                         if (name == "H") {
2289                                 lst parameter;
2290                                 if (is_a<lst>(e.op(0))) {
2291                                         parameter = ex_to<lst>(e.op(0));
2292                                 } else {
2293                                         parameter = lst(e.op(0));
2294                                 }
2295                                 ex arg = e.op(1);
2296                                 if (parameter.op(parameter.nops()-1) == 0) {
2297
2298                                         //
2299                                         if (parameter.nops() == 1) {
2300                                                 return log(arg);
2301                                         }
2302
2303                                         //
2304                                         lst::const_iterator it = parameter.begin();
2305                                         while ((it != parameter.end()) && (*it == 0)) {
2306                                                 it++;
2307                                         }
2308                                         if (it == parameter.end()) {
2309                                                 return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
2310                                         }
2311
2312                                         //
2313                                         parameter.remove_last();
2314                                         int lastentry = parameter.nops();
2315                                         while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
2316                                                 lastentry--;
2317                                         }
2318
2319                                         //
2320                                         ex result = log(arg) * H(parameter,arg).hold();
2321                                         ex acc = 0;
2322                                         for (ex i=0; i<lastentry; i++) {
2323                                                 if (parameter[i] > 0) {
2324                                                         parameter[i]++;
2325                                                         result -= (acc + parameter[i]-1) * H(parameter, arg).hold();
2326                                                         parameter[i]--;
2327                                                         acc = 0;
2328                                                 } else if (parameter[i] < 0) {
2329                                                         parameter[i]--;
2330                                                         result -= (acc + abs(parameter[i]+1)) * H(parameter, arg).hold();
2331                                                         parameter[i]++;
2332                                                         acc = 0;
2333                                                 } else {
2334                                                         acc++;
2335                                                 }
2336                                         }
2337
2338                                         if (lastentry < parameter.nops()) {
2339                                                 result = result / (parameter.nops()-lastentry+1);
2340                                                 return result.map(*this);
2341                                         } else {
2342                                                 return result;
2343                                         }
2344                                 }
2345                         }
2346                 }
2347                 return e;
2348         }
2349 };
2350
2351
2352 // returns an expression with zeta functions corresponding to the parameter list for H
2353 ex convert_H_to_zeta(const lst& m)
2354 {
2355         symbol xtemp("xtemp");
2356         map_trafo_H_reduce_trailing_zeros filter;
2357         map_trafo_H_convert_to_zeta filter2;
2358         return filter2(filter(H(m, xtemp).hold())).subs(xtemp == 1);
2359 }
2360
2361
2362 // convert signs form Li to H representation
2363 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf)
2364 {
2365         lst res;
2366         lst::const_iterator itm = m.begin();
2367         lst::const_iterator itx = ++x.begin();
2368         int signum = 1;
2369         pf = _ex1;
2370         res.append(*itm);
2371         itm++;
2372         while (itx != x.end()) {
2373                 signum *= (*itx > 0) ? 1 : -1;
2374                 pf *= signum;
2375                 res.append((*itm) * signum);
2376                 itm++;
2377                 itx++;
2378         }
2379         return res;
2380 }
2381
2382
2383 // multiplies an one-dimensional H with another H
2384 // [ReV] (18)
2385 ex trafo_H_mult(const ex& h1, const ex& h2)
2386 {
2387         ex res;
2388         ex hshort;
2389         lst hlong;
2390         ex h1nops = h1.op(0).nops();
2391         ex h2nops = h2.op(0).nops();
2392         if (h1nops > 1) {
2393                 hshort = h2.op(0).op(0);
2394                 hlong = ex_to<lst>(h1.op(0));
2395         } else {
2396                 hshort = h1.op(0).op(0);
2397                 if (h2nops > 1) {
2398                         hlong = ex_to<lst>(h2.op(0));
2399                 } else {
2400                         hlong = h2.op(0).op(0);
2401                 }
2402         }
2403         for (int i=0; i<=hlong.nops(); i++) {
2404                 lst newparameter;
2405                 int j=0;
2406                 for (; j<i; j++) {
2407                         newparameter.append(hlong[j]);
2408                 }
2409                 newparameter.append(hshort);
2410                 for (; j<hlong.nops(); j++) {
2411                         newparameter.append(hlong[j]);
2412                 }
2413                 res += H(newparameter, h1.op(1)).hold();
2414         }
2415         return res;
2416 }
2417
2418
2419 // applies trafo_H_mult recursively on expressions
2420 struct map_trafo_H_mult : public map_function
2421 {
2422         ex operator()(const ex& e)
2423         {
2424                 if (is_a<add>(e)) {
2425                         return e.map(*this);
2426                 }
2427
2428                 if (is_a<mul>(e)) {
2429
2430                         ex result = 1;
2431                         ex firstH;
2432                         lst Hlst;
2433                         for (int pos=0; pos<e.nops(); pos++) {
2434                                 if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
2435                                         std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
2436                                         if (name == "H") {
2437                                                 for (ex i=0; i<e.op(pos).op(1); i++) {
2438                                                         Hlst.append(e.op(pos).op(0));
2439                                                 }
2440                                                 continue;
2441                                         }
2442                                 } else if (is_a<function>(e.op(pos))) {
2443                                         std::string name = ex_to<function>(e.op(pos)).get_name();
2444                                         if (name == "H") {
2445                                                 if (e.op(pos).op(0).nops() > 1) {
2446                                                         firstH = e.op(pos);
2447                                                 } else {
2448                                                         Hlst.append(e.op(pos));
2449                                                 }
2450                                                 continue;
2451                                         }
2452                                 }
2453                                 result *= e.op(pos);
2454                         }
2455                         if (firstH == 0) {
2456                                 if (Hlst.nops() > 0) {
2457                                         firstH = Hlst[Hlst.nops()-1];
2458                                         Hlst.remove_last();
2459                                 } else {
2460                                         return e;
2461                                 }
2462                         }
2463
2464                         if (Hlst.nops() > 0) {
2465                                 ex buffer = trafo_H_mult(firstH, Hlst.op(0));
2466                                 result *= buffer;
2467                                 for (int i=1; i<Hlst.nops(); i++) {
2468                                         result *= Hlst.op(i);
2469                                 }
2470                                 result = result.expand();
2471                                 map_trafo_H_mult recursion;
2472                                 return recursion(result);
2473                         } else {
2474                                 return e;
2475                         }
2476
2477                 }
2478                 return e;
2479         }
2480 };
2481
2482
2483 // do integration [ReV] (55)
2484 // put parameter 0 in front of existing parameters
2485 ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
2486 {
2487         ex h;
2488         std::string name;
2489         if (is_a<function>(e)) {
2490                 name = ex_to<function>(e).get_name();
2491         }
2492         if (name == "H") {
2493                 h = e;
2494         } else {
2495                 for (int i=0; i<e.nops(); i++) {
2496                         if (is_a<function>(e.op(i))) {
2497                                 std::string name = ex_to<function>(e.op(i)).get_name();
2498                                 if (name == "H") {
2499                                         h = e.op(i);
2500                                 }
2501                         }
2502                 }
2503         }
2504         if (h != 0) {
2505                 lst newparameter = ex_to<lst>(h.op(0));
2506                 newparameter.prepend(0);
2507                 ex addzeta = convert_H_to_zeta(newparameter);
2508                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
2509         } else {
2510                 return e * (-H(lst(0),1/arg).hold());
2511         }
2512 }
2513
2514
2515 // do integration [ReV] (49)
2516 // put parameter 1 in front of existing parameters
2517 ex trafo_H_prepend_one(const ex& e, const ex& arg)
2518 {
2519         ex h;
2520         std::string name;
2521         if (is_a<function>(e)) {
2522                 name = ex_to<function>(e).get_name();
2523         }
2524         if (name == "H") {
2525                 h = e;
2526         } else {
2527                 for (int i=0; i<e.nops(); i++) {
2528                         if (is_a<function>(e.op(i))) {
2529                                 std::string name = ex_to<function>(e.op(i)).get_name();
2530                                 if (name == "H") {
2531                                         h = e.op(i);
2532                                 }
2533                         }
2534                 }
2535         }
2536         if (h != 0) {
2537                 lst newparameter = ex_to<lst>(h.op(0));
2538                 newparameter.prepend(1);
2539                 return e.subs(h == H(newparameter, h.op(1)).hold());
2540         } else {
2541                 return e * H(lst(1),1-arg).hold();
2542         }
2543 }
2544
2545
2546 // do integration [ReV] (55)
2547 // put parameter -1 in front of existing parameters
2548 ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
2549 {
2550         ex h;
2551         std::string name;
2552         if (is_a<function>(e)) {
2553                 name = ex_to<function>(e).get_name();
2554         }
2555         if (name == "H") {
2556                 h = e;
2557         } else {
2558                 for (int i=0; i<e.nops(); i++) {
2559                         if (is_a<function>(e.op(i))) {
2560                                 std::string name = ex_to<function>(e.op(i)).get_name();
2561                                 if (name == "H") {
2562                                         h = e.op(i);
2563                                 }
2564                         }
2565                 }
2566         }
2567         if (h != 0) {
2568                 lst newparameter = ex_to<lst>(h.op(0));
2569                 newparameter.prepend(-1);
2570                 ex addzeta = convert_H_to_zeta(newparameter);
2571                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
2572         } else {
2573                 ex addzeta = convert_H_to_zeta(lst(-1));
2574                 return (e * (addzeta - H(lst(-1),1/arg).hold())).expand();
2575         }
2576 }
2577
2578
2579 // do integration [ReV] (55)
2580 // put parameter -1 in front of existing parameters
2581 ex trafo_H_1mxt1px_prepend_minusone(const ex& e, const ex& arg)
2582 {
2583         ex h;
2584         std::string name;
2585         if (is_a<function>(e)) {
2586                 name = ex_to<function>(e).get_name();
2587         }
2588         if (name == "H") {
2589                 h = e;
2590         } else {
2591                 for (int i=0; i<e.nops(); i++) {
2592                         if (is_a<function>(e.op(i))) {
2593                                 std::string name = ex_to<function>(e.op(i)).get_name();
2594                                 if (name == "H") {
2595                                         h = e.op(i);
2596                                 }
2597                         }
2598                 }
2599         }
2600         if (h != 0) {
2601                 lst newparameter = ex_to<lst>(h.op(0));
2602                 newparameter.prepend(-1);
2603                 return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
2604         } else {
2605                 return (e * H(lst(-1),(1-arg)/(1+arg)).hold()).expand();
2606         }
2607 }
2608
2609
2610 // do integration [ReV] (55)
2611 // put parameter 1 in front of existing parameters
2612 ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
2613 {
2614         ex h;
2615         std::string name;
2616         if (is_a<function>(e)) {
2617                 name = ex_to<function>(e).get_name();
2618         }
2619         if (name == "H") {
2620                 h = e;
2621         } else {
2622                 for (int i=0; i<e.nops(); i++) {
2623                         if (is_a<function>(e.op(i))) {
2624                                 std::string name = ex_to<function>(e.op(i)).get_name();
2625                                 if (name == "H") {
2626                                         h = e.op(i);
2627                                 }
2628                         }
2629                 }
2630         }
2631         if (h != 0) {
2632                 lst newparameter = ex_to<lst>(h.op(0));
2633                 newparameter.prepend(1);
2634                 return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
2635         } else {
2636                 return (e * H(lst(1),(1-arg)/(1+arg)).hold()).expand();
2637         }
2638 }
2639
2640
2641 // do x -> 1-x transformation
2642 struct map_trafo_H_1mx : public map_function
2643 {
2644         ex operator()(const ex& e)
2645         {
2646                 if (is_a<add>(e) || is_a<mul>(e)) {
2647                         return e.map(*this);
2648                 }
2649
2650                 if (is_a<function>(e)) {
2651                         std::string name = ex_to<function>(e).get_name();
2652                         if (name == "H") {
2653
2654                                 lst parameter = ex_to<lst>(e.op(0));
2655                                 ex arg = e.op(1);
2656
2657                                 // special cases if all parameters are either 0, 1 or -1
2658                                 bool allthesame = true;
2659                                 if (parameter.op(0) == 0) {
2660                                         for (int i=1; i<parameter.nops(); i++) {
2661                                                 if (parameter.op(i) != 0) {
2662                                                         allthesame = false;
2663                                                         break;
2664                                                 }
2665                                         }
2666                                         if (allthesame) {
2667                                                 lst newparameter;
2668                                                 for (int i=parameter.nops(); i>0; i--) {
2669                                                         newparameter.append(1);
2670                                                 }
2671                                                 return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
2672                                         }
2673                                 } else if (parameter.op(0) == -1) {
2674                                         throw std::runtime_error("map_trafo_H_1mx: cannot handle weights equal -1!");
2675                                 } else {
2676                                         for (int i=1; i<parameter.nops(); i++) {
2677                                                 if (parameter.op(i) != 1) {
2678                                                         allthesame = false;
2679                                                         break;
2680                                                 }
2681                                         }
2682                                         if (allthesame) {
2683                                                 lst newparameter;
2684                                                 for (int i=parameter.nops(); i>0; i--) {
2685                                                         newparameter.append(0);
2686                                                 }
2687                                                 return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
2688                                         }
2689                                 }
2690
2691                                 lst newparameter = parameter;
2692                                 newparameter.remove_first();
2693
2694                                 if (parameter.op(0) == 0) {
2695
2696                                         // leading zero
2697                                         ex res = convert_H_to_zeta(parameter);
2698                                         //ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
2699                                         map_trafo_H_1mx recursion;
2700                                         ex buffer = recursion(H(newparameter, arg).hold());
2701                                         if (is_a<add>(buffer)) {
2702                                                 for (int i=0; i<buffer.nops(); i++) {
2703                                                         res -= trafo_H_prepend_one(buffer.op(i), arg);
2704                                                 }
2705                                         } else {
2706                                                 res -= trafo_H_prepend_one(buffer, arg);
2707                                         }
2708                                         return res;
2709
2710                                 } else {
2711
2712                                         // leading one
2713                                         map_trafo_H_1mx recursion;
2714                                         map_trafo_H_mult unify;
2715                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
2716                                         int firstzero = 0;
2717                                         while (parameter.op(firstzero) == 1) {
2718                                                 firstzero++;
2719                                         }
2720                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
2721                                                 lst newparameter;
2722                                                 int j=0;
2723                                                 for (; j<=i; j++) {
2724                                                         newparameter.append(parameter[j+1]);
2725                                                 }
2726                                                 newparameter.append(1);
2727                                                 for (; j<parameter.nops()-1; j++) {
2728                                                         newparameter.append(parameter[j+1]);
2729                                                 }
2730                                                 res -= H(newparameter, arg).hold();
2731                                         }
2732                                         res = recursion(res).expand() / firstzero;
2733                                         return unify(res);
2734                                 }
2735                         }
2736                 }
2737                 return e;
2738         }
2739 };
2740
2741
2742 // do x -> 1/x transformation
2743 struct map_trafo_H_1overx : public map_function
2744 {
2745         ex operator()(const ex& e)
2746         {
2747                 if (is_a<add>(e) || is_a<mul>(e)) {
2748                         return e.map(*this);
2749                 }
2750
2751                 if (is_a<function>(e)) {
2752                         std::string name = ex_to<function>(e).get_name();
2753                         if (name == "H") {
2754
2755                                 lst parameter = ex_to<lst>(e.op(0));
2756                                 ex arg = e.op(1);
2757
2758                                 // special cases if all parameters are either 0, 1 or -1
2759                                 bool allthesame = true;
2760                                 if (parameter.op(0) == 0) {
2761                                         for (int i=1; i<parameter.nops(); i++) {
2762                                                 if (parameter.op(i) != 0) {
2763                                                         allthesame = false;
2764                                                         break;
2765                                                 }
2766                                         }
2767                                         if (allthesame) {
2768                                                 return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
2769                                         }
2770                                 } else if (parameter.op(0) == -1) {
2771                                         for (int i=1; i<parameter.nops(); i++) {
2772                                                 if (parameter.op(i) != -1) {
2773                                                         allthesame = false;
2774                                                         break;
2775                                                 }
2776                                         }
2777                                         if (allthesame) {
2778                                                 map_trafo_H_mult unify;
2779                                                 return unify((pow(H(lst(-1),1/arg).hold() - H(lst(0),1/arg).hold(), parameter.nops())
2780                                                        / factorial(parameter.nops())).expand());
2781                                         }
2782                                 } else {
2783                                         for (int i=1; i<parameter.nops(); i++) {
2784                                                 if (parameter.op(i) != 1) {
2785                                                         allthesame = false;
2786                                                         break;
2787                                                 }
2788                                         }
2789                                         if (allthesame) {
2790                                                 map_trafo_H_mult unify;
2791                                                 return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() + H_polesign, parameter.nops())
2792                                                        / factorial(parameter.nops())).expand());
2793                                         }
2794                                 }
2795
2796                                 lst newparameter = parameter;
2797                                 newparameter.remove_first();
2798
2799                                 if (parameter.op(0) == 0) {
2800
2801                                         // leading zero
2802                                         ex res = convert_H_to_zeta(parameter);
2803                                         map_trafo_H_1overx recursion;
2804                                         ex buffer = recursion(H(newparameter, arg).hold());
2805                                         if (is_a<add>(buffer)) {
2806                                                 for (int i=0; i<buffer.nops(); i++) {
2807                                                         res += trafo_H_1tx_prepend_zero(buffer.op(i), arg);
2808                                                 }
2809                                         } else {
2810                                                 res += trafo_H_1tx_prepend_zero(buffer, arg);
2811                                         }
2812                                         return res;
2813
2814                                 } else if (parameter.op(0) == -1) {
2815
2816                                         // leading negative one
2817                                         ex res = convert_H_to_zeta(parameter);
2818                                         map_trafo_H_1overx recursion;
2819                                         ex buffer = recursion(H(newparameter, arg).hold());
2820                                         if (is_a<add>(buffer)) {
2821                                                 for (int i=0; i<buffer.nops(); i++) {
2822                                                         res += trafo_H_1tx_prepend_zero(buffer.op(i), arg) - trafo_H_1tx_prepend_minusone(buffer.op(i), arg);
2823                                                 }
2824                                         } else {
2825                                                 res += trafo_H_1tx_prepend_zero(buffer, arg) - trafo_H_1tx_prepend_minusone(buffer, arg);
2826                                         }
2827                                         return res;
2828
2829                                 } else {
2830
2831                                         // leading one
2832                                         map_trafo_H_1overx recursion;
2833                                         map_trafo_H_mult unify;
2834                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
2835                                         int firstzero = 0;
2836                                         while (parameter.op(firstzero) == 1) {
2837                                                 firstzero++;
2838                                         }
2839                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
2840                                                 lst newparameter;
2841                                                 int j=0;
2842                                                 for (; j<=i; j++) {
2843                                                         newparameter.append(parameter[j+1]);
2844                                                 }
2845                                                 newparameter.append(1);
2846                                                 for (; j<parameter.nops()-1; j++) {
2847                                                         newparameter.append(parameter[j+1]);
2848                                                 }
2849                                                 res -= H(newparameter, arg).hold();
2850                                         }
2851                                         res = recursion(res).expand() / firstzero;
2852                                         return unify(res);
2853
2854                                 }
2855
2856                         }
2857                 }
2858                 return e;
2859         }
2860 };
2861
2862
2863 // do x -> (1-x)/(1+x) transformation
2864 struct map_trafo_H_1mxt1px : public map_function
2865 {
2866         ex operator()(const ex& e)
2867         {
2868                 if (is_a<add>(e) || is_a<mul>(e)) {
2869                         return e.map(*this);
2870                 }
2871
2872                 if (is_a<function>(e)) {
2873                         std::string name = ex_to<function>(e).get_name();
2874                         if (name == "H") {
2875
2876                                 lst parameter = ex_to<lst>(e.op(0));
2877                                 ex arg = e.op(1);
2878
2879                                 // special cases if all parameters are either 0, 1 or -1
2880                                 bool allthesame = true;
2881                                 if (parameter.op(0) == 0) {
2882                                         for (int i=1; i<parameter.nops(); i++) {
2883                                                 if (parameter.op(i) != 0) {
2884                                                         allthesame = false;
2885                                                         break;
2886                                                 }
2887                                         }
2888                                         if (allthesame) {
2889                                                 map_trafo_H_mult unify;
2890                                                 return unify((pow(-H(lst(1),(1-arg)/(1+arg)).hold() - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
2891                                                        / factorial(parameter.nops())).expand());
2892                                         }
2893                                 } else if (parameter.op(0) == -1) {
2894                                         for (int i=1; i<parameter.nops(); i++) {
2895                                                 if (parameter.op(i) != -1) {
2896                                                         allthesame = false;
2897                                                         break;
2898                                                 }
2899                                         }
2900                                         if (allthesame) {
2901                                                 map_trafo_H_mult unify;
2902                                                 return unify((pow(log(2) - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
2903                                                        / factorial(parameter.nops())).expand());
2904                                         }
2905                                 } else {
2906                                         for (int i=1; i<parameter.nops(); i++) {
2907                                                 if (parameter.op(i) != 1) {
2908                                                         allthesame = false;
2909                                                         break;
2910                                                 }
2911                                         }
2912                                         if (allthesame) {
2913                                                 map_trafo_H_mult unify;
2914                                                 return unify((pow(-log(2) - H(lst(0),(1-arg)/(1+arg)).hold() + H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
2915                                                        / factorial(parameter.nops())).expand());
2916                                         }
2917                                 }
2918
2919                                 lst newparameter = parameter;
2920                                 newparameter.remove_first();
2921
2922                                 if (parameter.op(0) == 0) {
2923
2924                                         // leading zero
2925                                         ex res = convert_H_to_zeta(parameter);
2926                                         map_trafo_H_1mxt1px recursion;
2927                                         ex buffer = recursion(H(newparameter, arg).hold());
2928                                         if (is_a<add>(buffer)) {
2929                                                 for (int i=0; i<buffer.nops(); i++) {
2930                                                         res -= trafo_H_1mxt1px_prepend_one(buffer.op(i), arg) + trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
2931                                                 }
2932                                         } else {
2933                                                 res -= trafo_H_1mxt1px_prepend_one(buffer, arg) + trafo_H_1mxt1px_prepend_minusone(buffer, arg);
2934                                         }
2935                                         return res;
2936
2937                                 } else if (parameter.op(0) == -1) {
2938
2939                                         // leading negative one
2940                                         ex res = convert_H_to_zeta(parameter);
2941                                         map_trafo_H_1mxt1px recursion;
2942                                         ex buffer = recursion(H(newparameter, arg).hold());
2943                                         if (is_a<add>(buffer)) {
2944                                                 for (int i=0; i<buffer.nops(); i++) {
2945                                                         res -= trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
2946                                                 }
2947                                         } else {
2948                                                 res -= trafo_H_1mxt1px_prepend_minusone(buffer, arg);
2949                                         }
2950                                         return res;
2951
2952                                 } else {
2953
2954                                         // leading one
2955                                         map_trafo_H_1mxt1px recursion;
2956                                         map_trafo_H_mult unify;
2957                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
2958                                         int firstzero = 0;
2959                                         while (parameter.op(firstzero) == 1) {
2960                                                 firstzero++;
2961                                         }
2962                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
2963                                                 lst newparameter;
2964                                                 int j=0;
2965                                                 for (; j<=i; j++) {
2966                                                         newparameter.append(parameter[j+1]);
2967                                                 }
2968                                                 newparameter.append(1);
2969                                                 for (; j<parameter.nops()-1; j++) {
2970                                                         newparameter.append(parameter[j+1]);
2971                                                 }
2972                                                 res -= H(newparameter, arg).hold();
2973                                         }
2974                                         res = recursion(res).expand() / firstzero;
2975                                         return unify(res);
2976
2977                                 }
2978
2979                         }
2980                 }
2981                 return e;
2982         }
2983 };
2984
2985
2986 // do the actual summation.
2987 cln::cl_N H_do_sum(const std::vector<int>& m, const cln::cl_N& x)
2988 {
2989         const int j = m.size();
2990
2991         std::vector<cln::cl_N> t(j);
2992
2993         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
2994         cln::cl_N factor = cln::expt(x, j) * one;
2995         cln::cl_N t0buf;
2996         int q = 0;
2997         do {
2998                 t0buf = t[0];
2999                 q++;
3000                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),m[j-1]);
3001                 for (int k=j-2; k>=1; k--) {
3002                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), m[k]);
3003                 }
3004                 t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), m[0]);
3005                 factor = factor * x;
3006         } while (t[0] != t0buf);
3007
3008         return t[0];
3009 }
3010
3011
3012 } // end of anonymous namespace
3013
3014
3015 //////////////////////////////////////////////////////////////////////
3016 //
3017 // Harmonic polylogarithm  H(m,x)
3018 //
3019 // GiNaC function
3020 //
3021 //////////////////////////////////////////////////////////////////////
3022
3023
3024 static ex H_evalf(const ex& x1, const ex& x2)
3025 {
3026         if (is_a<lst>(x1)) {
3027
3028                 cln::cl_N x;
3029                 if (is_a<numeric>(x2)) {
3030                         x = ex_to<numeric>(x2).to_cl_N();
3031                 } else {
3032                         ex x2_val = x2.evalf();
3033                         if (is_a<numeric>(x2_val)) {
3034                                 x = ex_to<numeric>(x2_val).to_cl_N();
3035                         }
3036                 }
3037
3038                 for (int i=0; i<x1.nops(); i++) {
3039                         if (!x1.op(i).info(info_flags::integer)) {
3040                                 return H(x1, x2).hold();
3041                         }
3042                 }
3043                 if (x1.nops() < 1) {
3044                         return H(x1, x2).hold();
3045                 }
3046
3047                 const lst& morg = ex_to<lst>(x1);
3048                 // remove trailing zeros ...
3049                 if (*(--morg.end()) == 0) {
3050                         symbol xtemp("xtemp");
3051                         map_trafo_H_reduce_trailing_zeros filter;
3052                         return filter(H(x1, xtemp).hold()).subs(xtemp==x2).evalf();
3053                 }
3054                 // ... and expand parameter notation
3055                 bool has_minus_one = false;
3056                 lst m;
3057                 for (lst::const_iterator it = morg.begin(); it != morg.end(); it++) {
3058                         if (*it > 1) {
3059                                 for (ex count=*it-1; count > 0; count--) {
3060                                         m.append(0);
3061                                 }
3062                                 m.append(1);
3063                         } else if (*it <= -1) {
3064                                 for (ex count=*it+1; count < 0; count++) {
3065                                         m.append(0);
3066                                 }
3067                                 m.append(-1);
3068                                 has_minus_one = true;
3069                         } else {
3070                                 m.append(*it);
3071                         }
3072                 }
3073
3074                 // do summation
3075                 if (cln::abs(x) < 0.95) {
3076                         lst m_lst;
3077                         lst s_lst;
3078                         ex pf;
3079                         if (convert_parameter_H_to_Li(m, m_lst, s_lst, pf)) {
3080                                 // negative parameters -> s_lst is filled
3081                                 std::vector<int> m_int;
3082                                 std::vector<cln::cl_N> x_cln;
3083                                 for (lst::const_iterator it_int = m_lst.begin(), it_cln = s_lst.begin();
3084                                      it_int != m_lst.end(); it_int++, it_cln++) {
3085                                         m_int.push_back(ex_to<numeric>(*it_int).to_int());
3086                                         x_cln.push_back(ex_to<numeric>(*it_cln).to_cl_N());
3087                                 }
3088                                 x_cln.front() = x_cln.front() * x;
3089                                 return pf * numeric(multipleLi_do_sum(m_int, x_cln));
3090                         } else {
3091                                 // only positive parameters
3092                                 //TODO
3093                                 if (m_lst.nops() == 1) {
3094                                         return Li(m_lst.op(0), x2).evalf();
3095                                 }
3096                                 std::vector<int> m_int;
3097                                 for (lst::const_iterator it = m_lst.begin(); it != m_lst.end(); it++) {
3098                                         m_int.push_back(ex_to<numeric>(*it).to_int());
3099                                 }
3100                                 return numeric(H_do_sum(m_int, x));
3101                         }
3102                 }
3103
3104                 symbol xtemp("xtemp");
3105                 ex res = 1;
3106
3107                 // ensure that the realpart of the argument is positive
3108                 if (cln::realpart(x) < 0) {
3109                         x = -x;
3110                         for (int i=0; i<m.nops(); i++) {
3111                                 if (m.op(i) != 0) {
3112                                         m.let_op(i) = -m.op(i);
3113                                         res *= -1;
3114                                 }
3115                         }
3116                 }
3117
3118                 // x -> 1/x
3119                 if (cln::abs(x) >= 2.0) {
3120                         map_trafo_H_1overx trafo;
3121                         res *= trafo(H(m, xtemp));
3122                         if (cln::imagpart(x) <= 0) {
3123                                 res = res.subs(H_polesign == -I*Pi);
3124                         } else {
3125                                 res = res.subs(H_polesign == I*Pi);
3126                         }
3127                         return res.subs(xtemp == numeric(x)).evalf();
3128                 }
3129
3130                 // check transformations for 0.95 <= |x| < 2.0
3131
3132                 // |(1-x)/(1+x)| < 0.9 -> circular area with center=9.53+0i and radius=9.47
3133                 if (cln::abs(x-9.53) <= 9.47) {
3134                         // x -> (1-x)/(1+x)
3135                         map_trafo_H_1mxt1px trafo;
3136                         res *= trafo(H(m, xtemp));
3137                 } else {
3138                         // x -> 1-x
3139                         if (has_minus_one) {
3140                                 map_trafo_H_convert_to_Li filter;
3141                                 return filter(H(m, numeric(x)).hold()).evalf();
3142                         }
3143                         map_trafo_H_1mx trafo;
3144                         res *= trafo(H(m, xtemp));
3145                 }
3146
3147                 return res.subs(xtemp == numeric(x)).evalf();
3148         }
3149
3150         return H(x1,x2).hold();
3151 }
3152
3153
3154 static ex H_eval(const ex& m_, const ex& x)
3155 {
3156         lst m;
3157         if (is_a<lst>(m_)) {
3158                 m = ex_to<lst>(m_);
3159         } else {
3160                 m = lst(m_);
3161         }
3162         if (m.nops() == 0) {
3163                 return _ex1;
3164         }
3165         ex pos1;
3166         ex pos2;
3167         ex n;
3168         ex p;
3169         int step = 0;
3170         if (*m.begin() > _ex1) {
3171                 step++;
3172                 pos1 = _ex0;
3173                 pos2 = _ex1;
3174                 n = *m.begin()-1;
3175                 p = _ex1;
3176         } else if (*m.begin() < _ex_1) {
3177                 step++;
3178                 pos1 = _ex0;
3179                 pos2 = _ex_1;
3180                 n = -*m.begin()-1;
3181                 p = _ex1;
3182         } else if (*m.begin() == _ex0) {
3183                 pos1 = _ex0;
3184                 n = _ex1;
3185         } else {
3186                 pos1 = *m.begin();
3187                 p = _ex1;
3188         }
3189         for (lst::const_iterator it = ++m.begin(); it != m.end(); it++) {
3190                 if ((*it).info(info_flags::integer)) {
3191                         if (step == 0) {
3192                                 if (*it > _ex1) {
3193                                         if (pos1 == _ex0) {
3194                                                 step = 1;
3195                                                 pos2 = _ex1;
3196                                                 n += *it-1;
3197                                                 p = _ex1;
3198                                         } else {
3199                                                 step = 2;
3200                                         }
3201                                 } else if (*it < _ex_1) {
3202                                         if (pos1 == _ex0) {
3203                                                 step = 1;
3204                                                 pos2 = _ex_1;
3205                                                 n += -*it-1;
3206                                                 p = _ex1;
3207                                         } else {
3208                                                 step = 2;
3209                                         }
3210                                 } else {
3211                                         if (*it != pos1) {
3212                                                 step = 1;
3213                                                 pos2 = *it;
3214                                         }
3215                                         if (*it == _ex0) {
3216                                                 n++;
3217                                         } else {
3218                                                 p++;
3219                                         }
3220                                 }
3221                         } else if (step == 1) {
3222                                 if (*it != pos2) {
3223                                         step = 2;
3224                                 } else {
3225                                         if (*it == _ex0) {
3226                                                 n++;
3227                                         } else {
3228                                                 p++;
3229                                         }
3230                                 }
3231                         }
3232                 } else {
3233                         // if some m_i is not an integer
3234                         return H(m_, x).hold();
3235                 }
3236         }
3237         if ((x == _ex1) && (*(--m.end()) != _ex0)) {
3238                 return convert_H_to_zeta(m);
3239         }
3240         if (step == 0) {
3241                 if (pos1 == _ex0) {
3242                         // all zero
3243                         if (x == _ex0) {
3244                                 return H(m_, x).hold();
3245                         }
3246                         return pow(log(x), m.nops()) / factorial(m.nops());
3247                 } else {
3248                         // all (minus) one
3249                         return pow(-pos1*log(1-pos1*x), m.nops()) / factorial(m.nops());
3250                 }
3251         } else if ((step == 1) && (pos1 == _ex0)){
3252                 // convertible to S
3253                 if (pos2 == _ex1) {
3254                         return S(n, p, x);
3255                 } else {
3256                         return pow(-1, p) * S(n, p, -x);
3257                 }
3258         }
3259         if (x == _ex0) {
3260                 return _ex0;
3261         }
3262         if (x.info(info_flags::numeric) && (!x.info(info_flags::crational))) {
3263                 return H(m_, x).evalf();
3264         }
3265         return H(m_, x).hold();
3266 }
3267
3268
3269 static ex H_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
3270 {
3271         epvector seq;
3272         seq.push_back(expair(H(m, x), 0));
3273         return pseries(rel, seq);
3274 }
3275
3276
3277 static ex H_deriv(const ex& m_, const ex& x, unsigned deriv_param)
3278 {
3279         GINAC_ASSERT(deriv_param < 2);
3280         if (deriv_param == 0) {
3281                 return _ex0;
3282         }
3283         lst m;
3284         if (is_a<lst>(m_)) {
3285                 m = ex_to<lst>(m_);
3286         } else {
3287                 m = lst(m_);
3288         }
3289         ex mb = *m.begin();
3290         if (mb > _ex1) {
3291                 m[0]--;
3292                 return H(m, x) / x;
3293         }
3294         if (mb < _ex_1) {
3295                 m[0]++;
3296                 return H(m, x) / x;
3297         }
3298         m.remove_first();
3299         if (mb == _ex1) {
3300                 return 1/(1-x) * H(m, x);
3301         } else if (mb == _ex_1) {
3302                 return 1/(1+x) * H(m, x);
3303         } else {
3304                 return H(m, x) / x;
3305         }
3306 }
3307
3308
3309 static void H_print_latex(const ex& m_, const ex& x, const print_context& c)
3310 {
3311         lst m;
3312         if (is_a<lst>(m_)) {
3313                 m = ex_to<lst>(m_);
3314         } else {
3315                 m = lst(m_);
3316         }
3317         c.s << "\\mbox{H}_{";
3318         lst::const_iterator itm = m.begin();
3319         (*itm).print(c);
3320         itm++;
3321         for (; itm != m.end(); itm++) {
3322                 c.s << ",";
3323                 (*itm).print(c);
3324         }
3325         c.s << "}(";
3326         x.print(c);
3327         c.s << ")";
3328 }
3329
3330
3331 REGISTER_FUNCTION(H,
3332                   evalf_func(H_evalf).
3333                   eval_func(H_eval).
3334                   series_func(H_series).
3335                   derivative_func(H_deriv).
3336                   print_func<print_latex>(H_print_latex).
3337                   do_not_evalf_params());
3338
3339
3340 // takes a parameter list for H and returns an expression with corresponding multiple polylogarithms
3341 ex convert_H_to_Li(const ex& m, const ex& x)
3342 {
3343         map_trafo_H_reduce_trailing_zeros filter;
3344         map_trafo_H_convert_to_Li filter2;
3345         if (is_a<lst>(m)) {
3346                 return filter2(filter(H(m, x).hold()));
3347         } else {
3348                 return filter2(filter(H(lst(m), x).hold()));
3349         }
3350 }
3351
3352
3353 //////////////////////////////////////////////////////////////////////
3354 //
3355 // Multiple zeta values  zeta(x) and zeta(x,s)
3356 //
3357 // helper functions
3358 //
3359 //////////////////////////////////////////////////////////////////////
3360
3361
3362 // anonymous namespace for helper functions
3363 namespace {
3364
3365
3366 // parameters and data for [Cra] algorithm
3367 const cln::cl_N lambda = cln::cl_N("319/320");
3368 int L1;
3369 int L2;
3370 std::vector<std::vector<cln::cl_N> > f_kj;
3371 std::vector<cln::cl_N> crB;
3372 std::vector<std::vector<cln::cl_N> > crG;
3373 std::vector<cln::cl_N> crX;
3374
3375
3376 void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
3377 {
3378         const int size = a.size();
3379         for (int n=0; n<size; n++) {
3380                 c[n] = 0;
3381                 for (int m=0; m<=n; m++) {
3382                         c[n] = c[n] + a[m]*b[n-m];
3383                 }
3384         }
3385 }
3386
3387
3388 // [Cra] section 4
3389 void initcX(const std::vector<int>& s)
3390 {
3391         const int k = s.size();
3392
3393         crX.clear();
3394         crG.clear();
3395         crB.clear();
3396
3397         for (int i=0; i<=L2; i++) {
3398                 crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
3399         }
3400
3401         int Sm = 0;
3402         int Smp1 = 0;
3403         for (int m=0; m<k-1; m++) {
3404                 std::vector<cln::cl_N> crGbuf;
3405                 Sm = Sm + s[m];
3406                 Smp1 = Sm + s[m+1];
3407                 for (int i=0; i<=L2; i++) {
3408                         crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
3409                 }
3410                 crG.push_back(crGbuf);
3411         }
3412
3413         crX = crB;
3414
3415         for (int m=0; m<k-1; m++) {
3416                 std::vector<cln::cl_N> Xbuf;
3417                 for (int i=0; i<=L2; i++) {
3418                         Xbuf.push_back(crX[i] * crG[m][i]);
3419                 }
3420                 halfcyclic_convolute(Xbuf, crB, crX);
3421         }
3422 }
3423
3424
3425 // [Cra] section 4
3426 cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
3427 {
3428         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
3429         cln::cl_N factor = cln::expt(lambda, Sqk);
3430         cln::cl_N res = factor / Sqk * crX[0] * one;
3431         cln::cl_N resbuf;
3432         int N = 0;
3433         do {
3434                 resbuf = res;
3435                 factor = factor * lambda;
3436                 N++;
3437                 res = res + crX[N] * factor / (N+Sqk);
3438         } while ((res != resbuf) || cln::zerop(crX[N]));
3439         return res;
3440 }
3441
3442
3443 // [Cra] section 4
3444 void calc_f(int maxr)
3445 {
3446         f_kj.clear();
3447         f_kj.resize(L1);
3448
3449         cln::cl_N t0, t1, t2, t3, t4;
3450         int i, j, k;
3451         std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
3452         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
3453
3454         t0 = cln::exp(-lambda);
3455         t2 = 1;
3456         for (k=1; k<=L1; k++) {
3457                 t1 = k * lambda;
3458                 t2 = t0 * t2;
3459                 for (j=1; j<=maxr; j++) {
3460                         t3 = 1;
3461                         t4 = 1;
3462                         for (i=2; i<=j; i++) {
3463                                 t4 = t4 * (j-i+1);
3464                                 t3 = t1 * t3 + t4;
3465                         }
3466                         (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one);
3467                 }
3468                 it++;
3469         }
3470 }
3471
3472
3473 // [Cra] (3.1)
3474 cln::cl_N crandall_Z(const std::vector<int>& s)
3475 {
3476         const int j = s.size();
3477
3478         if (j == 1) {
3479                 cln::cl_N t0;
3480                 cln::cl_N t0buf;
3481                 int q = 0;
3482                 do {
3483                         t0buf = t0;
3484                         q++;
3485                         t0 = t0 + f_kj[q+j-2][s[0]-1];
3486                 } while (t0 != t0buf);
3487
3488                 return t0 / cln::factorial(s[0]-1);
3489         }
3490
3491         std::vector<cln::cl_N> t(j);
3492
3493         cln::cl_N t0buf;
3494         int q = 0;
3495         do {
3496                 t0buf = t[0];
3497                 q++;
3498                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
3499                 for (int k=j-2; k>=1; k--) {
3500                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
3501                 }
3502                 t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
3503         } while (t[0] != t0buf);
3504
3505         return t[0] / cln::factorial(s[0]-1);
3506 }
3507
3508
3509 // [Cra] (2.4)
3510 cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
3511 {
3512         std::vector<int> r = s;
3513         const int j = r.size();
3514
3515         // decide on maximal size of f_kj for crandall_Z
3516         if (Digits < 50) {
3517                 L1 = 150;
3518         } else {
3519                 L1 = Digits * 3 + j*2;
3520         }
3521
3522         // decide on maximal size of crX for crandall_Y
3523         if (Digits < 38) {
3524                 L2 = 63;
3525         } else if (Digits < 86) {
3526                 L2 = 127;
3527         } else if (Digits < 192) {
3528                 L2 = 255;
3529         } else if (Digits < 394) {
3530                 L2 = 511;
3531         } else if (Digits < 808) {
3532                 L2 = 1023;
3533         } else {
3534                 L2 = 2047;
3535         }
3536
3537         cln::cl_N res;
3538
3539         int maxr = 0;
3540         int S = 0;
3541         for (int i=0; i<j; i++) {
3542                 S += r[i];
3543                 if (r[i] > maxr) {
3544                         maxr = r[i];
3545                 }
3546         }
3547
3548         calc_f(maxr);
3549
3550         const cln::cl_N r0factorial = cln::factorial(r[0]-1);
3551
3552         std::vector<int> rz;
3553         int skp1buf;
3554         int Srun = S;
3555         for (int k=r.size()-1; k>0; k--) {
3556
3557                 rz.insert(rz.begin(), r.back());
3558                 skp1buf = rz.front();
3559                 Srun -= skp1buf;
3560                 r.pop_back();
3561
3562                 initcX(r);
3563
3564                 for (int q=0; q<skp1buf; q++) {
3565
3566                         cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
3567                         cln::cl_N pp2 = crandall_Z(rz);
3568
3569                         rz.front()--;
3570
3571                         if (q & 1) {
3572                                 res = res - pp1 * pp2 / cln::factorial(q);
3573                         } else {
3574                                 res = res + pp1 * pp2 / cln::factorial(q);
3575                         }
3576                 }
3577                 rz.front() = skp1buf;
3578         }
3579         rz.insert(rz.begin(), r.back());
3580
3581         initcX(rz);
3582
3583         res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
3584
3585         return res;
3586 }
3587
3588
3589 cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
3590 {
3591         const int j = r.size();
3592
3593         // buffer for subsums
3594         std::vector<cln::cl_N> t(j);
3595         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
3596
3597         cln::cl_N t0buf;
3598         int q = 0;
3599         do {
3600                 t0buf = t[0];
3601                 q++;
3602                 t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
3603                 for (int k=j-2; k>=0; k--) {
3604                         t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
3605                 }
3606         } while (t[0] != t0buf);
3607
3608         return t[0];
3609 }
3610
3611
3612 // does Hoelder convolution. see [BBB] (7.0)
3613 cln::cl_N zeta_do_Hoelder_convolution(const std::vector<int>& m_, const std::vector<int>& s_)
3614 {
3615         // prepare parameters
3616         // holds Li arguments in [BBB] notation
3617         std::vector<int> s = s_;
3618         std::vector<int> m_p = m_;
3619         std::vector<int> m_q;
3620         // holds Li arguments in nested sums notation
3621         std::vector<cln::cl_N> s_p(s.size(), cln::cl_N(1));
3622         s_p[0] = s_p[0] * cln::cl_N("1/2");
3623         // convert notations
3624         int sig = 1;
3625         for (int i=0; i<s_.size(); i++) {
3626                 if (s_[i] < 0) {
3627                         sig = -sig;
3628                         s_p[i] = -s_p[i];
3629                 }
3630                 s[i] = sig * std::abs(s[i]);
3631         }
3632         std::vector<cln::cl_N> s_q;
3633         cln::cl_N signum = 1;
3634
3635         // first term
3636         cln::cl_N res = multipleLi_do_sum(m_p, s_p);
3637
3638         // middle terms
3639         do {
3640
3641                 // change parameters
3642                 if (s.front() > 0) {
3643                         if (m_p.front() == 1) {
3644                                 m_p.erase(m_p.begin());
3645                                 s_p.erase(s_p.begin());
3646                                 if (s_p.size() > 0) {
3647                                         s_p.front() = s_p.front() * cln::cl_N("1/2");
3648                                 }
3649                                 s.erase(s.begin());
3650                                 m_q.front()++;
3651                         } else {
3652                                 m_p.front()--;
3653                                 m_q.insert(m_q.begin(), 1);
3654                                 if (s_q.size() > 0) {
3655                                         s_q.front() = s_q.front() * 2;
3656                                 }
3657                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
3658                         }
3659                 } else {
3660                         if (m_p.front() == 1) {
3661                                 m_p.erase(m_p.begin());
3662                                 cln::cl_N spbuf = s_p.front();
3663                                 s_p.erase(s_p.begin());
3664                                 if (s_p.size() > 0) {
3665                                         s_p.front() = s_p.front() * spbuf;
3666                                 }
3667                                 s.erase(s.begin());
3668                                 m_q.insert(m_q.begin(), 1);
3669                                 if (s_q.size() > 0) {
3670                                         s_q.front() = s_q.front() * 4;
3671                                 }
3672                                 s_q.insert(s_q.begin(), cln::cl_N("1/4"));
3673                                 signum = -signum;
3674                         } else {
3675                                 m_p.front()--;
3676                                 m_q.insert(m_q.begin(), 1);
3677                                 if (s_q.size() > 0) {
3678                                         s_q.front() = s_q.front() * 2;
3679                                 }
3680                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
3681                         }
3682                 }
3683
3684                 // exiting the loop
3685                 if (m_p.size() == 0) break;
3686
3687                 res = res + signum * multipleLi_do_sum(m_p, s_p) * multipleLi_do_sum(m_q, s_q);
3688
3689         } while (true);
3690
3691         // last term
3692         res = res + signum * multipleLi_do_sum(m_q, s_q);
3693
3694         return res;
3695 }
3696
3697
3698 } // end of anonymous namespace
3699
3700
3701 //////////////////////////////////////////////////////////////////////
3702 //
3703 // Multiple zeta values  zeta(x)
3704 //
3705 // GiNaC function
3706 //
3707 //////////////////////////////////////////////////////////////////////
3708
3709
3710 static ex zeta1_evalf(const ex& x)
3711 {
3712         if (is_exactly_a<lst>(x) && (x.nops()>1)) {
3713
3714                 // multiple zeta value
3715                 const int count = x.nops();
3716                 const lst& xlst = ex_to<lst>(x);
3717                 std::vector<int> r(count);
3718
3719                 // check parameters and convert them
3720                 lst::const_iterator it1 = xlst.begin();
3721                 std::vector<int>::iterator it2 = r.begin();
3722                 do {
3723                         if (!(*it1).info(info_flags::posint)) {
3724                                 return zeta(x).hold();
3725                         }
3726                         *it2 = ex_to<numeric>(*it1).to_int();
3727                         it1++;
3728                         it2++;
3729                 } while (it2 != r.end());
3730
3731                 // check for divergence
3732                 if (r[0] == 1) {
3733                         return zeta(x).hold();
3734                 }
3735
3736                 // decide on summation algorithm
3737                 // this is still a bit clumsy
3738                 int limit = (Digits>17) ? 10 : 6;
3739                 if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) {
3740                         return numeric(zeta_do_sum_Crandall(r));
3741                 } else {
3742                         return numeric(zeta_do_sum_simple(r));
3743                 }
3744         }
3745
3746         // single zeta value
3747         if (is_exactly_a<numeric>(x) && (x != 1)) {
3748                 try {
3749                         return zeta(ex_to<numeric>(x));
3750                 } catch (const dunno &e) { }
3751         }
3752
3753         return zeta(x).hold();
3754 }
3755
3756
3757 static ex zeta1_eval(const ex& m)
3758 {
3759         if (is_exactly_a<lst>(m)) {
3760                 if (m.nops() == 1) {
3761                         return zeta(m.op(0));
3762                 }
3763                 return zeta(m).hold();
3764         }
3765
3766         if (m.info(info_flags::numeric)) {
3767                 const numeric& y = ex_to<numeric>(m);
3768                 // trap integer arguments:
3769                 if (y.is_integer()) {
3770                         if (y.is_zero()) {
3771                                 return _ex_1_2;
3772                         }
3773                         if (y.is_equal(*_num1_p)) {
3774                                 return zeta(m).hold();
3775                         }
3776                         if (y.info(info_flags::posint)) {
3777                                 if (y.info(info_flags::odd)) {
3778                                         return zeta(m).hold();
3779                                 } else {
3780                                         return abs(bernoulli(y)) * pow(Pi, y) * pow(*_num2_p, y-(*_num1_p)) / factorial(y);
3781                                 }
3782                         } else {
3783                                 if (y.info(info_flags::odd)) {
3784                                         return -bernoulli((*_num1_p)-y) / ((*_num1_p)-y);
3785                                 } else {
3786                                         return _ex0;
3787                                 }
3788                         }
3789                 }
3790                 // zeta(float)
3791                 if (y.info(info_flags::numeric) && !y.info(info_flags::crational)) {
3792                         return zeta1_evalf(m);
3793                 }
3794         }
3795         return zeta(m).hold();
3796 }
3797
3798
3799 static ex zeta1_deriv(const ex& m, unsigned deriv_param)
3800 {
3801         GINAC_ASSERT(deriv_param==0);
3802
3803         if (is_exactly_a<lst>(m)) {
3804                 return _ex0;
3805         } else {
3806                 return zetaderiv(_ex1, m);
3807         }
3808 }
3809
3810
3811 static void zeta1_print_latex(const ex& m_, const print_context& c)
3812 {
3813         c.s << "\\zeta(";
3814         if (is_a<lst>(m_)) {
3815                 const lst& m = ex_to<lst>(m_);
3816                 lst::const_iterator it = m.begin();
3817                 (*it).print(c);
3818                 it++;
3819                 for (; it != m.end(); it++) {
3820                         c.s << ",";
3821                         (*it).print(c);
3822                 }
3823         } else {
3824                 m_.print(c);
3825         }
3826         c.s << ")";
3827 }
3828
3829
3830 unsigned zeta1_SERIAL::serial = function::register_new(function_options("zeta", 1).
3831                                 evalf_func(zeta1_evalf).
3832                                 eval_func(zeta1_eval).
3833                                 derivative_func(zeta1_deriv).
3834                                 print_func<print_latex>(zeta1_print_latex).
3835                                 do_not_evalf_params().
3836                                 overloaded(2));
3837
3838
3839 //////////////////////////////////////////////////////////////////////
3840 //
3841 // Alternating Euler sum  zeta(x,s)
3842 //
3843 // GiNaC function
3844 //
3845 //////////////////////////////////////////////////////////////////////
3846
3847
3848 static ex zeta2_evalf(const ex& x, const ex& s)
3849 {
3850         if (is_exactly_a<lst>(x)) {
3851
3852                 // alternating Euler sum
3853                 const int count = x.nops();
3854                 const lst& xlst = ex_to<lst>(x);
3855                 const lst& slst = ex_to<lst>(s);
3856                 std::vector<int> xi(count);
3857                 std::vector<int> si(count);
3858
3859                 // check parameters and convert them
3860                 lst::const_iterator it_xread = xlst.begin();
3861                 lst::const_iterator it_sread = slst.begin();
3862                 std::vector<int>::iterator it_xwrite = xi.begin();
3863                 std::vector<int>::iterator it_swrite = si.begin();
3864                 do {
3865                         if (!(*it_xread).info(info_flags::posint)) {
3866                                 return zeta(x, s).hold();
3867                         }
3868                         *it_xwrite = ex_to<numeric>(*it_xread).to_int();
3869                         if (*it_sread > 0) {
3870                                 *it_swrite = 1;
3871                         } else {
3872                                 *it_swrite = -1;
3873                         }
3874                         it_xread++;
3875                         it_sread++;
3876                         it_xwrite++;
3877                         it_swrite++;
3878                 } while (it_xwrite != xi.end());
3879
3880                 // check for divergence
3881                 if ((xi[0] == 1) && (si[0] == 1)) {
3882                         return zeta(x, s).hold();
3883                 }
3884
3885                 // use Hoelder convolution
3886                 return numeric(zeta_do_Hoelder_convolution(xi, si));
3887         }
3888
3889         return zeta(x, s).hold();
3890 }
3891
3892
3893 static ex zeta2_eval(const ex& m, const ex& s_)
3894 {
3895         if (is_exactly_a<lst>(s_)) {
3896                 const lst& s = ex_to<lst>(s_);
3897                 for (lst::const_iterator it = s.begin(); it != s.end(); it++) {
3898                         if ((*it).info(info_flags::positive)) {
3899                                 continue;
3900                         }
3901                         return zeta(m, s_).hold();
3902                 }
3903                 return zeta(m);
3904         } else if (s_.info(info_flags::positive)) {
3905                 return zeta(m);
3906         }
3907
3908         return zeta(m, s_).hold();
3909 }
3910
3911
3912 static ex zeta2_deriv(const ex& m, const ex& s, unsigned deriv_param)
3913 {
3914         GINAC_ASSERT(deriv_param==0);
3915
3916         if (is_exactly_a<lst>(m)) {
3917                 return _ex0;
3918         } else {
3919                 if ((is_exactly_a<lst>(s) && s.op(0).info(info_flags::positive)) || s.info(info_flags::positive)) {
3920                         return zetaderiv(_ex1, m);
3921                 }
3922                 return _ex0;
3923         }
3924 }
3925
3926
3927 static void zeta2_print_latex(const ex& m_, const ex& s_, const print_context& c)
3928 {
3929         lst m;
3930         if (is_a<lst>(m_)) {
3931                 m = ex_to<lst>(m_);
3932         } else {
3933                 m = lst(m_);
3934         }
3935         lst s;
3936         if (is_a<lst>(s_)) {
3937                 s = ex_to<lst>(s_);
3938         } else {
3939                 s = lst(s_);
3940         }
3941         c.s << "\\zeta(";
3942         lst::const_iterator itm = m.begin();
3943         lst::const_iterator its = s.begin();
3944         if (*its < 0) {
3945                 c.s << "\\overline{";
3946                 (*itm).print(c);
3947                 c.s << "}";
3948         } else {
3949                 (*itm).print(c);
3950         }
3951         its++;
3952         itm++;
3953         for (; itm != m.end(); itm++, its++) {
3954                 c.s << ",";
3955                 if (*its < 0) {
3956                         c.s << "\\overline{";
3957                         (*itm).print(c);
3958                         c.s << "}";
3959                 } else {
3960                         (*itm).print(c);
3961                 }
3962         }
3963         c.s << ")";
3964 }
3965
3966
3967 unsigned zeta2_SERIAL::serial = function::register_new(function_options("zeta", 2).
3968                                 evalf_func(zeta2_evalf).
3969                                 eval_func(zeta2_eval).
3970                                 derivative_func(zeta2_deriv).
3971                                 print_func<print_latex>(zeta2_print_latex).
3972                                 do_not_evalf_params().
3973                                 overloaded(2));
3974
3975
3976 } // namespace GiNaC
3977