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