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