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  *    nielsen's generalized polylogarithm  S(n,p,x)
   9  *    harmonic polylogarithm               H(m,x) or H(lst(m_1,...,m_k),x)
  10  *    multiple zeta value                  zeta(m) or zeta(lst(m_1,...,m_k))
  11  *    alternating Euler sum                zeta(m,s) or zeta(lst(m_1,...,m_k),lst(s_1,...,s_k))
  12  *
  13  *  Some remarks:
  14  *
  15  *    - All formulae used can be looked up in the following publications:
  16  *      [Kol] Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
  17  *      [Cra] Fast Evaluation of Multiple Zeta Sums, R.E.Crandall, Math.Comp. 67 (1998), pp. 1163-1172.
  18  *      [ReV] Harmonic Polylogarithms, E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
  19  *      [BBB] Special Values of Multiple Polylogarithms, J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001), pp. 907-941
  20  *
  21  *    - The order of parameters and arguments of Li and zeta is defined according to the nested sums
  22  *      representation. The parameters for H are understood as in [ReV]. They can be in expanded --- only
  23  *      0, 1 and -1 --- or in compactified --- a string with zeros in front of 1 or -1 is written as a single
  24  *      number --- notation.
  25  *
  26  *    - Except for the multiple polylogarithm all functions can be nummerically evaluated with arguments in
  27  *      the whole complex plane. Multiple polylogarithms evaluate only if for each argument x_i the product
  28  *      x_1 * x_2 * ... * x_i is smaller than one. The parameters for Li, zeta and S must be positive integers.
  29  *      If you want to have an alternating Euler sum, you have to give the signs of the parameters as a
  30  *      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.
  35  *
  36  *    - The functions have no series expansion into nested sums. To do this, you have to convert these functions
  37  *      into the appropriate objects from the nestedsums library, do the expansion and convert the
  38  *      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.
  45  *
  46  */
  47
  48 /*
  49  *  GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
  50  *
  51  *  This program is free software; you can redistribute it and/or modify
  52  *  it under the terms of the GNU General Public License as published by
  53  *  the Free Software Foundation; either version 2 of the License, or
  54  *  (at your option) any later version.
  55  *
  56  *  This program is distributed in the hope that it will be useful,
  57  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  58  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  59  *  GNU General Public License for more details.
  60  *
  61  *  You should have received a copy of the GNU General Public License
  62  *  along with this program; if not, write to the Free Software
  63  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  64  */
  65
  66 #include <stdexcept>
  67 #include <vector>
  68 #include <cln/cln.h>
  69
  70 #include "inifcns.h"
  71
  72 #include "add.h"
  73 #include "constant.h"
  74 #include "lst.h"
  75 #include "mul.h"
  76 #include "numeric.h"
  77 #include "operators.h"
  78 #include "power.h"
  79 #include "pseries.h"
  80 #include "relational.h"
  81 #include "symbol.h"
  82 #include "utils.h"
  83 #include "wildcard.h"
  84
  85
  86 namespace GiNaC {
  87
  88
  89 //////////////////////////////////////////////////////////////////////
  90 //
  91 // Classical polylogarithm  Li(n,x)
  92 //
  93 // helper functions
  94 //
  95 //////////////////////////////////////////////////////////////////////
  96
  97
  98 // anonymous namespace for helper functions
  99 namespace {
 100
 101
 102 // lookup table for factors built from Bernoulli numbers
 103 // see fill_Xn()
 104 std::vector<std::vector<cln::cl_N> > Xn;
 105 // initial size of Xn that should suffice for 32bit machines (must be even)
 106 const int xninitsizestep = 26;
 107 int xninitsize = xninitsizestep;
 108 int xnsize = 0;
 109
 110
 111 // This function calculates the X_n. The X_n are needed for speed up of classical polylogarithms.
 112 // With these numbers the polylogs can be calculated as follows:
 113 //   Li_p (x)  =  \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with  u = -log(1-x)
 114 //   X_0(n) = B_n (Bernoulli numbers)
 115 //   X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
 116 // The calculation of Xn depends on X0 and X{n-1}.
 117 // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
 118 // This results in a slightly more complicated algorithm for the X_n.
 119 // The first index in Xn corresponds to the index of the polylog minus 2.
 120 // The second index in Xn corresponds to the index from the actual sum.
 121 void fill_Xn(int n)
 122 {
 123         if (n>1) {
 124                 // calculate X_2 and higher (corresponding to Li_4 and higher)
 125                 std::vector<cln::cl_N> buf(xninitsize);
 126                 std::vector<cln::cl_N>::iterator it = buf.begin();
 127                 cln::cl_N result;
 128                 *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
 129                 it++;
 130                 for (int i=2; i<=xninitsize; i++) {
 131                         if (i&1) {
 132                                 result = 0; // k == 0
 133                         } else {
 134                                 result = Xn[0][i/2-1]; // k == 0
 135                         }
 136                         for (int k=1; k<i-1; k++) {
 137                                 if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
 138                                         result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
 139                                 }
 140                         }
 141                         result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
 142                         result = result + Xn[n-1][i-1] / (i+1); // k == i
 143
 144                         *it = result;
 145                         it++;
 146                 }
 147                 Xn.push_back(buf);
 148         } else if (n==1) {
 149                 // special case to handle the X_0 correct
 150                 std::vector<cln::cl_N> buf(xninitsize);
 151                 std::vector<cln::cl_N>::iterator it = buf.begin();
 152                 cln::cl_N result;
 153                 *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
 154                 it++;
 155                 *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
 156                 it++;
 157                 for (int i=3; i<=xninitsize; i++) {
 158                         if (i & 1) {
 159                                 result = -Xn[0][(i-3)/2]/2;
 160                                 *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
 161                                 it++;
 162                         } else {
 163                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
 164                                 for (int k=1; k<i/2; k++) {
 165                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
 166                                 }
 167                                 *it = result;
 168                                 it++;
 169                         }
 170                 }
 171                 Xn.push_back(buf);
 172         } else {
 173                 // calculate X_0
 174                 std::vector<cln::cl_N> buf(xninitsize/2);
 175                 std::vector<cln::cl_N>::iterator it = buf.begin();
 176                 for (int i=1; i<=xninitsize/2; i++) {
 177                         *it = bernoulli(i*2).to_cl_N();
 178                         it++;
 179                 }
 180                 Xn.push_back(buf);
 181         }
 182
 183         xnsize++;
 184 }
 185
 186
 187 // doubles the number of entries in each Xn[]
 188 void double_Xn()
 189 {
 190         const int pos0 = xninitsize / 2;
 191         // X_0
 192         for (int i=1; i<=xninitsizestep/2; ++i) {
 193                 Xn[0].push_back(bernoulli((i+pos0)*2).to_cl_N());
 194         }
 195         if (Xn.size() > 0) {
 196                 int xend = xninitsize + xninitsizestep;
 197                 cln::cl_N result;
 198                 // X_1
 199                 for (int i=xninitsize+1; i<=xend; ++i) {
 200                         if (i & 1) {
 201                                 result = -Xn[0][(i-3)/2]/2;
 202                                 Xn[1].push_back((cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result);
 203                         } else {
 204                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
 205                                 for (int k=1; k<i/2; k++) {
 206                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
 207                                 }
 208                                 Xn[1].push_back(result);
 209                         }
 210                 }
 211                 // X_n
 212                 for (int n=2; n<Xn.size(); ++n) {
 213                         for (int i=xninitsize+1; i<=xend; ++i) {
 214                                 if (i & 1) {
 215                                         result = 0; // k == 0
 216                                 } else {
 217                                         result = Xn[0][i/2-1]; // k == 0
 218                                 }
 219                                 for (int k=1; k<i-1; ++k) {
 220                                         if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
 221                                                 result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
 222                                         }
 223                                 }
 224                                 result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
 225                                 result = result + Xn[n-1][i-1] / (i+1); // k == i
 226                                 Xn[n].push_back(result);
 227                         }
 228                 }
 229         }
 230         xninitsize += xninitsizestep;
 231 }
 232
 233
 234 // calculates Li(2,x) without Xn
 235 cln::cl_N Li2_do_sum(const cln::cl_N& x)
 236 {
 237         cln::cl_N res = x;
 238         cln::cl_N resbuf;
 239         cln::cl_N num = x * cln::cl_float(1, cln::float_format(Digits));
 240         cln::cl_I den = 1; // n^2 = 1
 241         unsigned i = 3;
 242         do {
 243                 resbuf = res;
 244                 num = num * x;
 245                 den = den + i;  // n^2 = 4, 9, 16, ...
 246                 i += 2;
 247                 res = res + num / den;
 248         } while (res != resbuf);
 249         return res;
 250 }
 251
 252
 253 // calculates Li(2,x) with Xn
 254 cln::cl_N Li2_do_sum_Xn(const cln::cl_N& x)
 255 {
 256         std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
 257         std::vector<cln::cl_N>::const_iterator xend = Xn[0].end();
 258         cln::cl_N u = -cln::log(1-x);
 259         cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
 260         cln::cl_N uu = cln::square(u);
 261         cln::cl_N res = u - uu/4;
 262         cln::cl_N resbuf;
 263         unsigned i = 1;
 264         do {
 265                 resbuf = res;
 266                 factor = factor * uu / (2*i * (2*i+1));
 267                 res = res + (*it) * factor;
 268                 i++;
 269                 if (++it == xend) {
 270                         double_Xn();
 271                         it = Xn[0].begin() + (i-1);
 272                         xend = Xn[0].end();
 273                 }
 274         } while (res != resbuf);
 275         return res;
 276 }
 277
 278
 279 // calculates Li(n,x), n>2 without Xn
 280 cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
 281 {
 282         cln::cl_N factor = x * cln::cl_float(1, cln::float_format(Digits));
 283         cln::cl_N res = x;
 284         cln::cl_N resbuf;
 285         int i=2;
 286         do {
 287                 resbuf = res;
 288                 factor = factor * x;
 289                 res = res + factor / cln::expt(cln::cl_I(i),n);
 290                 i++;
 291         } while (res != resbuf);
 292         return res;
 293 }
 294
 295
 296 // calculates Li(n,x), n>2 with Xn
 297 cln::cl_N Lin_do_sum_Xn(int n, const cln::cl_N& x)
 298 {
 299         std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
 300         std::vector<cln::cl_N>::const_iterator xend = Xn[n-2].end();
 301         cln::cl_N u = -cln::log(1-x);
 302         cln::cl_N factor = u * cln::cl_float(1, cln::float_format(Digits));
 303         cln::cl_N res = u;
 304         cln::cl_N resbuf;
 305         unsigned i=2;
 306         do {
 307                 resbuf = res;
 308                 factor = factor * u / i;
 309                 res = res + (*it) * factor;
 310                 i++;
 311                 if (++it == xend) {
 312                         double_Xn();
 313                         it = Xn[n-2].begin() + (i-2);
 314                         xend = Xn[n-2].end();
 315                 }
 316         } while (res != resbuf);
 317         return res;
 318 }
 319
 320
 321 // forward declaration needed by function Li_projection and C below
 322 numeric S_num(int n, int p, const numeric& x);
 323
 324
 325 // helper function for classical polylog Li
 326 cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
 327 {
 328         // treat n=2 as special case
 329         if (n == 2) {
 330                 // check if precalculated X0 exists
 331                 if (xnsize == 0) {
 332                         fill_Xn(0);
 333                 }
 334
 335                 if (cln::realpart(x) < 0.5) {
 336                         // choose the faster algorithm
 337                         // the switching point was empirically determined. the optimal point
 338                         // depends on hardware, Digits, ... so an approx value is okay.
 339                         // it solves also the problem with precision due to the u=-log(1-x) transformation
 340                         if (cln::abs(cln::realpart(x)) < 0.25) {
 341
 342                                 return Li2_do_sum(x);
 343                         } else {
 344                                 return Li2_do_sum_Xn(x);
 345                         }
 346                 } else {
 347                         // choose the faster algorithm
 348                         if (cln::abs(cln::realpart(x)) > 0.75) {
 349                                 return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
 350                         } else {
 351                                 return -Li2_do_sum_Xn(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
 352                         }
 353                 }
 354         } else {
 355                 // check if precalculated Xn exist
 356                 if (n > xnsize+1) {
 357                         for (int i=xnsize; i<n-1; i++) {
 358                                 fill_Xn(i);
 359                         }
 360                 }
 361
 362                 if (cln::realpart(x) < 0.5) {
 363                         // choose the faster algorithm
 364                         // with n>=12 the "normal" summation always wins against the method with Xn
 365                         if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
 366                                 return Lin_do_sum(n, x);
 367                         } else {
 368                                 return Lin_do_sum_Xn(n, x);
 369                         }
 370                 } else {
 371                         cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
 372                         for (int j=0; j<n-1; j++) {
 373                                 result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
 374                                                   * cln::expt(cln::log(x), j) / cln::factorial(j);
 375                         }
 376                         return result;
 377                 }
 378         }
 379 }
 380
 381
 382 // helper function for classical polylog Li
 383 numeric Li_num(int n, const numeric& x)
 384 {
 385         if (n == 1) {
 386                 // just a log
 387                 return -cln::log(1-x.to_cl_N());
 388         }
 389         if (x.is_zero()) {
 390                 return 0;
 391         }
 392         if (x == 1) {
 393                 // [Kol] (2.22)
 394                 return cln::zeta(n);
 395         }
 396         else if (x == -1) {
 397                 // [Kol] (2.22)
 398                 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
 399         }
 400
 401         // what is the desired float format?
 402         // first guess: default format
 403         cln::float_format_t prec = cln::default_float_format;
 404         const cln::cl_N value = x.to_cl_N();
 405         // second guess: the argument's format
 406         if (!x.real().is_rational())
 407                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
 408         else if (!x.imag().is_rational())
 409                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 410
 411         // [Kol] (5.15)
 412         if (cln::abs(value) > 1) {
 413                 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
 414                 // check if argument is complex. if it is real, the new polylog has to be conjugated.
 415                 if (cln::zerop(cln::imagpart(value))) {
 416                         if (n & 1) {
 417                                 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
 418                         }
 419                         else {
 420                                 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
 421                         }
 422                 }
 423                 else {
 424                         if (n & 1) {
 425                                 result = result + Li_projection(n, cln::recip(value), prec);
 426                         }
 427                         else {
 428                                 result = result - Li_projection(n, cln::recip(value), prec);
 429                         }
 430                 }
 431                 cln::cl_N add;
 432                 for (int j=0; j<n-1; j++) {
 433                         add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
 434                                     * Li_num(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
 435                 }
 436                 result = result - add;
 437                 return result;
 438         }
 439         else {
 440                 return Li_projection(n, value, prec);
 441         }
 442 }
 443
 444
 445 } // end of anonymous namespace
 446
 447
 448 //////////////////////////////////////////////////////////////////////
 449 //
 450 // Multiple polylogarithm  Li(n,x)
 451 //
 452 // helper function
 453 //
 454 //////////////////////////////////////////////////////////////////////
 455
 456
 457 // anonymous namespace for helper function
 458 namespace {
 459
 460
 461 cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
 462 {
 463         const int j = s.size();
 464
 465         std::vector<cln::cl_N> t(j);
 466         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
 467
 468         cln::cl_N t0buf;
 469         int q = 0;
 470         do {
 471                 t0buf = t[0];
 472                 // do it once ...
 473                 q++;
 474                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
 475                 for (int k=j-2; k>=0; k--) {
 476                         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]);
 477                 }
 478                 // ... and do it again (to avoid premature drop out due to special arguments)
 479                 q++;
 480                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
 481                 for (int k=j-2; k>=0; k--) {
 482                         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]);
 483                 }
 484         } while (t[0] != t0buf);
 485
 486         return t[0];
 487 }
 488
 489 // forward declaration for Li_eval()
 490 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf);
 491
 492
 493 } // end of anonymous namespace
 494
 495
 496 //////////////////////////////////////////////////////////////////////
 497 //
 498 // Classical polylogarithm and multiple polylogarithm  Li(n,x)
 499 //
 500 // GiNaC function
 501 //
 502 //////////////////////////////////////////////////////////////////////
 503
 504
 505 static ex Li_evalf(const ex& x1, const ex& x2)
 506 {
 507         // classical polylogs
 508         if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
 509                 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
 510         }
 511         if (is_a<numeric>(x1) && !is_a<lst>(x2)) {
 512                 // try to numerically evaluate second argument
 513                 ex x2_val = x2.evalf();
 514                 if (is_a<numeric>(x2_val)) {
 515                         return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2_val));
 516                 } else {
 517                         return Li(x1, x2).hold();
 518                 }
 519         }
 520         // multiple polylogs
 521         else if (is_a<lst>(x1) && is_a<lst>(x2)) {
 522                 ex conv = 1;
 523                 for (int i=0; i<x1.nops(); i++) {
 524                         if (!x1.op(i).info(info_flags::posint)) {
 525                                 return Li(x1, x2).hold();
 526                         }
 527                         if (!is_a<numeric>(x2.op(i))) {
 528                                 return Li(x1, x2).hold();
 529                         }
 530                         conv *= x2.op(i);
 531                         if (abs(conv) >= 1) {
 532                                 return Li(x1, x2).hold();
 533                         }
 534                 }
 535
 536                 std::vector<int> m;
 537                 std::vector<cln::cl_N> x;
 538                 for (int i=0; i<ex_to<numeric>(x1.nops()).to_int(); i++) {
 539                         m.push_back(ex_to<numeric>(x1.op(i)).to_int());
 540                         x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
 541                 }
 542
 543                 return numeric(multipleLi_do_sum(m, x));
 544         }
 545
 546         return Li(x1,x2).hold();
 547 }
 548
 549
 550 static ex Li_eval(const ex& m_, const ex& x_)
 551 {
 552         if (m_.nops() < 2) {
 553                 ex m;
 554                 if (is_a<lst>(m_)) {
 555                         m = m_.op(0);
 556                 } else {
 557                         m = m_;
 558                 }
 559                 ex x;
 560                 if (is_a<lst>(x_)) {
 561                         x = x_.op(0);
 562                 } else {
 563                         x = x_;
 564                 }
 565                 if (x == _ex0) {
 566                         return _ex0;
 567                 }
 568                 if (x == _ex1) {
 569                         return zeta(m);
 570                 }
 571                 if (x == _ex_1) {
 572                         return (pow(2,1-m)-1) * zeta(m);
 573                 }
 574                 if (m == _ex1) {
 575                         return -log(1-x);
 576                 }
 577                 if (m.info(info_flags::posint) && x.info(info_flags::numeric) && (!x.info(info_flags::crational))) {
 578                         return Li_num(ex_to<numeric>(m).to_int(), ex_to<numeric>(x));
 579                 }
 580         } else {
 581                 bool ish = true;
 582                 bool iszeta = true;
 583                 bool iszero = false;
 584                 bool doevalf = false;
 585                 bool doevalfveto = true;
 586                 const lst& m = ex_to<lst>(m_);
 587                 const lst& x = ex_to<lst>(x_);
 588                 lst::const_iterator itm = m.begin();
 589                 lst::const_iterator itx = x.begin();
 590                 for (; itm != m.end(); itm++, itx++) {
 591                         if (!(*itm).info(info_flags::posint)) {
 592                                 return Li(m_, x_).hold();
 593                         }
 594                         if ((*itx != _ex1) && (*itx != _ex_1)) {
 595                                 if (itx != x.begin()) {
 596                                         ish = false;
 597                                 }
 598                                 iszeta = false;
 599                         }
 600                         if (*itx == _ex0) {
 601                                 iszero = true;
 602                         }
 603                         if (!(*itx).info(info_flags::numeric)) {
 604                                 doevalfveto = false;
 605                         }
 606                         if (!(*itx).info(info_flags::crational)) {
 607                                 doevalf = true;
 608                         }
 609                 }
 610                 if (iszeta) {
 611                         return zeta(m_, x_);
 612                 }
 613                 if (iszero) {
 614                         return _ex0;
 615                 }
 616                 if (ish) {
 617                         ex pf;
 618                         lst newm = convert_parameter_Li_to_H(m, x, pf);
 619                         return pf * H(newm, x[0]);
 620                 }
 621                 if (doevalfveto && doevalf) {
 622                         return Li(m_, x_).evalf();
 623                 }
 624         }
 625         return Li(m_, x_).hold();
 626 }
 627
 628
 629 static ex Li_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
 630 {
 631         epvector seq;
 632         seq.push_back(expair(Li(m, x), 0));
 633         return pseries(rel, seq);
 634 }
 635
 636
 637 static ex Li_deriv(const ex& m_, const ex& x_, unsigned deriv_param)
 638 {
 639         GINAC_ASSERT(deriv_param < 2);
 640         if (deriv_param == 0) {
 641                 return _ex0;
 642         }
 643         if (m_.nops() > 1) {
 644                 throw std::runtime_error("don't know how to derivate multiple polylogarithm!");
 645         }
 646         ex m;
 647         if (is_a<lst>(m_)) {
 648                 m = m_.op(0);
 649         } else {
 650                 m = m_;
 651         }
 652         ex x;
 653         if (is_a<lst>(x_)) {
 654                 x = x_.op(0);
 655         } else {
 656                 x = x_;
 657         }
 658         if (m > 0) {
 659                 return Li(m-1, x) / x;
 660         } else {
 661                 return 1/(1-x);
 662         }
 663 }
 664
 665
 666 static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c)
 667 {
 668         lst m;
 669         if (is_a<lst>(m_)) {
 670                 m = ex_to<lst>(m_);
 671         } else {
 672                 m = lst(m_);
 673         }
 674         lst x;
 675         if (is_a<lst>(x_)) {
 676                 x = ex_to<lst>(x_);
 677         } else {
 678                 x = lst(x_);
 679         }
 680         c.s << "\\mbox{Li}_{";
 681         lst::const_iterator itm = m.begin();
 682         (*itm).print(c);
 683         itm++;
 684         for (; itm != m.end(); itm++) {
 685                 c.s << ",";
 686                 (*itm).print(c);
 687         }
 688         c.s << "}(";
 689         lst::const_iterator itx = x.begin();
 690         (*itx).print(c);
 691         itx++;
 692         for (; itx != x.end(); itx++) {
 693                 c.s << ",";
 694                 (*itx).print(c);
 695         }
 696         c.s << ")";
 697 }
 698
 699
 700 REGISTER_FUNCTION(Li,
 701                   evalf_func(Li_evalf).
 702                   eval_func(Li_eval).
 703                   series_func(Li_series).
 704                   derivative_func(Li_deriv).
 705                   print_func<print_latex>(Li_print_latex).
 706                   do_not_evalf_params());
 707
 708
 709 //////////////////////////////////////////////////////////////////////
 710 //
 711 // Nielsen's generalized polylogarithm  S(n,p,x)
 712 //
 713 // helper functions
 714 //
 715 //////////////////////////////////////////////////////////////////////
 716
 717
 718 // anonymous namespace for helper functions
 719 namespace {
 720
 721
 722 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
 723 // see fill_Yn()
 724 std::vector<std::vector<cln::cl_N> > Yn;
 725 int ynsize = 0; // number of Yn[]
 726 int ynlength = 100; // initial length of all Yn[i]
 727
 728
 729 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
 730 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
 731 // representing S_{n,p}(x).
 732 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
 733 // equivalent Z-sum.
 734 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
 735 // representing S_{n,p}(x).
 736 // The calculation of Y_n uses the values from Y_{n-1}.
 737 void fill_Yn(int n, const cln::float_format_t& prec)
 738 {
 739         const int initsize = ynlength;
 740         //const int initsize = initsize_Yn;
 741         cln::cl_N one = cln::cl_float(1, prec);
 742
 743         if (n) {
 744                 std::vector<cln::cl_N> buf(initsize);
 745                 std::vector<cln::cl_N>::iterator it = buf.begin();
 746                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
 747                 *it = (*itprev) / cln::cl_N(n+1) * one;
 748                 it++;
 749                 itprev++;
 750                 // sums with an index smaller than the depth are zero and need not to be calculated.
 751                 // calculation starts with depth, which is n+2)
 752                 for (int i=n+2; i<=initsize+n; i++) {
 753                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
 754                         it++;
 755                         itprev++;
 756                 }
 757                 Yn.push_back(buf);
 758         } else {
 759                 std::vector<cln::cl_N> buf(initsize);
 760                 std::vector<cln::cl_N>::iterator it = buf.begin();
 761                 *it = 1 * one;
 762                 it++;
 763                 for (int i=2; i<=initsize; i++) {
 764                         *it = *(it-1) + 1 / cln::cl_N(i) * one;
 765                         it++;
 766                 }
 767                 Yn.push_back(buf);
 768         }
 769         ynsize++;
 770 }
 771
 772
 773 // make Yn longer ...
 774 void make_Yn_longer(int newsize, const cln::float_format_t& prec)
 775 {
 776
 777         cln::cl_N one = cln::cl_float(1, prec);
 778
 779         Yn[0].resize(newsize);
 780         std::vector<cln::cl_N>::iterator it = Yn[0].begin();
 781         it += ynlength;
 782         for (int i=ynlength+1; i<=newsize; i++) {
 783                 *it = *(it-1) + 1 / cln::cl_N(i) * one;
 784                 it++;
 785         }
 786
 787         for (int n=1; n<ynsize; n++) {
 788                 Yn[n].resize(newsize);
 789                 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
 790                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
 791                 it += ynlength;
 792                 itprev += ynlength;
 793                 for (int i=ynlength+n+1; i<=newsize+n; i++) {
 794                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
 795                         it++;
 796                         itprev++;
 797                 }
 798         }
 799
 800         ynlength = newsize;
 801 }
 802
 803
 804 // helper function for S(n,p,x)
 805 // [Kol] (7.2)
 806 cln::cl_N C(int n, int p)
 807 {
 808         cln::cl_N result;
 809
 810         for (int k=0; k<p; k++) {
 811                 for (int j=0; j<=(n+k-1)/2; j++) {
 812                         if (k == 0) {
 813                                 if (n & 1) {
 814                                         if (j & 1) {
 815                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 816                                         }
 817                                         else {
 818                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 819                                         }
 820                                 }
 821                         }
 822                         else {
 823                                 if (k & 1) {
 824                                         if (j & 1) {
 825                                                 result = result + cln::factorial(n+k-1)
 826                                                                   * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 827                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 828                                         }
 829                                         else {
 830                                                 result = result - cln::factorial(n+k-1)
 831                                                                   * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 832                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 833                                         }
 834                                 }
 835                                 else {
 836                                         if (j & 1) {
 837                                                 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()
 838                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 839                                         }
 840                                         else {
 841                                                 result = result + cln::factorial(n+k-1)
 842                                                                   * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 843                                                                   / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 844                                         }
 845                                 }
 846                         }
 847                 }
 848         }
 849         int np = n+p;
 850         if ((np-1) & 1) {
 851                 if (((np)/2+n) & 1) {
 852                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 853                 }
 854                 else {
 855                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 856                 }
 857         }
 858
 859         return result;
 860 }
 861
 862
 863 // helper function for S(n,p,x)
 864 // [Kol] remark to (9.1)
 865 cln::cl_N a_k(int k)
 866 {
 867         cln::cl_N result;
 868
 869         if (k == 0) {
 870                 return 1;
 871         }
 872
 873         result = result;
 874         for (int m=2; m<=k; m++) {
 875                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
 876         }
 877
 878         return -result / k;
 879 }
 880
 881
 882 // helper function for S(n,p,x)
 883 // [Kol] remark to (9.1)
 884 cln::cl_N b_k(int k)
 885 {
 886         cln::cl_N result;
 887
 888         if (k == 0) {
 889                 return 1;
 890         }
 891
 892         result = result;
 893         for (int m=2; m<=k; m++) {
 894                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
 895         }
 896
 897         return result / k;
 898 }
 899
 900
 901 // helper function for S(n,p,x)
 902 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 903 {
 904         if (p==1) {
 905                 return Li_projection(n+1, x, prec);
 906         }
 907
 908         // check if precalculated values are sufficient
 909         if (p > ynsize+1) {
 910                 for (int i=ynsize; i<p-1; i++) {
 911                         fill_Yn(i, prec);
 912                 }
 913         }
 914
 915         // should be done otherwise
 916         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
 917         cln::cl_N xf = x * one;
 918         //cln::cl_N xf = x * cln::cl_float(1, prec);
 919
 920         cln::cl_N res;
 921         cln::cl_N resbuf;
 922         cln::cl_N factor = cln::expt(xf, p);
 923         int i = p;
 924         do {
 925                 resbuf = res;
 926                 if (i-p >= ynlength) {
 927                         // make Yn longer
 928                         make_Yn_longer(ynlength*2, prec);
 929                 }
 930                 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? ...
 931                 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
 932                 factor = factor * xf;
 933                 i++;
 934         } while (res != resbuf);
 935
 936         return res;
 937 }
 938
 939
 940 // helper function for S(n,p,x)
 941 cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 942 {
 943         // [Kol] (5.3)
 944         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
 945
 946                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
 947                                    * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
 948
 949                 for (int s=0; s<n; s++) {
 950                         cln::cl_N res2;
 951                         for (int r=0; r<p; r++) {
 952                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
 953                                               * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
 954                         }
 955                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
 956                 }
 957
 958                 return result;
 959         }
 960
 961         return S_do_sum(n, p, x, prec);
 962 }
 963
 964
 965 // helper function for S(n,p,x)
 966 numeric S_num(int n, int p, const numeric& x)
 967 {
 968         if (x == 1) {
 969                 if (n == 1) {
 970                     // [Kol] (2.22) with (2.21)
 971                         return cln::zeta(p+1);
 972                 }
 973
 974                 if (p == 1) {
 975                     // [Kol] (2.22)
 976                         return cln::zeta(n+1);
 977                 }
 978
 979                 // [Kol] (9.1)
 980                 cln::cl_N result;
 981                 for (int nu=0; nu<n; nu++) {
 982                         for (int rho=0; rho<=p; rho++) {
 983                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
 984                                                   * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
 985                         }
 986                 }
 987                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
 988
 989                 return result;
 990         }
 991         else if (x == -1) {
 992                 // [Kol] (2.22)
 993                 if (p == 1) {
 994                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
 995                 }
 996 //              throw std::runtime_error("don't know how to evaluate this function!");
 997         }
 998
 999         // what is the desired float format?
1000         // first guess: default format
1001         cln::float_format_t prec = cln::default_float_format;
1002         const cln::cl_N value = x.to_cl_N();
1003         // second guess: the argument's format
1004         if (!x.real().is_rational())
1005                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
1006         else if (!x.imag().is_rational())
1007                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
1008
1009         // [Kol] (5.3)
1010         if ((cln::realpart(value) < -0.5) || (n == 0)) {
1011
1012                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
1013                                    * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
1014
1015                 for (int s=0; s<n; s++) {
1016                         cln::cl_N res2;
1017                         for (int r=0; r<p; r++) {
1018                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
1019                                               * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
1020                         }
1021                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
1022                 }
1023
1024                 return result;
1025
1026         }
1027         // [Kol] (5.12)
1028         if (cln::abs(value) > 1) {
1029
1030                 cln::cl_N result;
1031
1032                 for (int s=0; s<p; s++) {
1033                         for (int r=0; r<=s; r++) {
1034                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
1035                                                   / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
1036                                                   * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
1037                         }
1038                 }
1039                 result = result * cln::expt(cln::cl_I(-1),n);
1040
1041                 cln::cl_N res2;
1042                 for (int r=0; r<n; r++) {
1043                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
1044                 }
1045                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
1046
1047                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
1048
1049                 return result;
1050         }
1051         else {
1052                 return S_projection(n, p, value, prec);
1053         }
1054 }
1055
1056
1057 } // end of anonymous namespace
1058
1059
1060 //////////////////////////////////////////////////////////////////////
1061 //
1062 // Nielsen's generalized polylogarithm  S(n,p,x)
1063 //
1064 // GiNaC function
1065 //
1066 //////////////////////////////////////////////////////////////////////
1067
1068
1069 static ex S_evalf(const ex& n, const ex& p, const ex& x)
1070 {
1071         if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
1072                 if (is_a<numeric>(x)) {
1073                         return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
1074                 } else {
1075                         ex x_val = x.evalf();
1076                         if (is_a<numeric>(x_val)) {
1077                                 return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x_val));
1078                         }
1079                 }
1080         }
1081         return S(n, p, x).hold();
1082 }
1083
1084
1085 static ex S_eval(const ex& n, const ex& p, const ex& x)
1086 {
1087         if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
1088                 if (x == 0) {
1089                         return _ex0;
1090                 }
1091                 if (x == 1) {
1092                         lst m(n+1);
1093                         for (int i=ex_to<numeric>(p).to_int()-1; i>0; i--) {
1094                                 m.append(1);
1095                         }
1096                         return zeta(m);
1097                 }
1098                 if (p == 1) {
1099                         return Li(n+1, x);
1100                 }
1101                 if (x.info(info_flags::numeric) && (!x.info(info_flags::crational))) {
1102                         return S_num(ex_to<numeric>(n).to_int(), ex_to<numeric>(p).to_int(), ex_to<numeric>(x));
1103                 }
1104         }
1105         if (n.is_zero()) {
1106                 // [Kol] (5.3)
1107                 return pow(-log(1-x), p) / factorial(p);
1108         }
1109         return S(n, p, x).hold();
1110 }
1111
1112
1113 static ex S_series(const ex& n, const ex& p, const ex& x, const relational& rel, int order, unsigned options)
1114 {
1115         epvector seq;
1116         seq.push_back(expair(S(n, p, x), 0));
1117         return pseries(rel, seq);
1118 }
1119
1120
1121 static ex S_deriv(const ex& n, const ex& p, const ex& x, unsigned deriv_param)
1122 {
1123         GINAC_ASSERT(deriv_param < 3);
1124         if (deriv_param < 2) {
1125                 return _ex0;
1126         }
1127         if (n > 0) {
1128                 return S(n-1, p, x) / x;
1129         } else {
1130                 return S(n, p-1, x) / (1-x);
1131         }
1132 }
1133
1134
1135 static void S_print_latex(const ex& n, const ex& p, const ex& x, const print_context& c)
1136 {
1137         c.s << "\\mbox{S}_{";
1138         n.print(c);
1139         c.s << ",";
1140         p.print(c);
1141         c.s << "}(";
1142         x.print(c);
1143         c.s << ")";
1144 }
1145
1146
1147 REGISTER_FUNCTION(S,
1148                   evalf_func(S_evalf).
1149                   eval_func(S_eval).
1150                   series_func(S_series).
1151                   derivative_func(S_deriv).
1152                   print_func<print_latex>(S_print_latex).
1153                   do_not_evalf_params());
1154
1155
1156 //////////////////////////////////////////////////////////////////////
1157 //
1158 // Harmonic polylogarithm  H(m,x)
1159 //
1160 // helper functions
1161 //
1162 //////////////////////////////////////////////////////////////////////
1163
1164
1165 // anonymous namespace for helper functions
1166 namespace {
1167
1168
1169 // regulates the pole (used by 1/x-transformation)
1170 symbol H_polesign("IMSIGN");
1171
1172
1173 // convert parameters from H to Li representation
1174 // parameters are expected to be in expanded form, i.e. only 0, 1 and -1
1175 // returns true if some parameters are negative
1176 bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
1177 {
1178         // expand parameter list
1179         lst mexp;
1180         for (lst::const_iterator it = l.begin(); it != l.end(); it++) {
1181                 if (*it > 1) {
1182                         for (ex count=*it-1; count > 0; count--) {
1183                                 mexp.append(0);
1184                         }
1185                         mexp.append(1);
1186                 } else if (*it < -1) {
1187                         for (ex count=*it+1; count < 0; count++) {
1188                                 mexp.append(0);
1189                         }
1190                         mexp.append(-1);
1191                 } else {
1192                         mexp.append(*it);
1193                 }
1194         }
1195
1196         ex signum = 1;
1197         pf = 1;
1198         bool has_negative_parameters = false;
1199         ex acc = 1;
1200         for (lst::const_iterator it = mexp.begin(); it != mexp.end(); it++) {
1201                 if (*it == 0) {
1202                         acc++;
1203                         continue;
1204                 }
1205                 if (*it > 0) {
1206                         m.append((*it+acc-1) * signum);
1207                 } else {
1208                         m.append((*it-acc+1) * signum);
1209                 }
1210                 acc = 1;
1211                 signum = *it;
1212                 pf *= *it;
1213                 if (pf < 0) {
1214                         has_negative_parameters = true;
1215                 }
1216         }
1217         if (has_negative_parameters) {
1218                 for (int i=0; i<m.nops(); i++) {
1219                         if (m.op(i) < 0) {
1220                                 m.let_op(i) = -m.op(i);
1221                                 s.append(-1);
1222                         } else {
1223                                 s.append(1);
1224                         }
1225                 }
1226         }
1227
1228         return has_negative_parameters;
1229 }
1230
1231
1232 // recursivly transforms H to corresponding multiple polylogarithms
1233 struct map_trafo_H_convert_to_Li : public map_function
1234 {
1235         ex operator()(const ex& e)
1236         {
1237                 if (is_a<add>(e) || is_a<mul>(e)) {
1238                         return e.map(*this);
1239                 }
1240                 if (is_a<function>(e)) {
1241                         std::string name = ex_to<function>(e).get_name();
1242                         if (name == "H") {
1243                                 lst parameter;
1244                                 if (is_a<lst>(e.op(0))) {
1245                                                 parameter = ex_to<lst>(e.op(0));
1246                                 } else {
1247                                         parameter = lst(e.op(0));
1248                                 }
1249                                 ex arg = e.op(1);
1250
1251                                 lst m;
1252                                 lst s;
1253                                 ex pf;
1254                                 if (convert_parameter_H_to_Li(parameter, m, s, pf)) {
1255                                         s.let_op(0) = s.op(0) * arg;
1256                                         return pf * Li(m, s).hold();
1257                                 } else {
1258                                         for (int i=0; i<m.nops(); i++) {
1259                                                 s.append(1);
1260                                         }
1261                                         s.let_op(0) = s.op(0) * arg;
1262                                         return Li(m, s).hold();
1263                                 }
1264                         }
1265                 }
1266                 return e;
1267         }
1268 };
1269
1270
1271 // recursivly transforms H to corresponding zetas
1272 struct map_trafo_H_convert_to_zeta : public map_function
1273 {
1274         ex operator()(const ex& e)
1275         {
1276                 if (is_a<add>(e) || is_a<mul>(e)) {
1277                         return e.map(*this);
1278                 }
1279                 if (is_a<function>(e)) {
1280                         std::string name = ex_to<function>(e).get_name();
1281                         if (name == "H") {
1282                                 lst parameter;
1283                                 if (is_a<lst>(e.op(0))) {
1284                                                 parameter = ex_to<lst>(e.op(0));
1285                                 } else {
1286                                         parameter = lst(e.op(0));
1287                                 }
1288
1289                                 lst m;
1290                                 lst s;
1291                                 ex pf;
1292                                 if (convert_parameter_H_to_Li(parameter, m, s, pf)) {
1293                                         return pf * zeta(m, s);
1294                                 } else {
1295                                         return zeta(m);
1296                                 }
1297                         }
1298                 }
1299                 return e;
1300         }
1301 };
1302
1303
1304 // remove trailing zeros from H-parameters
1305 struct map_trafo_H_reduce_trailing_zeros : public map_function
1306 {
1307         ex operator()(const ex& e)
1308         {
1309                 if (is_a<add>(e) || is_a<mul>(e)) {
1310                         return e.map(*this);
1311                 }
1312                 if (is_a<function>(e)) {
1313                         std::string name = ex_to<function>(e).get_name();
1314                         if (name == "H") {
1315                                 lst parameter;
1316                                 if (is_a<lst>(e.op(0))) {
1317                                         parameter = ex_to<lst>(e.op(0));
1318                                 } else {
1319                                         parameter = lst(e.op(0));
1320                                 }
1321                                 ex arg = e.op(1);
1322                                 if (parameter.op(parameter.nops()-1) == 0) {
1323
1324                                         //
1325                                         if (parameter.nops() == 1) {
1326                                                 return log(arg);
1327                                         }
1328
1329                                         //
1330                                         lst::const_iterator it = parameter.begin();
1331                                         while ((it != parameter.end()) && (*it == 0)) {
1332                                                 it++;
1333                                         }
1334                                         if (it == parameter.end()) {
1335                                                 return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
1336                                         }
1337
1338                                         //
1339                                         parameter.remove_last();
1340                                         int lastentry = parameter.nops();
1341                                         while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
1342                                                 lastentry--;
1343                                         }
1344
1345                                         //
1346                                         ex result = log(arg) * H(parameter,arg).hold();
1347                                         ex acc = 0;
1348                                         for (ex i=0; i<lastentry; i++) {
1349                                                 if (parameter[i] > 0) {
1350                                                         parameter[i]++;
1351                                                         result -= (acc + parameter[i]-1) * H(parameter, arg).hold();
1352                                                         parameter[i]--;
1353                                                         acc = 0;
1354                                                 } else if (parameter[i] < 0) {
1355                                                         parameter[i]--;
1356                                                         result -= (acc + abs(parameter[i]+1)) * H(parameter, arg).hold();
1357                                                         parameter[i]++;
1358                                                         acc = 0;
1359                                                 } else {
1360                                                         acc++;
1361                                                 }
1362                                         }
1363
1364                                         if (lastentry < parameter.nops()) {
1365                                                 result = result / (parameter.nops()-lastentry+1);
1366                                                 return result.map(*this);
1367                                         } else {
1368                                                 return result;
1369                                         }
1370                                 }
1371                         }
1372                 }
1373                 return e;
1374         }
1375 };
1376
1377
1378 // returns an expression with zeta functions corresponding to the parameter list for H
1379 ex convert_H_to_zeta(const lst& m)
1380 {
1381         symbol xtemp("xtemp");
1382         map_trafo_H_reduce_trailing_zeros filter;
1383         map_trafo_H_convert_to_zeta filter2;
1384         return filter2(filter(H(m, xtemp).hold())).subs(xtemp == 1);
1385 }
1386
1387
1388 // convert signs form Li to H representation
1389 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf)
1390 {
1391         lst res;
1392         lst::const_iterator itm = m.begin();
1393         lst::const_iterator itx = ++x.begin();
1394         ex signum = _ex1;
1395         pf = _ex1;
1396         res.append(*itm);
1397         itm++;
1398         while (itx != x.end()) {
1399                 signum *= *itx;
1400                 pf *= signum;
1401                 res.append((*itm) * signum);
1402                 itm++;
1403                 itx++;
1404         }
1405         return res;
1406 }
1407
1408
1409 // multiplies an one-dimensional H with another H
1410 // [ReV] (18)
1411 ex trafo_H_mult(const ex& h1, const ex& h2)
1412 {
1413         ex res;
1414         ex hshort;
1415         lst hlong;
1416         ex h1nops = h1.op(0).nops();
1417         ex h2nops = h2.op(0).nops();
1418         if (h1nops > 1) {
1419                 hshort = h2.op(0).op(0);
1420                 hlong = ex_to<lst>(h1.op(0));
1421         } else {
1422                 hshort = h1.op(0).op(0);
1423                 if (h2nops > 1) {
1424                         hlong = ex_to<lst>(h2.op(0));
1425                 } else {
1426                         hlong = h2.op(0).op(0);
1427                 }
1428         }
1429         for (int i=0; i<=hlong.nops(); i++) {
1430                 lst newparameter;
1431                 int j=0;
1432                 for (; j<i; j++) {
1433                         newparameter.append(hlong[j]);
1434                 }
1435                 newparameter.append(hshort);
1436                 for (; j<hlong.nops(); j++) {
1437                         newparameter.append(hlong[j]);
1438                 }
1439                 res += H(newparameter, h1.op(1)).hold();
1440         }
1441         return res;
1442 }
1443
1444
1445 // applies trafo_H_mult recursively on expressions
1446 struct map_trafo_H_mult : public map_function
1447 {
1448         ex operator()(const ex& e)
1449         {
1450                 if (is_a<add>(e)) {
1451                         return e.map(*this);
1452                 }
1453
1454                 if (is_a<mul>(e)) {
1455
1456                         ex result = 1;
1457                         ex firstH;
1458                         lst Hlst;
1459                         for (int pos=0; pos<e.nops(); pos++) {
1460                                 if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
1461                                         std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
1462                                         if (name == "H") {
1463                                                 for (ex i=0; i<e.op(pos).op(1); i++) {
1464                                                         Hlst.append(e.op(pos).op(0));
1465                                                 }
1466                                                 continue;
1467                                         }
1468                                 } else if (is_a<function>(e.op(pos))) {
1469                                         std::string name = ex_to<function>(e.op(pos)).get_name();
1470                                         if (name == "H") {
1471                                                 if (e.op(pos).op(0).nops() > 1) {
1472                                                         firstH = e.op(pos);
1473                                                 } else {
1474                                                         Hlst.append(e.op(pos));
1475                                                 }
1476                                                 continue;
1477                                         }
1478                                 }
1479                                 result *= e.op(pos);
1480                         }
1481                         if (firstH == 0) {
1482                                 if (Hlst.nops() > 0) {
1483                                         firstH = Hlst[Hlst.nops()-1];
1484                                         Hlst.remove_last();
1485                                 } else {
1486                                         return e;
1487                                 }
1488                         }
1489
1490                         if (Hlst.nops() > 0) {
1491                                 ex buffer = trafo_H_mult(firstH, Hlst.op(0));
1492                                 result *= buffer;
1493                                 for (int i=1; i<Hlst.nops(); i++) {
1494                                         result *= Hlst.op(i);
1495                                 }
1496                                 result = result.expand();
1497                                 map_trafo_H_mult recursion;
1498                                 return recursion(result);
1499                         } else {
1500                                 return e;
1501                         }
1502
1503                 }
1504                 return e;
1505         }
1506 };
1507
1508
1509 // do integration [ReV] (55)
1510 // put parameter 0 in front of existing parameters
1511 ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
1512 {
1513         ex h;
1514         std::string name;
1515         if (is_a<function>(e)) {
1516                 name = ex_to<function>(e).get_name();
1517         }
1518         if (name == "H") {
1519                 h = e;
1520         } else {
1521                 for (int i=0; i<e.nops(); i++) {
1522                         if (is_a<function>(e.op(i))) {
1523                                 std::string name = ex_to<function>(e.op(i)).get_name();
1524                                 if (name == "H") {
1525                                         h = e.op(i);
1526                                 }
1527                         }
1528                 }
1529         }
1530         if (h != 0) {
1531                 lst newparameter = ex_to<lst>(h.op(0));
1532                 newparameter.prepend(0);
1533                 ex addzeta = convert_H_to_zeta(newparameter);
1534                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
1535         } else {
1536                 return e * (-H(lst(0),1/arg).hold());
1537         }
1538 }
1539
1540
1541 // do integration [ReV] (55)
1542 // put parameter -1 in front of existing parameters
1543 ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
1544 {
1545         ex h;
1546         std::string name;
1547         if (is_a<function>(e)) {
1548                 name = ex_to<function>(e).get_name();
1549         }
1550         if (name == "H") {
1551                 h = e;
1552         } else {
1553                 for (int i=0; i<e.nops(); i++) {
1554                         if (is_a<function>(e.op(i))) {
1555                                 std::string name = ex_to<function>(e.op(i)).get_name();
1556                                 if (name == "H") {
1557                                         h = e.op(i);
1558                                 }
1559                         }
1560                 }
1561         }
1562         if (h != 0) {
1563                 lst newparameter = ex_to<lst>(h.op(0));
1564                 newparameter.prepend(-1);
1565                 ex addzeta = convert_H_to_zeta(newparameter);
1566                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
1567         } else {
1568                 ex addzeta = convert_H_to_zeta(lst(-1));
1569                 return (e * (addzeta - H(lst(-1),1/arg).hold())).expand();
1570         }
1571 }
1572
1573
1574 // do integration [ReV] (55)
1575 // put parameter -1 in front of existing parameters
1576 ex trafo_H_1mxt1px_prepend_minusone(const ex& e, const ex& arg)
1577 {
1578         ex h;
1579         std::string name;
1580         if (is_a<function>(e)) {
1581                 name = ex_to<function>(e).get_name();
1582         }
1583         if (name == "H") {
1584                 h = e;
1585         } else {
1586                 for (int i=0; i<e.nops(); i++) {
1587                         if (is_a<function>(e.op(i))) {
1588                                 std::string name = ex_to<function>(e.op(i)).get_name();
1589                                 if (name == "H") {
1590                                         h = e.op(i);
1591                                 }
1592                         }
1593                 }
1594         }
1595         if (h != 0) {
1596                 lst newparameter = ex_to<lst>(h.op(0));
1597                 newparameter.prepend(-1);
1598                 return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
1599         } else {
1600                 return (e * H(lst(-1),(1-arg)/(1+arg)).hold()).expand();
1601         }
1602 }
1603
1604
1605 // do integration [ReV] (55)
1606 // put parameter 1 in front of existing parameters
1607 ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
1608 {
1609         ex h;
1610         std::string name;
1611         if (is_a<function>(e)) {
1612                 name = ex_to<function>(e).get_name();
1613         }
1614         if (name == "H") {
1615                 h = e;
1616         } else {
1617                 for (int i=0; i<e.nops(); i++) {
1618                         if (is_a<function>(e.op(i))) {
1619                                 std::string name = ex_to<function>(e.op(i)).get_name();
1620                                 if (name == "H") {
1621                                         h = e.op(i);
1622                                 }
1623                         }
1624                 }
1625         }
1626         if (h != 0) {
1627                 lst newparameter = ex_to<lst>(h.op(0));
1628                 newparameter.prepend(1);
1629                 return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
1630         } else {
1631                 return (e * H(lst(1),(1-arg)/(1+arg)).hold()).expand();
1632         }
1633 }
1634
1635
1636 // do x -> 1/x transformation
1637 struct map_trafo_H_1overx : public map_function
1638 {
1639         ex operator()(const ex& e)
1640         {
1641                 if (is_a<add>(e) || is_a<mul>(e)) {
1642                         return e.map(*this);
1643                 }
1644
1645                 if (is_a<function>(e)) {
1646                         std::string name = ex_to<function>(e).get_name();
1647                         if (name == "H") {
1648
1649                                 lst parameter = ex_to<lst>(e.op(0));
1650                                 ex arg = e.op(1);
1651
1652                                 // special cases if all parameters are either 0, 1 or -1
1653                                 bool allthesame = true;
1654                                 if (parameter.op(0) == 0) {
1655                                         for (int i=1; i<parameter.nops(); i++) {
1656                                                 if (parameter.op(i) != 0) {
1657                                                         allthesame = false;
1658                                                         break;
1659                                                 }
1660                                         }
1661                                         if (allthesame) {
1662                                                 return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
1663                                         }
1664                                 } else if (parameter.op(0) == -1) {
1665                                         for (int i=1; i<parameter.nops(); i++) {
1666                                                 if (parameter.op(i) != -1) {
1667                                                         allthesame = false;
1668                                                         break;
1669                                                 }
1670                                         }
1671                                         if (allthesame) {
1672                                                 map_trafo_H_mult unify;
1673                                                 return unify((pow(H(lst(-1),1/arg).hold() - H(lst(0),1/arg).hold(), parameter.nops())
1674                                                        / factorial(parameter.nops())).expand());
1675                                         }
1676                                 } else {
1677                                         for (int i=1; i<parameter.nops(); i++) {
1678                                                 if (parameter.op(i) != 1) {
1679                                                         allthesame = false;
1680                                                         break;
1681                                                 }
1682                                         }
1683                                         if (allthesame) {
1684                                                 map_trafo_H_mult unify;
1685                                                 return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() + H_polesign, parameter.nops())
1686                                                        / factorial(parameter.nops())).expand());
1687                                         }
1688                                 }
1689
1690                                 lst newparameter = parameter;
1691                                 newparameter.remove_first();
1692
1693                                 if (parameter.op(0) == 0) {
1694
1695                                         // leading zero
1696                                         ex res = convert_H_to_zeta(parameter);
1697                                         map_trafo_H_1overx recursion;
1698                                         ex buffer = recursion(H(newparameter, arg).hold());
1699                                         if (is_a<add>(buffer)) {
1700                                                 for (int i=0; i<buffer.nops(); i++) {
1701                                                         res += trafo_H_1tx_prepend_zero(buffer.op(i), arg);
1702                                                 }
1703                                         } else {
1704                                                 res += trafo_H_1tx_prepend_zero(buffer, arg);
1705                                         }
1706                                         return res;
1707
1708                                 } else if (parameter.op(0) == -1) {
1709
1710                                         // leading negative one
1711                                         ex res = convert_H_to_zeta(parameter);
1712                                         map_trafo_H_1overx recursion;
1713                                         ex buffer = recursion(H(newparameter, arg).hold());
1714                                         if (is_a<add>(buffer)) {
1715                                                 for (int i=0; i<buffer.nops(); i++) {
1716                                                         res += trafo_H_1tx_prepend_zero(buffer.op(i), arg) - trafo_H_1tx_prepend_minusone(buffer.op(i), arg);
1717                                                 }
1718                                         } else {
1719                                                 res += trafo_H_1tx_prepend_zero(buffer, arg) - trafo_H_1tx_prepend_minusone(buffer, arg);
1720                                         }
1721                                         return res;
1722
1723                                 } else {
1724
1725                                         // leading one
1726                                         map_trafo_H_1overx recursion;
1727                                         map_trafo_H_mult unify;
1728                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
1729                                         int firstzero = 0;
1730                                         while (parameter.op(firstzero) == 1) {
1731                                                 firstzero++;
1732                                         }
1733                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1734                                                 lst newparameter;
1735                                                 int j=0;
1736                                                 for (; j<=i; j++) {
1737                                                         newparameter.append(parameter[j+1]);
1738                                                 }
1739                                                 newparameter.append(1);
1740                                                 for (; j<parameter.nops()-1; j++) {
1741                                                         newparameter.append(parameter[j+1]);
1742                                                 }
1743                                                 res -= H(newparameter, arg).hold();
1744                                         }
1745                                         res = recursion(res).expand() / firstzero;
1746                                         return unify(res);
1747
1748                                 }
1749
1750                         }
1751                 }
1752                 return e;
1753         }
1754 };
1755
1756
1757 // do x -> (1-x)/(1+x) transformation
1758 struct map_trafo_H_1mxt1px : public map_function
1759 {
1760         ex operator()(const ex& e)
1761         {
1762                 if (is_a<add>(e) || is_a<mul>(e)) {
1763                         return e.map(*this);
1764                 }
1765
1766                 if (is_a<function>(e)) {
1767                         std::string name = ex_to<function>(e).get_name();
1768                         if (name == "H") {
1769
1770                                 lst parameter = ex_to<lst>(e.op(0));
1771                                 ex arg = e.op(1);
1772
1773                                 // special cases if all parameters are either 0, 1 or -1
1774                                 bool allthesame = true;
1775                                 if (parameter.op(0) == 0) {
1776                                         for (int i=1; i<parameter.nops(); i++) {
1777                                                 if (parameter.op(i) != 0) {
1778                                                         allthesame = false;
1779                                                         break;
1780                                                 }
1781                                         }
1782                                         if (allthesame) {
1783                                                 map_trafo_H_mult unify;
1784                                                 return unify((pow(-H(lst(1),(1-arg)/(1+arg)).hold() - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
1785                                                        / factorial(parameter.nops())).expand());
1786                                         }
1787                                 } else if (parameter.op(0) == -1) {
1788                                         for (int i=1; i<parameter.nops(); i++) {
1789                                                 if (parameter.op(i) != -1) {
1790                                                         allthesame = false;
1791                                                         break;
1792                                                 }
1793                                         }
1794                                         if (allthesame) {
1795                                                 map_trafo_H_mult unify;
1796                                                 return unify((pow(log(2) - H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
1797                                                        / factorial(parameter.nops())).expand());
1798                                         }
1799                                 } else {
1800                                         for (int i=1; i<parameter.nops(); i++) {
1801                                                 if (parameter.op(i) != 1) {
1802                                                         allthesame = false;
1803                                                         break;
1804                                                 }
1805                                         }
1806                                         if (allthesame) {
1807                                                 map_trafo_H_mult unify;
1808                                                 return unify((pow(-log(2) - H(lst(0),(1-arg)/(1+arg)).hold() + H(lst(-1),(1-arg)/(1+arg)).hold(), parameter.nops())
1809                                                        / factorial(parameter.nops())).expand());
1810                                         }
1811                                 }
1812
1813                                 lst newparameter = parameter;
1814                                 newparameter.remove_first();
1815
1816                                 if (parameter.op(0) == 0) {
1817
1818                                         // leading zero
1819                                         ex res = convert_H_to_zeta(parameter);
1820                                         map_trafo_H_1mxt1px recursion;
1821                                         ex buffer = recursion(H(newparameter, arg).hold());
1822                                         if (is_a<add>(buffer)) {
1823                                                 for (int i=0; i<buffer.nops(); i++) {
1824                                                         res -= trafo_H_1mxt1px_prepend_one(buffer.op(i), arg) + trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
1825                                                 }
1826                                         } else {
1827                                                 res -= trafo_H_1mxt1px_prepend_one(buffer, arg) + trafo_H_1mxt1px_prepend_minusone(buffer, arg);
1828                                         }
1829                                         return res;
1830
1831                                 } else if (parameter.op(0) == -1) {
1832
1833                                         // leading negative one
1834                                         ex res = convert_H_to_zeta(parameter);
1835                                         map_trafo_H_1mxt1px recursion;
1836                                         ex buffer = recursion(H(newparameter, arg).hold());
1837                                         if (is_a<add>(buffer)) {
1838                                                 for (int i=0; i<buffer.nops(); i++) {
1839                                                         res -= trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
1840                                                 }
1841                                         } else {
1842                                                 res -= trafo_H_1mxt1px_prepend_minusone(buffer, arg);
1843                                         }
1844                                         return res;
1845
1846                                 } else {
1847
1848                                         // leading one
1849                                         map_trafo_H_1mxt1px recursion;
1850                                         map_trafo_H_mult unify;
1851                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
1852                                         int firstzero = 0;
1853                                         while (parameter.op(firstzero) == 1) {
1854                                                 firstzero++;
1855                                         }
1856                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1857                                                 lst newparameter;
1858                                                 int j=0;
1859                                                 for (; j<=i; j++) {
1860                                                         newparameter.append(parameter[j+1]);
1861                                                 }
1862                                                 newparameter.append(1);
1863                                                 for (; j<parameter.nops()-1; j++) {
1864                                                         newparameter.append(parameter[j+1]);
1865                                                 }
1866                                                 res -= H(newparameter, arg).hold();
1867                                         }
1868                                         res = recursion(res).expand() / firstzero;
1869                                         return unify(res);
1870
1871                                 }
1872
1873                         }
1874                 }
1875                 return e;
1876         }
1877 };
1878
1879
1880 // do the actual summation.
1881 cln::cl_N H_do_sum(const std::vector<int>& m, const cln::cl_N& x)
1882 {
1883         const int j = m.size();
1884
1885         std::vector<cln::cl_N> t(j);
1886
1887         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1888         cln::cl_N factor = cln::expt(x, j) * one;
1889         cln::cl_N t0buf;
1890         int q = 0;
1891         do {
1892                 t0buf = t[0];
1893                 q++;
1894                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),m[j-1]);
1895                 for (int k=j-2; k>=1; k--) {
1896                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), m[k]);
1897                 }
1898                 t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), m[0]);
1899                 factor = factor * x;
1900         } while (t[0] != t0buf);
1901
1902         return t[0];
1903 }
1904
1905
1906 } // end of anonymous namespace
1907
1908
1909 //////////////////////////////////////////////////////////////////////
1910 //
1911 // Harmonic polylogarithm  H(m,x)
1912 //
1913 // GiNaC function
1914 //
1915 //////////////////////////////////////////////////////////////////////
1916
1917
1918 static ex H_evalf(const ex& x1, const ex& x2)
1919 {
1920         if (is_a<lst>(x1)) {
1921
1922                 cln::cl_N x;
1923                 if (is_a<numeric>(x2)) {
1924                         x = ex_to<numeric>(x2).to_cl_N();
1925                 } else {
1926                         ex x2_val = x2.evalf();
1927                         if (is_a<numeric>(x2_val)) {
1928                                 x = ex_to<numeric>(x2_val).to_cl_N();
1929                         }
1930                 }
1931
1932                 for (int i=0; i<x1.nops(); i++) {
1933                         if (!x1.op(i).info(info_flags::integer)) {
1934                                 return H(x1, x2).hold();
1935                         }
1936                 }
1937                 if (x1.nops() < 1) {
1938                         return H(x1, x2).hold();
1939                 }
1940
1941                 const lst& morg = ex_to<lst>(x1);
1942                 // remove trailing zeros ...
1943                 if (*(--morg.end()) == 0) {
1944                         symbol xtemp("xtemp");
1945                         map_trafo_H_reduce_trailing_zeros filter;
1946                         return filter(H(x1, xtemp).hold()).subs(xtemp==x2).evalf();
1947                 }
1948                 // ... and expand parameter notation
1949                 lst m;
1950                 for (lst::const_iterator it = morg.begin(); it != morg.end(); it++) {
1951                         if (*it > 1) {
1952                                 for (ex count=*it-1; count > 0; count--) {
1953                                         m.append(0);
1954                                 }
1955                                 m.append(1);
1956                         } else if (*it < -1) {
1957                                 for (ex count=*it+1; count < 0; count++) {
1958                                         m.append(0);
1959                                 }
1960                                 m.append(-1);
1961                         } else {
1962                                 m.append(*it);
1963                         }
1964                 }
1965
1966                 // since the transformations produce a lot of terms, they are only efficient for
1967                 // argument near one.
1968                 // no transformation needed -> do summation
1969                 if (cln::abs(x) < 0.95) {
1970                         lst m_lst;
1971                         lst s_lst;
1972                         ex pf;
1973                         if (convert_parameter_H_to_Li(m, m_lst, s_lst, pf)) {
1974                                 // negative parameters -> s_lst is filled
1975                                 std::vector<int> m_int;
1976                                 std::vector<cln::cl_N> x_cln;
1977                                 for (lst::const_iterator it_int = m_lst.begin(), it_cln = s_lst.begin();
1978                                      it_int != m_lst.end(); it_int++, it_cln++) {
1979                                         m_int.push_back(ex_to<numeric>(*it_int).to_int());
1980                                         x_cln.push_back(ex_to<numeric>(*it_cln).to_cl_N());
1981                                 }
1982                                 x_cln.front() = x_cln.front() * x;
1983                                 return pf * numeric(multipleLi_do_sum(m_int, x_cln));
1984                         } else {
1985                                 // only positive parameters
1986                                 //TODO
1987                                 if (m_lst.nops() == 1) {
1988                                         return Li(m_lst.op(0), x2).evalf();
1989                                 }
1990                                 std::vector<int> m_int;
1991                                 for (lst::const_iterator it = m_lst.begin(); it != m_lst.end(); it++) {
1992                                         m_int.push_back(ex_to<numeric>(*it).to_int());
1993                                 }
1994                                 return numeric(H_do_sum(m_int, x));
1995                         }
1996                 }
1997
1998                 ex res = 1;
1999
2000                 // ensure that the realpart of the argument is positive
2001                 if (cln::realpart(x) < 0) {
2002                         x = -x;
2003                         for (int i=0; i<m.nops(); i++) {
2004                                 if (m.op(i) != 0) {
2005                                         m.let_op(i) = -m.op(i);
2006                                         res *= -1;
2007                                 }
2008                         }
2009                 }
2010
2011                 // choose transformations
2012                 symbol xtemp("xtemp");
2013                 if (cln::abs(x-1) < 1.4142) {
2014                         // x -> (1-x)/(1+x)
2015                         map_trafo_H_1mxt1px trafo;
2016                         res *= trafo(H(m, xtemp));
2017                 } else {
2018                         // x -> 1/x
2019                         map_trafo_H_1overx trafo;
2020                         res *= trafo(H(m, xtemp));
2021                         if (cln::imagpart(x) <= 0) {
2022                                 res = res.subs(H_polesign == -I*Pi);
2023                         } else {
2024                                 res = res.subs(H_polesign == I*Pi);
2025                         }
2026                 }
2027
2028                 // simplify result
2029 // TODO
2030 //              map_trafo_H_convert converter;
2031 //              res = converter(res).expand();
2032 //              lst ll;
2033 //              res.find(H(wild(1),wild(2)), ll);
2034 //              res.find(zeta(wild(1)), ll);
2035 //              res.find(zeta(wild(1),wild(2)), ll);
2036 //              res = res.collect(ll);
2037
2038                 return res.subs(xtemp == numeric(x)).evalf();
2039         }
2040
2041         return H(x1,x2).hold();
2042 }
2043
2044
2045 static ex H_eval(const ex& m_, const ex& x)
2046 {
2047         lst m;
2048         if (is_a<lst>(m_)) {
2049                 m = ex_to<lst>(m_);
2050         } else {
2051                 m = lst(m_);
2052         }
2053         if (m.nops() == 0) {
2054                 return _ex1;
2055         }
2056         ex pos1;
2057         ex pos2;
2058         ex n;
2059         ex p;
2060         int step = 0;
2061         if (*m.begin() > _ex1) {
2062                 step++;
2063                 pos1 = _ex0;
2064                 pos2 = _ex1;
2065                 n = *m.begin()-1;
2066                 p = _ex1;
2067         } else if (*m.begin() < _ex_1) {
2068                 step++;
2069                 pos1 = _ex0;
2070                 pos2 = _ex_1;
2071                 n = -*m.begin()-1;
2072                 p = _ex1;
2073         } else if (*m.begin() == _ex0) {
2074                 pos1 = _ex0;
2075                 n = _ex1;
2076         } else {
2077                 pos1 = *m.begin();
2078                 p = _ex1;
2079         }
2080         for (lst::const_iterator it = ++m.begin(); it != m.end(); it++) {
2081                 if ((*it).info(info_flags::integer)) {
2082                         if (step == 0) {
2083                                 if (*it > _ex1) {
2084                                         if (pos1 == _ex0) {
2085                                                 step = 1;
2086                                                 pos2 = _ex1;
2087                                                 n += *it-1;
2088                                                 p = _ex1;
2089                                         } else {
2090                                                 step = 2;
2091                                         }
2092                                 } else if (*it < _ex_1) {
2093                                         if (pos1 == _ex0) {
2094                                                 step = 1;
2095                                                 pos2 = _ex_1;
2096                                                 n += -*it-1;
2097                                                 p = _ex1;
2098                                         } else {
2099                                                 step = 2;
2100                                         }
2101                                 } else {
2102                                         if (*it != pos1) {
2103                                                 step = 1;
2104                                                 pos2 = *it;
2105                                         }
2106                                         if (*it == _ex0) {
2107                                                 n++;
2108                                         } else {
2109                                                 p++;
2110                                         }
2111                                 }
2112                         } else if (step == 1) {
2113                                 if (*it != pos2) {
2114                                         step = 2;
2115                                 } else {
2116                                         if (*it == _ex0) {
2117                                                 n++;
2118                                         } else {
2119                                                 p++;
2120                                         }
2121                                 }
2122                         }
2123                 } else {
2124                         // if some m_i is not an integer
2125                         return H(m_, x).hold();
2126                 }
2127         }
2128         if ((x == _ex1) && (*(--m.end()) != _ex0)) {
2129                 return convert_H_to_zeta(m);
2130         }
2131         if (step == 0) {
2132                 if (pos1 == _ex0) {
2133                         // all zero
2134                         if (x == _ex0) {
2135                                 return H(m_, x).hold();
2136                         }
2137                         return pow(log(x), m.nops()) / factorial(m.nops());
2138                 } else {
2139                         // all (minus) one
2140                         return pow(-pos1*log(1-pos1*x), m.nops()) / factorial(m.nops());
2141                 }
2142         } else if ((step == 1) && (pos1 == _ex0)){
2143                 // convertible to S
2144                 if (pos2 == _ex1) {
2145                         return S(n, p, x);
2146                 } else {
2147                         return pow(-1, p) * S(n, p, -x);
2148                 }
2149         }
2150         if (x == _ex0) {
2151                 return _ex0;
2152         }
2153         if (x.info(info_flags::numeric) && (!x.info(info_flags::crational))) {
2154                 return H(m_, x).evalf();
2155         }
2156         return H(m_, x).hold();
2157 }
2158
2159
2160 static ex H_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
2161 {
2162         epvector seq;
2163         seq.push_back(expair(H(m, x), 0));
2164         return pseries(rel, seq);
2165 }
2166
2167
2168 static ex H_deriv(const ex& m_, const ex& x, unsigned deriv_param)
2169 {
2170         GINAC_ASSERT(deriv_param < 2);
2171         if (deriv_param == 0) {
2172                 return _ex0;
2173         }
2174         lst m;
2175         if (is_a<lst>(m_)) {
2176                 m = ex_to<lst>(m_);
2177         } else {
2178                 m = lst(m_);
2179         }
2180         ex mb = *m.begin();
2181         if (mb > _ex1) {
2182                 m[0]--;
2183                 return H(m, x) / x;
2184         }
2185         if (mb < _ex_1) {
2186                 m[0]++;
2187                 return H(m, x) / x;
2188         }
2189         m.remove_first();
2190         if (mb == _ex1) {
2191                 return 1/(1-x) * H(m, x);
2192         } else if (mb == _ex_1) {
2193                 return 1/(1+x) * H(m, x);
2194         } else {
2195                 return H(m, x) / x;
2196         }
2197 }
2198
2199
2200 static void H_print_latex(const ex& m_, const ex& x, const print_context& c)
2201 {
2202         lst m;
2203         if (is_a<lst>(m_)) {
2204                 m = ex_to<lst>(m_);
2205         } else {
2206                 m = lst(m_);
2207         }
2208         c.s << "\\mbox{H}_{";
2209         lst::const_iterator itm = m.begin();
2210         (*itm).print(c);
2211         itm++;
2212         for (; itm != m.end(); itm++) {
2213                 c.s << ",";
2214                 (*itm).print(c);
2215         }
2216         c.s << "}(";
2217         x.print(c);
2218         c.s << ")";
2219 }
2220
2221
2222 REGISTER_FUNCTION(H,
2223                   evalf_func(H_evalf).
2224                   eval_func(H_eval).
2225                   series_func(H_series).
2226                   derivative_func(H_deriv).
2227                   print_func<print_latex>(H_print_latex).
2228                   do_not_evalf_params());
2229
2230
2231 // takes a parameter list for H and returns an expression with corresponding multiple polylogarithms
2232 ex convert_H_to_Li(const ex& m, const ex& x)
2233 {
2234         map_trafo_H_reduce_trailing_zeros filter;
2235         map_trafo_H_convert_to_Li filter2;
2236         if (is_a<lst>(m)) {
2237                 return filter2(filter(H(m, x).hold()));
2238         } else {
2239                 return filter2(filter(H(lst(m), x).hold()));
2240         }
2241 }
2242
2243
2244 //////////////////////////////////////////////////////////////////////
2245 //
2246 // Multiple zeta values  zeta(x) and zeta(x,s)
2247 //
2248 // helper functions
2249 //
2250 //////////////////////////////////////////////////////////////////////
2251
2252
2253 // anonymous namespace for helper functions
2254 namespace {
2255
2256
2257 // parameters and data for [Cra] algorithm
2258 const cln::cl_N lambda = cln::cl_N("319/320");
2259 int L1;
2260 int L2;
2261 std::vector<std::vector<cln::cl_N> > f_kj;
2262 std::vector<cln::cl_N> crB;
2263 std::vector<std::vector<cln::cl_N> > crG;
2264 std::vector<cln::cl_N> crX;
2265
2266
2267 void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
2268 {
2269         const int size = a.size();
2270         for (int n=0; n<size; n++) {
2271                 c[n] = 0;
2272                 for (int m=0; m<=n; m++) {
2273                         c[n] = c[n] + a[m]*b[n-m];
2274                 }
2275         }
2276 }
2277
2278
2279 // [Cra] section 4
2280 void initcX(const std::vector<int>& s)
2281 {
2282         const int k = s.size();
2283
2284         crX.clear();
2285         crG.clear();
2286         crB.clear();
2287
2288         for (int i=0; i<=L2; i++) {
2289                 crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
2290         }
2291
2292         int Sm = 0;
2293         int Smp1 = 0;
2294         for (int m=0; m<k-1; m++) {
2295                 std::vector<cln::cl_N> crGbuf;
2296                 Sm = Sm + s[m];
2297                 Smp1 = Sm + s[m+1];
2298                 for (int i=0; i<=L2; i++) {
2299                         crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
2300                 }
2301                 crG.push_back(crGbuf);
2302         }
2303
2304         crX = crB;
2305
2306         for (int m=0; m<k-1; m++) {
2307                 std::vector<cln::cl_N> Xbuf;
2308                 for (int i=0; i<=L2; i++) {
2309                         Xbuf.push_back(crX[i] * crG[m][i]);
2310                 }
2311                 halfcyclic_convolute(Xbuf, crB, crX);
2312         }
2313 }
2314
2315
2316 // [Cra] section 4
2317 cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
2318 {
2319         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
2320         cln::cl_N factor = cln::expt(lambda, Sqk);
2321         cln::cl_N res = factor / Sqk * crX[0] * one;
2322         cln::cl_N resbuf;
2323         int N = 0;
2324         do {
2325                 resbuf = res;
2326                 factor = factor * lambda;
2327                 N++;
2328                 res = res + crX[N] * factor / (N+Sqk);
2329         } while ((res != resbuf) || cln::zerop(crX[N]));
2330         return res;
2331 }
2332
2333
2334 // [Cra] section 4
2335 void calc_f(int maxr)
2336 {
2337         f_kj.clear();
2338         f_kj.resize(L1);
2339
2340         cln::cl_N t0, t1, t2, t3, t4;
2341         int i, j, k;
2342         std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
2343         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
2344
2345         t0 = cln::exp(-lambda);
2346         t2 = 1;
2347         for (k=1; k<=L1; k++) {
2348                 t1 = k * lambda;
2349                 t2 = t0 * t2;
2350                 for (j=1; j<=maxr; j++) {
2351                         t3 = 1;
2352                         t4 = 1;
2353                         for (i=2; i<=j; i++) {
2354                                 t4 = t4 * (j-i+1);
2355                                 t3 = t1 * t3 + t4;
2356                         }
2357                         (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one);
2358                 }
2359                 it++;
2360         }
2361 }
2362
2363
2364 // [Cra] (3.1)
2365 cln::cl_N crandall_Z(const std::vector<int>& s)
2366 {
2367         const int j = s.size();
2368
2369         if (j == 1) {
2370                 cln::cl_N t0;
2371                 cln::cl_N t0buf;
2372                 int q = 0;
2373                 do {
2374                         t0buf = t0;
2375                         q++;
2376                         t0 = t0 + f_kj[q+j-2][s[0]-1];
2377                 } while (t0 != t0buf);
2378
2379                 return t0 / cln::factorial(s[0]-1);
2380         }
2381
2382         std::vector<cln::cl_N> t(j);
2383
2384         cln::cl_N t0buf;
2385         int q = 0;
2386         do {
2387                 t0buf = t[0];
2388                 q++;
2389                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
2390                 for (int k=j-2; k>=1; k--) {
2391                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
2392                 }
2393                 t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
2394         } while (t[0] != t0buf);
2395
2396         return t[0] / cln::factorial(s[0]-1);
2397 }
2398
2399
2400 // [Cra] (2.4)
2401 cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
2402 {
2403         std::vector<int> r = s;
2404         const int j = r.size();
2405
2406         // decide on maximal size of f_kj for crandall_Z
2407         if (Digits < 50) {
2408                 L1 = 150;
2409         } else {
2410                 L1 = Digits * 3 + j*2;
2411         }
2412
2413         // decide on maximal size of crX for crandall_Y
2414         if (Digits < 38) {
2415                 L2 = 63;
2416         } else if (Digits < 86) {
2417                 L2 = 127;
2418         } else if (Digits < 192) {
2419                 L2 = 255;
2420         } else if (Digits < 394) {
2421                 L2 = 511;
2422         } else if (Digits < 808) {
2423                 L2 = 1023;
2424         } else {
2425                 L2 = 2047;
2426         }
2427
2428         cln::cl_N res;
2429
2430         int maxr = 0;
2431         int S = 0;
2432         for (int i=0; i<j; i++) {
2433                 S += r[i];
2434                 if (r[i] > maxr) {
2435                         maxr = r[i];
2436                 }
2437         }
2438
2439         calc_f(maxr);
2440
2441         const cln::cl_N r0factorial = cln::factorial(r[0]-1);
2442
2443         std::vector<int> rz;
2444         int skp1buf;
2445         int Srun = S;
2446         for (int k=r.size()-1; k>0; k--) {
2447
2448                 rz.insert(rz.begin(), r.back());
2449                 skp1buf = rz.front();
2450                 Srun -= skp1buf;
2451                 r.pop_back();
2452
2453                 initcX(r);
2454
2455                 for (int q=0; q<skp1buf; q++) {
2456
2457                         cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
2458                         cln::cl_N pp2 = crandall_Z(rz);
2459
2460                         rz.front()--;
2461
2462                         if (q & 1) {
2463                                 res = res - pp1 * pp2 / cln::factorial(q);
2464                         } else {
2465                                 res = res + pp1 * pp2 / cln::factorial(q);
2466                         }
2467                 }
2468                 rz.front() = skp1buf;
2469         }
2470         rz.insert(rz.begin(), r.back());
2471
2472         initcX(rz);
2473
2474         res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
2475
2476         return res;
2477 }
2478
2479
2480 cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
2481 {
2482         const int j = r.size();
2483
2484         // buffer for subsums
2485         std::vector<cln::cl_N> t(j);
2486         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
2487
2488         cln::cl_N t0buf;
2489         int q = 0;
2490         do {
2491                 t0buf = t[0];
2492                 q++;
2493                 t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
2494                 for (int k=j-2; k>=0; k--) {
2495                         t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
2496                 }
2497         } while (t[0] != t0buf);
2498
2499         return t[0];
2500 }
2501
2502
2503 // does Hoelder convolution. see [BBB] (7.0)
2504 cln::cl_N zeta_do_Hoelder_convolution(const std::vector<int>& m_, const std::vector<int>& s_)
2505 {
2506         // prepare parameters
2507         // holds Li arguments in [BBB] notation
2508         std::vector<int> s = s_;
2509         std::vector<int> m_p = m_;
2510         std::vector<int> m_q;
2511         // holds Li arguments in nested sums notation
2512         std::vector<cln::cl_N> s_p(s.size(), cln::cl_N(1));
2513         s_p[0] = s_p[0] * cln::cl_N("1/2");
2514         // convert notations
2515         int sig = 1;
2516         for (int i=0; i<s_.size(); i++) {
2517                 if (s_[i] < 0) {
2518                         sig = -sig;
2519                         s_p[i] = -s_p[i];
2520                 }
2521                 s[i] = sig * std::abs(s[i]);
2522         }
2523         std::vector<cln::cl_N> s_q;
2524         cln::cl_N signum = 1;
2525
2526         // first term
2527         cln::cl_N res = multipleLi_do_sum(m_p, s_p);
2528
2529         // middle terms
2530         do {
2531
2532                 // change parameters
2533                 if (s.front() > 0) {
2534                         if (m_p.front() == 1) {
2535                                 m_p.erase(m_p.begin());
2536                                 s_p.erase(s_p.begin());
2537                                 if (s_p.size() > 0) {
2538                                         s_p.front() = s_p.front() * cln::cl_N("1/2");
2539                                 }
2540                                 s.erase(s.begin());
2541                                 m_q.front()++;
2542                         } else {
2543                                 m_p.front()--;
2544                                 m_q.insert(m_q.begin(), 1);
2545                                 if (s_q.size() > 0) {
2546                                         s_q.front() = s_q.front() * 2;
2547                                 }
2548                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
2549                         }
2550                 } else {
2551                         if (m_p.front() == 1) {
2552                                 m_p.erase(m_p.begin());
2553                                 cln::cl_N spbuf = s_p.front();
2554                                 s_p.erase(s_p.begin());
2555                                 if (s_p.size() > 0) {
2556                                         s_p.front() = s_p.front() * spbuf;
2557                                 }
2558                                 s.erase(s.begin());
2559                                 m_q.insert(m_q.begin(), 1);
2560                                 if (s_q.size() > 0) {
2561                                         s_q.front() = s_q.front() * 4;
2562                                 }
2563                                 s_q.insert(s_q.begin(), cln::cl_N("1/4"));
2564                                 signum = -signum;
2565                         } else {
2566                                 m_p.front()--;
2567                                 m_q.insert(m_q.begin(), 1);
2568                                 if (s_q.size() > 0) {
2569                                         s_q.front() = s_q.front() * 2;
2570                                 }
2571                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
2572                         }
2573                 }
2574
2575                 // exiting the loop
2576                 if (m_p.size() == 0) break;
2577
2578                 res = res + signum * multipleLi_do_sum(m_p, s_p) * multipleLi_do_sum(m_q, s_q);
2579
2580         } while (true);
2581
2582         // last term
2583         res = res + signum * multipleLi_do_sum(m_q, s_q);
2584
2585         return res;
2586 }
2587
2588
2589 } // end of anonymous namespace
2590
2591
2592 //////////////////////////////////////////////////////////////////////
2593 //
2594 // Multiple zeta values  zeta(x)
2595 //
2596 // GiNaC function
2597 //
2598 //////////////////////////////////////////////////////////////////////
2599
2600
2601 static ex zeta1_evalf(const ex& x)
2602 {
2603         if (is_exactly_a<lst>(x) && (x.nops()>1)) {
2604
2605                 // multiple zeta value
2606                 const int count = x.nops();
2607                 const lst& xlst = ex_to<lst>(x);
2608                 std::vector<int> r(count);
2609
2610                 // check parameters and convert them
2611                 lst::const_iterator it1 = xlst.begin();
2612                 std::vector<int>::iterator it2 = r.begin();
2613                 do {
2614                         if (!(*it1).info(info_flags::posint)) {
2615                                 return zeta(x).hold();
2616                         }
2617                         *it2 = ex_to<numeric>(*it1).to_int();
2618                         it1++;
2619                         it2++;
2620                 } while (it2 != r.end());
2621
2622                 // check for divergence
2623                 if (r[0] == 1) {
2624                         return zeta(x).hold();
2625                 }
2626
2627                 // decide on summation algorithm
2628                 // this is still a bit clumsy
2629                 int limit = (Digits>17) ? 10 : 6;
2630                 if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) {
2631                         return numeric(zeta_do_sum_Crandall(r));
2632                 } else {
2633                         return numeric(zeta_do_sum_simple(r));
2634                 }
2635         }
2636
2637         // single zeta value
2638         if (is_exactly_a<numeric>(x) && (x != 1)) {
2639                 try {
2640                         return zeta(ex_to<numeric>(x));
2641                 } catch (const dunno &e) { }
2642         }
2643
2644         return zeta(x).hold();
2645 }
2646
2647
2648 static ex zeta1_eval(const ex& m)
2649 {
2650         if (is_exactly_a<lst>(m)) {
2651                 if (m.nops() == 1) {
2652                         return zeta(m.op(0));
2653                 }
2654                 return zeta(m).hold();
2655         }
2656
2657         if (m.info(info_flags::numeric)) {
2658                 const numeric& y = ex_to<numeric>(m);
2659                 // trap integer arguments:
2660                 if (y.is_integer()) {
2661                         if (y.is_zero()) {
2662                                 return _ex_1_2;
2663                         }
2664                         if (y.is_equal(_num1)) {
2665                                 return zeta(m).hold();
2666                         }
2667                         if (y.info(info_flags::posint)) {
2668                                 if (y.info(info_flags::odd)) {
2669                                         return zeta(m).hold();
2670                                 } else {
2671                                         return abs(bernoulli(y)) * pow(Pi, y) * pow(_num2, y-_num1) / factorial(y);
2672                                 }
2673                         } else {
2674                                 if (y.info(info_flags::odd)) {
2675                                         return -bernoulli(_num1-y) / (_num1-y);
2676                                 } else {
2677                                         return _ex0;
2678                                 }
2679                         }
2680                 }
2681                 // zeta(float)
2682                 if (y.info(info_flags::numeric) && !y.info(info_flags::crational)) {
2683                         return zeta1_evalf(m);
2684                 }
2685         }
2686         return zeta(m).hold();
2687 }
2688
2689
2690 static ex zeta1_deriv(const ex& m, unsigned deriv_param)
2691 {
2692         GINAC_ASSERT(deriv_param==0);
2693
2694         if (is_exactly_a<lst>(m)) {
2695                 return _ex0;
2696         } else {
2697                 return zetaderiv(_ex1, m);
2698         }
2699 }
2700
2701
2702 static void zeta1_print_latex(const ex& m_, const print_context& c)
2703 {
2704         c.s << "\\zeta(";
2705         if (is_a<lst>(m_)) {
2706                 const lst& m = ex_to<lst>(m_);
2707                 lst::const_iterator it = m.begin();
2708                 (*it).print(c);
2709                 it++;
2710                 for (; it != m.end(); it++) {
2711                         c.s << ",";
2712                         (*it).print(c);
2713                 }
2714         } else {
2715                 m_.print(c);
2716         }
2717         c.s << ")";
2718 }
2719
2720
2721 unsigned zeta1_SERIAL::serial = function::register_new(function_options("zeta", 1).
2722                                 evalf_func(zeta1_evalf).
2723                                 eval_func(zeta1_eval).
2724                                 derivative_func(zeta1_deriv).
2725                                 print_func<print_latex>(zeta1_print_latex).
2726                                 do_not_evalf_params().
2727                                 overloaded(2));
2728
2729
2730 //////////////////////////////////////////////////////////////////////
2731 //
2732 // Alternating Euler sum  zeta(x,s)
2733 //
2734 // GiNaC function
2735 //
2736 //////////////////////////////////////////////////////////////////////
2737
2738
2739 static ex zeta2_evalf(const ex& x, const ex& s)
2740 {
2741         if (is_exactly_a<lst>(x)) {
2742
2743                 // alternating Euler sum
2744                 const int count = x.nops();
2745                 const lst& xlst = ex_to<lst>(x);
2746                 const lst& slst = ex_to<lst>(s);
2747                 std::vector<int> xi(count);
2748                 std::vector<int> si(count);
2749
2750                 // check parameters and convert them
2751                 lst::const_iterator it_xread = xlst.begin();
2752                 lst::const_iterator it_sread = slst.begin();
2753                 std::vector<int>::iterator it_xwrite = xi.begin();
2754                 std::vector<int>::iterator it_swrite = si.begin();
2755                 do {
2756                         if (!(*it_xread).info(info_flags::posint)) {
2757                                 return zeta(x, s).hold();
2758                         }
2759                         *it_xwrite = ex_to<numeric>(*it_xread).to_int();
2760                         if (*it_sread > 0) {
2761                                 *it_swrite = 1;
2762                         } else {
2763                                 *it_swrite = -1;
2764                         }
2765                         it_xread++;
2766                         it_sread++;
2767                         it_xwrite++;
2768                         it_swrite++;
2769                 } while (it_xwrite != xi.end());
2770
2771                 // check for divergence
2772                 if ((xi[0] == 1) && (si[0] == 1)) {
2773                         return zeta(x, s).hold();
2774                 }
2775
2776                 // use Hoelder convolution
2777                 return numeric(zeta_do_Hoelder_convolution(xi, si));
2778         }
2779
2780         return zeta(x, s).hold();
2781 }
2782
2783
2784 static ex zeta2_eval(const ex& m, const ex& s_)
2785 {
2786         if (is_exactly_a<lst>(s_)) {
2787                 const lst& s = ex_to<lst>(s_);
2788                 for (lst::const_iterator it = s.begin(); it != s.end(); it++) {
2789                         if ((*it).info(info_flags::positive)) {
2790                                 continue;
2791                         }
2792                         return zeta(m, s_).hold();
2793                 }
2794                 return zeta(m);
2795         } else if (s_.info(info_flags::positive)) {
2796                 return zeta(m);
2797         }
2798
2799         return zeta(m, s_).hold();
2800 }
2801
2802
2803 static ex zeta2_deriv(const ex& m, const ex& s, unsigned deriv_param)
2804 {
2805         GINAC_ASSERT(deriv_param==0);
2806
2807         if (is_exactly_a<lst>(m)) {
2808                 return _ex0;
2809         } else {
2810                 if ((is_exactly_a<lst>(s) && s.op(0).info(info_flags::positive)) || s.info(info_flags::positive)) {
2811                         return zetaderiv(_ex1, m);
2812                 }
2813                 return _ex0;
2814         }
2815 }
2816
2817
2818 static void zeta2_print_latex(const ex& m_, const ex& s_, const print_context& c)
2819 {
2820         lst m;
2821         if (is_a<lst>(m_)) {
2822                 m = ex_to<lst>(m_);
2823         } else {
2824                 m = lst(m_);
2825         }
2826         lst s;
2827         if (is_a<lst>(s_)) {
2828                 s = ex_to<lst>(s_);
2829         } else {
2830                 s = lst(s_);
2831         }
2832         c.s << "\\zeta(";
2833         lst::const_iterator itm = m.begin();
2834         lst::const_iterator its = s.begin();
2835         if (*its < 0) {
2836                 c.s << "\\overline{";
2837                 (*itm).print(c);
2838                 c.s << "}";
2839         } else {
2840                 (*itm).print(c);
2841         }
2842         its++;
2843         itm++;
2844         for (; itm != m.end(); itm++, its++) {
2845                 c.s << ",";
2846                 if (*its < 0) {
2847                         c.s << "\\overline{";
2848                         (*itm).print(c);
2849                         c.s << "}";
2850                 } else {
2851                         (*itm).print(c);
2852                 }
2853         }
2854         c.s << ")";
2855 }
2856
2857
2858 unsigned zeta2_SERIAL::serial = function::register_new(function_options("zeta", 2).
2859                                 evalf_func(zeta2_evalf).
2860                                 eval_func(zeta2_eval).
2861                                 derivative_func(zeta2_deriv).
2862                                 print_func<print_latex>(zeta2_print_latex).
2863                                 do_not_evalf_params().
2864                                 overloaded(2));
2865
2866
2867 } // namespace GiNaC
2868