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                 q++;
 411                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
 412                 for (int k=j-2; k>=0; k--) {
 413                         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]);
 414                 }
 415         } while ((t[0] != t0buf) || (q<10));
 416
 417         return t[0];
 418 }
 419
 420
 421 } // end of anonymous namespace
 422
 423
 424 //////////////////////////////////////////////////////////////////////
 425 //
 426 // Classical polylogarithm and multiple polylogarithm  Li(n,x)
 427 //
 428 // GiNaC function
 429 //
 430 //////////////////////////////////////////////////////////////////////
 431
 432
 433 static ex Li_eval(const ex& x1, const ex& x2)
 434 {
 435         if (x2.is_zero()) {
 436                 return _ex0;
 437         }
 438         else {
 439                 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational)))
 440                         return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
 441                 if (is_a<lst>(x2)) {
 442                         for (int i=0; i<x2.nops(); i++) {
 443                                 if (!is_a<numeric>(x2.op(i))) {
 444                                         return Li(x1,x2).hold();
 445                                 }
 446                         }
 447                         return Li(x1,x2).evalf();
 448                 }
 449                 return Li(x1,x2).hold();
 450         }
 451 }
 452
 453
 454 static ex Li_evalf(const ex& x1, const ex& x2)
 455 {
 456         // classical polylogs
 457         if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
 458                 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
 459         }
 460         // multiple polylogs
 461         else if (is_a<lst>(x1) && is_a<lst>(x2)) {
 462                 ex conv = 1;
 463                 for (int i=0; i<x1.nops(); i++) {
 464                         if (!x1.op(i).info(info_flags::posint)) {
 465                                 return Li(x1,x2).hold();
 466                         }
 467                         if (!is_a<numeric>(x2.op(i))) {
 468                                 return Li(x1,x2).hold();
 469                         }
 470                         conv *= x2.op(i);
 471                         if ((conv > 1) || ((conv == 1) && (x1.op(0) == 1))) {
 472                                 return Li(x1,x2).hold();
 473                         }
 474                 }
 475
 476                 std::vector<int> m;
 477                 std::vector<cln::cl_N> x;
 478                 for (int i=0; i<ex_to<numeric>(x1.nops()).to_int(); i++) {
 479                         m.push_back(ex_to<numeric>(x1.op(i)).to_int());
 480                         x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
 481                 }
 482
 483                 return numeric(multipleLi_do_sum(m, x));
 484         }
 485
 486         return Li(x1,x2).hold();
 487 }
 488
 489
 490 static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
 491 {
 492         epvector seq;
 493         seq.push_back(expair(Li(x1,x2), 0));
 494         return pseries(rel,seq);
 495 }
 496
 497
 498 static ex Li_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
 499 {
 500         GINAC_ASSERT(deriv_param < 2);
 501         if (deriv_param == 0) {
 502                 return _ex0;
 503         }
 504         if (x1 > 0) {
 505                 return Li(x1-1, x2) / x2;
 506         } else {
 507                 return 1/(1-x2);
 508         }
 509 }
 510
 511
 512 REGISTER_FUNCTION(Li,
 513                 eval_func(Li_eval).
 514                 evalf_func(Li_evalf).
 515                 do_not_evalf_params().
 516                 series_func(Li_series).
 517                 derivative_func(Li_deriv));
 518
 519
 520 //////////////////////////////////////////////////////////////////////
 521 //
 522 // Nielsen's generalized polylogarithm  S(n,p,x)
 523 //
 524 // helper functions
 525 //
 526 //////////////////////////////////////////////////////////////////////
 527
 528
 529 // anonymous namespace for helper functions
 530 namespace {
 531
 532
 533 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
 534 // see fill_Yn()
 535 std::vector<std::vector<cln::cl_N> > Yn;
 536 int ynsize = 0; // number of Yn[]
 537 int ynlength = 100; // initial length of all Yn[i]
 538
 539
 540 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
 541 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
 542 // representing S_{n,p}(x).
 543 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
 544 // equivalent Z-sum.
 545 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
 546 // representing S_{n,p}(x).
 547 // The calculation of Y_n uses the values from Y_{n-1}.
 548 void fill_Yn(int n, const cln::float_format_t& prec)
 549 {
 550         const int initsize = ynlength;
 551         //const int initsize = initsize_Yn;
 552         cln::cl_N one = cln::cl_float(1, prec);
 553
 554         if (n) {
 555                 std::vector<cln::cl_N> buf(initsize);
 556                 std::vector<cln::cl_N>::iterator it = buf.begin();
 557                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
 558                 *it = (*itprev) / cln::cl_N(n+1) * one;
 559                 it++;
 560                 itprev++;
 561                 // sums with an index smaller than the depth are zero and need not to be calculated.
 562                 // calculation starts with depth, which is n+2)
 563                 for (int i=n+2; i<=initsize+n; i++) {
 564                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
 565                         it++;
 566                         itprev++;
 567                 }
 568                 Yn.push_back(buf);
 569         } else {
 570                 std::vector<cln::cl_N> buf(initsize);
 571                 std::vector<cln::cl_N>::iterator it = buf.begin();
 572                 *it = 1 * one;
 573                 it++;
 574                 for (int i=2; i<=initsize; i++) {
 575                         *it = *(it-1) + 1 / cln::cl_N(i) * one;
 576                         it++;
 577                 }
 578                 Yn.push_back(buf);
 579         }
 580         ynsize++;
 581 }
 582
 583
 584 // make Yn longer ...
 585 void make_Yn_longer(int newsize, const cln::float_format_t& prec)
 586 {
 587
 588         cln::cl_N one = cln::cl_float(1, prec);
 589
 590         Yn[0].resize(newsize);
 591         std::vector<cln::cl_N>::iterator it = Yn[0].begin();
 592         it += ynlength;
 593         for (int i=ynlength+1; i<=newsize; i++) {
 594                 *it = *(it-1) + 1 / cln::cl_N(i) * one;
 595                 it++;
 596         }
 597
 598         for (int n=1; n<ynsize; n++) {
 599                 Yn[n].resize(newsize);
 600                 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
 601                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
 602                 it += ynlength;
 603                 itprev += ynlength;
 604                 for (int i=ynlength+n+1; i<=newsize+n; i++) {
 605                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
 606                         it++;
 607                         itprev++;
 608                 }
 609         }
 610
 611         ynlength = newsize;
 612 }
 613
 614
 615 // helper function for S(n,p,x)
 616 // [Kol] (7.2)
 617 cln::cl_N C(int n, int p)
 618 {
 619         cln::cl_N result;
 620
 621         for (int k=0; k<p; k++) {
 622                 for (int j=0; j<=(n+k-1)/2; j++) {
 623                         if (k == 0) {
 624                                 if (n & 1) {
 625                                         if (j & 1) {
 626                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 627                                         }
 628                                         else {
 629                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 630                                         }
 631                                 }
 632                         }
 633                         else {
 634                                 if (k & 1) {
 635                                         if (j & 1) {
 636                                                 result = result + cln::factorial(n+k-1)
 637                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 638                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 639                                         }
 640                                         else {
 641                                                 result = result - cln::factorial(n+k-1)
 642                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 643                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 644                                         }
 645                                 }
 646                                 else {
 647                                         if (j & 1) {
 648                                                 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()
 649                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 650                                         }
 651                                         else {
 652                                                 result = result + cln::factorial(n+k-1)
 653                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 654                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 655                                         }
 656                                 }
 657                         }
 658                 }
 659         }
 660         int np = n+p;
 661         if ((np-1) & 1) {
 662                 if (((np)/2+n) & 1) {
 663                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 664                 }
 665                 else {
 666                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 667                 }
 668         }
 669
 670         return result;
 671 }
 672
 673
 674 // helper function for S(n,p,x)
 675 // [Kol] remark to (9.1)
 676 cln::cl_N a_k(int k)
 677 {
 678         cln::cl_N result;
 679
 680         if (k == 0) {
 681                 return 1;
 682         }
 683
 684         result = result;
 685         for (int m=2; m<=k; m++) {
 686                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
 687         }
 688
 689         return -result / k;
 690 }
 691
 692
 693 // helper function for S(n,p,x)
 694 // [Kol] remark to (9.1)
 695 cln::cl_N b_k(int k)
 696 {
 697         cln::cl_N result;
 698
 699         if (k == 0) {
 700                 return 1;
 701         }
 702
 703         result = result;
 704         for (int m=2; m<=k; m++) {
 705                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
 706         }
 707
 708         return result / k;
 709 }
 710
 711
 712 // helper function for S(n,p,x)
 713 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 714 {
 715         if (p==1) {
 716                 return Li_projection(n+1, x, prec);
 717         }
 718
 719         // check if precalculated values are sufficient
 720         if (p > ynsize+1) {
 721                 for (int i=ynsize; i<p-1; i++) {
 722                         fill_Yn(i, prec);
 723                 }
 724         }
 725
 726         // should be done otherwise
 727         cln::cl_N xf = x * cln::cl_float(1, prec);
 728
 729         cln::cl_N res;
 730         cln::cl_N resbuf;
 731         cln::cl_N factor = cln::expt(xf, p);
 732         int i = p;
 733         do {
 734                 resbuf = res;
 735                 if (i-p >= ynlength) {
 736                         // make Yn longer
 737                         make_Yn_longer(ynlength*2, prec);
 738                 }
 739                 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? ...
 740                 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
 741                 factor = factor * xf;
 742                 i++;
 743         } while (res != resbuf);
 744
 745         return res;
 746 }
 747
 748
 749 // helper function for S(n,p,x)
 750 cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 751 {
 752         // [Kol] (5.3)
 753         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
 754
 755                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
 756                         * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
 757
 758                 for (int s=0; s<n; s++) {
 759                         cln::cl_N res2;
 760                         for (int r=0; r<p; r++) {
 761                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
 762                                         * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
 763                         }
 764                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
 765                 }
 766
 767                 return result;
 768         }
 769
 770         return S_do_sum(n, p, x, prec);
 771 }
 772
 773
 774 // helper function for S(n,p,x)
 775 numeric S_num(int n, int p, const numeric& x)
 776 {
 777         if (x == 1) {
 778                 if (n == 1) {
 779                     // [Kol] (2.22) with (2.21)
 780                         return cln::zeta(p+1);
 781                 }
 782
 783                 if (p == 1) {
 784                     // [Kol] (2.22)
 785                         return cln::zeta(n+1);
 786                 }
 787
 788                 // [Kol] (9.1)
 789                 cln::cl_N result;
 790                 for (int nu=0; nu<n; nu++) {
 791                         for (int rho=0; rho<=p; rho++) {
 792                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
 793                                         * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
 794                         }
 795                 }
 796                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
 797
 798                 return result;
 799         }
 800         else if (x == -1) {
 801                 // [Kol] (2.22)
 802                 if (p == 1) {
 803                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
 804                 }
 805 //              throw std::runtime_error("don't know how to evaluate this function!");
 806         }
 807
 808         // what is the desired float format?
 809         // first guess: default format
 810         cln::float_format_t prec = cln::default_float_format;
 811         const cln::cl_N value = x.to_cl_N();
 812         // second guess: the argument's format
 813         if (!x.real().is_rational())
 814                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
 815         else if (!x.imag().is_rational())
 816                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 817
 818
 819         // [Kol] (5.3)
 820         if (cln::realpart(value) < -0.5) {
 821
 822                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
 823                         * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
 824
 825                 for (int s=0; s<n; s++) {
 826                         cln::cl_N res2;
 827                         for (int r=0; r<p; r++) {
 828                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
 829                                         * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
 830                         }
 831                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
 832                 }
 833
 834                 return result;
 835
 836         }
 837         // [Kol] (5.12)
 838         if (cln::abs(value) > 1) {
 839
 840                 cln::cl_N result;
 841
 842                 for (int s=0; s<p; s++) {
 843                         for (int r=0; r<=s; r++) {
 844                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
 845                                         / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
 846                                         * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
 847                         }
 848                 }
 849                 result = result * cln::expt(cln::cl_I(-1),n);
 850
 851                 cln::cl_N res2;
 852                 for (int r=0; r<n; r++) {
 853                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
 854                 }
 855                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
 856
 857                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
 858
 859                 return result;
 860         }
 861         else {
 862                 return S_projection(n, p, value, prec);
 863         }
 864 }
 865
 866
 867 } // end of anonymous namespace
 868
 869
 870 //////////////////////////////////////////////////////////////////////
 871 //
 872 // Nielsen's generalized polylogarithm  S(n,p,x)
 873 //
 874 // GiNaC function
 875 //
 876 //////////////////////////////////////////////////////////////////////
 877
 878
 879 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
 880 {
 881         if (x2 == 1) {
 882                 return Li(x1+1,x3);
 883         }
 884         if (x3.info(info_flags::numeric) && (!x3.info(info_flags::crational)) &&
 885                         x1.info(info_flags::posint) && x2.info(info_flags::posint)) {
 886                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
 887         }
 888         return S(x1,x2,x3).hold();
 889 }
 890
 891
 892 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
 893 {
 894         if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
 895                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
 896         }
 897         return S(x1,x2,x3).hold();
 898 }
 899
 900
 901 static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options)
 902 {
 903         epvector seq;
 904         seq.push_back(expair(S(x1,x2,x3), 0));
 905         return pseries(rel,seq);
 906 }
 907
 908
 909 static ex S_deriv(const ex& x1, const ex& x2, const ex& x3, unsigned deriv_param)
 910 {
 911         GINAC_ASSERT(deriv_param < 3);
 912         if (deriv_param < 2) {
 913                 return _ex0;
 914         }
 915         if (x1 > 0) {
 916                 return S(x1-1, x2, x3) / x3;
 917         } else {
 918                 return S(x1, x2-1, x3) / (1-x3);
 919         }
 920 }
 921
 922
 923 REGISTER_FUNCTION(S,
 924                 eval_func(S_eval).
 925                 evalf_func(S_evalf).
 926                 do_not_evalf_params().
 927                 series_func(S_series).
 928                 derivative_func(S_deriv));
 929
 930
 931 //////////////////////////////////////////////////////////////////////
 932 //
 933 // Harmonic polylogarithm  H(m,x)
 934 //
 935 // helper function
 936 //
 937 //////////////////////////////////////////////////////////////////////
 938
 939
 940 // anonymous namespace for helper functions
 941 namespace {
 942
 943
 944 // forward declaration
 945 ex convert_from_RV(const lst& parameterlst, const ex& arg);
 946
 947
 948 // multiplies an one-dimensional H with another H
 949 // [ReV] (18)
 950 ex trafo_H_mult(const ex& h1, const ex& h2)
 951 {
 952         ex res;
 953         ex hshort;
 954         lst hlong;
 955         ex h1nops = h1.op(0).nops();
 956         ex h2nops = h2.op(0).nops();
 957         if (h1nops > 1) {
 958                 hshort = h2.op(0).op(0);
 959                 hlong = ex_to<lst>(h1.op(0));
 960         } else {
 961                 hshort = h1.op(0).op(0);
 962                 if (h2nops > 1) {
 963                         hlong = ex_to<lst>(h2.op(0));
 964                 } else {
 965                         hlong = h2.op(0).op(0);
 966                 }
 967         }
 968         for (int i=0; i<=hlong.nops(); i++) {
 969                 lst newparameter;
 970                 int j=0;
 971                 for (; j<i; j++) {
 972                         newparameter.append(hlong[j]);
 973                 }
 974                 newparameter.append(hshort);
 975                 for (; j<hlong.nops(); j++) {
 976                         newparameter.append(hlong[j]);
 977                 }
 978                 res += H(newparameter, h1.op(1)).hold();
 979         }
 980         return res;
 981 }
 982
 983
 984 // applies trafo_H_mult recursively on expressions
 985 struct map_trafo_H_mult : public map_function
 986 {
 987         ex operator()(const ex& e)
 988         {
 989                 if (is_a<add>(e)) {
 990                         return e.map(*this);
 991                 }
 992
 993                 if (is_a<mul>(e)) {
 994
 995                         ex result = 1;
 996                         ex firstH;
 997                         lst Hlst;
 998                         for (int pos=0; pos<e.nops(); pos++) {
 999                                 if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
1000                                         std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
1001                                         if (name == "H") {
1002                                                 for (ex i=0; i<e.op(pos).op(1); i++) {
1003                                                         Hlst.append(e.op(pos).op(0));
1004                                                 }
1005                                                 continue;
1006                                         }
1007                                 } else if (is_a<function>(e.op(pos))) {
1008                                         std::string name = ex_to<function>(e.op(pos)).get_name();
1009                                         if (name == "H") {
1010                                                 if (e.op(pos).op(0).nops() > 1) {
1011                                                         firstH = e.op(pos);
1012                                                 } else {
1013                                                         Hlst.append(e.op(pos));
1014                                                 }
1015                                                 continue;
1016                                         }
1017                                 }
1018                                 result *= e.op(pos);
1019                         }
1020                         if (firstH == 0) {
1021                                 if (Hlst.nops() > 0) {
1022                                         firstH = Hlst[Hlst.nops()-1];
1023                                         Hlst.remove_last();
1024                                 } else {
1025                                         return e;
1026                                 }
1027                         }
1028
1029                         if (Hlst.nops() > 0) {
1030                                 ex buffer = trafo_H_mult(firstH, Hlst.op(0));
1031                                 result *= buffer;
1032                                 for (int i=1; i<Hlst.nops(); i++) {
1033                                         result *= Hlst.op(i);
1034                                 }
1035                                 result = result.expand();
1036                                 map_trafo_H_mult recursion;
1037                                 return recursion(result);
1038                         } else {
1039                                 return e;
1040                         }
1041
1042                 }
1043                 return e;
1044         }
1045 };
1046
1047
1048 // do integration [ReV] (49)
1049 // put parameter 1 in front of existing parameters
1050 ex trafo_H_prepend_one(const ex& e, const ex& arg)
1051 {
1052         ex h;
1053         std::string name;
1054         if (is_a<function>(e)) {
1055                 name = ex_to<function>(e).get_name();
1056         }
1057         if (name == "H") {
1058                 h = e;
1059         } else {
1060                 for (int i=0; i<e.nops(); i++) {
1061                         if (is_a<function>(e.op(i))) {
1062                                 std::string name = ex_to<function>(e.op(i)).get_name();
1063                                 if (name == "H") {
1064                                         h = e.op(i);
1065                                 }
1066                         }
1067                 }
1068         }
1069         if (h != 0) {
1070                 lst newparameter = ex_to<lst>(h.op(0));
1071                 newparameter.prepend(1);
1072                 return e.subs(h == H(newparameter, h.op(1)).hold());
1073         } else {
1074                 return e * H(lst(1),1-arg).hold();
1075         }
1076 }
1077
1078
1079 // do integration [ReV] (55)
1080 // put parameter 0 in front of existing parameters
1081 ex trafo_H_prepend_zero(const ex& e, const ex& arg)
1082 {
1083         ex h;
1084         std::string name;
1085         if (is_a<function>(e)) {
1086                 name = ex_to<function>(e).get_name();
1087         }
1088         if (name == "H") {
1089                 h = e;
1090         } else {
1091                 for (int i=0; i<e.nops(); i++) {
1092                         if (is_a<function>(e.op(i))) {
1093                                 std::string name = ex_to<function>(e.op(i)).get_name();
1094                                 if (name == "H") {
1095                                         h = e.op(i);
1096                                 }
1097                         }
1098                 }
1099         }
1100         if (h != 0) {
1101                 lst newparameter = ex_to<lst>(h.op(0));
1102                 newparameter.prepend(0);
1103                 ex addzeta = convert_from_RV(newparameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1104                 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
1105         } else {
1106                 return e * (-H(lst(0),1/arg).hold());
1107         }
1108 }
1109
1110
1111 // do x -> 1-x transformation
1112 struct map_trafo_H_1mx : public map_function
1113 {
1114         ex operator()(const ex& e)
1115         {
1116                 if (is_a<add>(e) || is_a<mul>(e)) {
1117                         return e.map(*this);
1118                 }
1119
1120                 if (is_a<function>(e)) {
1121                         std::string name = ex_to<function>(e).get_name();
1122                         if (name == "H") {
1123
1124                                 lst parameter = ex_to<lst>(e.op(0));
1125                                 ex arg = e.op(1);
1126
1127                                 // if all parameters are either zero or one return the transformed function
1128                                 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1129                                         lst newparameter;
1130                                         for (int i=parameter.nops(); i>0; i--) {
1131                                                 newparameter.append(0);
1132                                         }
1133                                         return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1134                                 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1135                                         lst newparameter;
1136                                         for (int i=parameter.nops(); i>0; i--) {
1137                                                 newparameter.append(1);
1138                                         }
1139                                         return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
1140                                 }
1141
1142                                 lst newparameter = parameter;
1143                                 newparameter.remove_first();
1144
1145                                 if (parameter.op(0) == 0) {
1146
1147                                         // leading zero
1148                                         ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1149                                         map_trafo_H_1mx recursion;
1150                                         ex buffer = recursion(H(newparameter, arg).hold());
1151                                         if (is_a<add>(buffer)) {
1152                                                 for (int i=0; i<buffer.nops(); i++) {
1153                                                         res -= trafo_H_prepend_one(buffer.op(i), arg);
1154                                                 }
1155                                         } else {
1156                                                 res -= trafo_H_prepend_one(buffer, arg);
1157                                         }
1158                                         return res;
1159
1160                                 } else {
1161
1162                                         // leading one
1163                                         map_trafo_H_1mx recursion;
1164                                         map_trafo_H_mult unify;
1165                                         ex res;
1166                                         int firstzero = 0;
1167                                         while (parameter.op(firstzero) == 1) {
1168                                                 firstzero++;
1169                                         }
1170                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1171                                                 lst newparameter;
1172                                                 int j=0;
1173                                                 for (; j<=i; j++) {
1174                                                         newparameter.append(parameter[j+1]);
1175                                                 }
1176                                                 newparameter.append(1);
1177                                                 for (; j<parameter.nops()-1; j++) {
1178                                                         newparameter.append(parameter[j+1]);
1179                                                 }
1180                                                 res -= H(newparameter, arg).hold();
1181                                         }
1182                                         return (unify((-H(lst(0), 1-arg).hold() * recursion(H(newparameter, arg).hold())).expand()) +
1183                                                         recursion(res)) / firstzero;
1184
1185                                 }
1186
1187                         }
1188                 }
1189                 return e;
1190         }
1191 };
1192
1193
1194 // do x -> 1/x transformation
1195 struct map_trafo_H_1overx : public map_function
1196 {
1197         ex operator()(const ex& e)
1198         {
1199                 if (is_a<add>(e) || is_a<mul>(e)) {
1200                         return e.map(*this);
1201                 }
1202
1203                 if (is_a<function>(e)) {
1204                         std::string name = ex_to<function>(e).get_name();
1205                         if (name == "H") {
1206
1207                                 lst parameter = ex_to<lst>(e.op(0));
1208                                 ex arg = e.op(1);
1209
1210                                 // if all parameters are either zero or one return the transformed function
1211                                 if (find(parameter.begin(), parameter.end(), 0) == parameter.end()) {
1212                                         map_trafo_H_mult unify;
1213                                         return unify((pow(H(lst(1),1/arg).hold() + H(lst(0),1/arg).hold() - I*Pi, parameter.nops()) /
1214                                                                 factorial(parameter.nops())).expand());
1215                                 } else if (find(parameter.begin(), parameter.end(), 1) == parameter.end()) {
1216                                         return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
1217                                 }
1218
1219                                 lst newparameter = parameter;
1220                                 newparameter.remove_first();
1221
1222                                 if (parameter.op(0) == 0) {
1223
1224                                         // leading zero
1225                                         ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
1226                                         map_trafo_H_1overx recursion;
1227                                         ex buffer = recursion(H(newparameter, arg).hold());
1228                                         if (is_a<add>(buffer)) {
1229                                                 for (int i=0; i<buffer.nops(); i++) {
1230                                                         res += trafo_H_prepend_zero(buffer.op(i), arg);
1231                                                 }
1232                                         } else {
1233                                                 res += trafo_H_prepend_zero(buffer, arg);
1234                                         }
1235                                         return res;
1236
1237                                 } else {
1238
1239                                         // leading one
1240                                         map_trafo_H_1overx recursion;
1241                                         map_trafo_H_mult unify;
1242                                         ex res = H(lst(1), arg).hold() * H(newparameter, arg).hold();
1243                                         int firstzero = 0;
1244                                         while (parameter.op(firstzero) == 1) {
1245                                                 firstzero++;
1246                                         }
1247                                         for (int i=firstzero-1; i<parameter.nops()-1; i++) {
1248                                                 lst newparameter;
1249                                                 int j=0;
1250                                                 for (; j<=i; j++) {
1251                                                         newparameter.append(parameter[j+1]);
1252                                                 }
1253                                                 newparameter.append(1);
1254                                                 for (; j<parameter.nops()-1; j++) {
1255                                                         newparameter.append(parameter[j+1]);
1256                                                 }
1257                                                 res -= H(newparameter, arg).hold();
1258                                         }
1259                                         res = recursion(res).expand() / firstzero;
1260                                         return unify(res);
1261
1262                                 }
1263
1264                         }
1265                 }
1266                 return e;
1267         }
1268 };
1269
1270
1271 // remove trailing zeros from H-parameters
1272 struct map_trafo_H_reduce_trailing_zeros : 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                                 ex arg = e.op(1);
1289                                 if (parameter.op(parameter.nops()-1) == 0) {
1290
1291                                         //
1292                                         if (parameter.nops() == 1) {
1293                                                 return log(arg);
1294                                         }
1295
1296                                         //
1297                                         lst::const_iterator it = parameter.begin();
1298                                         while ((it != parameter.end()) && (*it == 0)) {
1299                                                 it++;
1300                                         }
1301                                         if (it == parameter.end()) {
1302                                                 return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
1303                                         }
1304
1305                                         //
1306                                         parameter.remove_last();
1307                                         int lastentry = parameter.nops();
1308                                         while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
1309                                                 lastentry--;
1310                                         }
1311
1312                                         //
1313                                         ex result = log(arg) * H(parameter,arg).hold();
1314                                         for (ex i=0; i<lastentry; i++) {
1315                                                 parameter[i]++;
1316                                                 result -= (parameter[i]-1) * H(parameter, arg).hold();
1317                                                 parameter[i]--;
1318                                         }
1319
1320                                         if (lastentry < parameter.nops()) {
1321                                                 result = result / (parameter.nops()-lastentry+1);
1322                                                 return result.map(*this);
1323                                         } else {
1324                                                 return result;
1325                                         }
1326                                 }
1327                         }
1328                 }
1329                 return e;
1330         }
1331 };
1332
1333
1334 // recursively call convert_from_RV on expression
1335 struct map_trafo_H_convert : public map_function
1336 {
1337         ex operator()(const ex& e)
1338         {
1339                 if (is_a<add>(e) || is_a<mul>(e) || is_a<power>(e)) {
1340                         return e.map(*this);
1341                 }
1342                 if (is_a<function>(e)) {
1343                         std::string name = ex_to<function>(e).get_name();
1344                         if (name == "H") {
1345                                 lst parameter = ex_to<lst>(e.op(0));
1346                                 ex arg = e.op(1);
1347                                 return convert_from_RV(parameter, arg);
1348                         }
1349                 }
1350                 return e;
1351         }
1352 };
1353
1354
1355 // translate notation from nested sums to Remiddi/Vermaseren
1356 lst convert_to_RV(const lst& o)
1357 {
1358         lst res;
1359         for (lst::const_iterator it = o.begin(); it != o.end(); it++) {
1360                 for (ex i=0; i<(*it)-1; i++) {
1361                         res.append(0);
1362                 }
1363                 res.append(1);
1364         }
1365         return res;
1366 }
1367
1368
1369 // translate notation from Remiddi/Vermaseren to nested sums
1370 ex convert_from_RV(const lst& parameterlst, const ex& arg)
1371 {
1372         lst newparameterlst;
1373
1374         lst::const_iterator it = parameterlst.begin();
1375         int count = 1;
1376         while (it != parameterlst.end()) {
1377                 if (*it == 0) {
1378                         count++;
1379                 } else {
1380                         newparameterlst.append((*it>0) ? count : -count);
1381                         count = 1;
1382                 }
1383                 it++;
1384         }
1385         for (int i=1; i<count; i++) {
1386                 newparameterlst.append(0);
1387         }
1388
1389         map_trafo_H_reduce_trailing_zeros filter;
1390         return filter(H(newparameterlst, arg).hold());
1391 }
1392
1393
1394 // do the actual summation.
1395 cln::cl_N H_do_sum(const std::vector<int>& s, const cln::cl_N& x)
1396 {
1397         const int j = s.size();
1398
1399         std::vector<cln::cl_N> t(j);
1400
1401         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1402         cln::cl_N factor = cln::expt(x, j) * one;
1403         cln::cl_N t0buf;
1404         int q = 0;
1405         do {
1406                 t0buf = t[0];
1407                 q++;
1408                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1409                 for (int k=j-2; k>=1; k--) {
1410                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1411                 }
1412                 t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), s[0]);
1413                 factor = factor * x;
1414         } while (t[0] != t0buf);
1415
1416         return t[0];
1417 }
1418
1419
1420 } // end of anonymous namespace
1421
1422
1423 //////////////////////////////////////////////////////////////////////
1424 //
1425 // Harmonic polylogarithm  H(m,x)
1426 //
1427 // GiNaC function
1428 //
1429 //////////////////////////////////////////////////////////////////////
1430
1431
1432 static ex H_eval(const ex& x1, const ex& x2)
1433 {
1434         if (x2 == 0) {
1435                 return 0;
1436         }
1437         if (x2 == 1) {
1438                 return zeta(x1);
1439         }
1440         if (x1.nops() == 1) {
1441                 return Li(x1.op(0), x2);
1442         }
1443         if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational))) {
1444                 return H(x1,x2).evalf();
1445         }
1446         return H(x1,x2).hold();
1447 }
1448
1449
1450 static ex H_evalf(const ex& x1, const ex& x2)
1451 {
1452         if (is_a<lst>(x1) && is_a<numeric>(x2)) {
1453                 for (int i=0; i<x1.nops(); i++) {
1454                         if (!x1.op(i).info(info_flags::posint)) {
1455                                 return H(x1,x2).hold();
1456                         }
1457                 }
1458                 if (x1.nops() < 1) {
1459                         return _ex1;
1460                 }
1461                 if (x1.nops() == 1) {
1462                         return Li(x1.op(0), x2).evalf();
1463                 }
1464                 cln::cl_N x = ex_to<numeric>(x2).to_cl_N();
1465                 if (x == 1) {
1466                         return zeta(x1).evalf();
1467                 }
1468
1469                 // choose trafo
1470                 if (cln::abs(x) > 1) {
1471                         symbol xtemp("xtemp");
1472                         map_trafo_H_1overx trafo;
1473                         ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1474                         map_trafo_H_convert converter;
1475                         res = converter(res);
1476                         return res.subs(xtemp==x2).evalf();
1477                 }
1478
1479                 // since the x->1-x transformation produces a lot of terms, it is only
1480                 // efficient for argument near one.
1481                 if (cln::realpart(x) > 0.95) {
1482                         symbol xtemp("xtemp");
1483                         map_trafo_H_1mx trafo;
1484                         ex res = trafo(H(convert_to_RV(ex_to<lst>(x1)), xtemp));
1485                         map_trafo_H_convert converter;
1486                         res = converter(res);
1487                         return res.subs(xtemp==x2).evalf();
1488                 }
1489
1490                 // no trafo -> do summation
1491                 int count = x1.nops();
1492                 std::vector<int> r(count);
1493                 for (int i=0; i<count; i++) {
1494                         r[i] = ex_to<numeric>(x1.op(i)).to_int();
1495                 }
1496
1497                 return numeric(H_do_sum(r,x));
1498         }
1499
1500         return H(x1,x2).hold();
1501 }
1502
1503
1504 static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
1505 {
1506         epvector seq;
1507         seq.push_back(expair(H(x1,x2), 0));
1508         return pseries(rel,seq);
1509 }
1510
1511
1512 static ex H_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
1513 {
1514         GINAC_ASSERT(deriv_param < 2);
1515         if (deriv_param == 0) {
1516                 return _ex0;
1517         }
1518         if (is_a<lst>(x1)) {
1519                 lst newparameter = ex_to<lst>(x1);
1520                 if (x1.op(0) == 1) {
1521                         newparameter.remove_first();
1522                         return 1/(1-x2) * H(newparameter, x2);
1523                 } else {
1524                         newparameter[0]--;
1525                         return H(newparameter, x2).hold() / x2;
1526                 }
1527         } else {
1528                 if (x1 == 1) {
1529                         return 1/(1-x2);
1530                 } else {
1531                         return H(x1-1, x2).hold() / x2;
1532                 }
1533         }
1534 }
1535
1536
1537 REGISTER_FUNCTION(H,
1538                 eval_func(H_eval).
1539                 evalf_func(H_evalf).
1540                 do_not_evalf_params().
1541                 series_func(H_series).
1542                 derivative_func(H_deriv));
1543
1544
1545 //////////////////////////////////////////////////////////////////////
1546 //
1547 // Multiple zeta values  zeta(x) and zeta(x,s)
1548 //
1549 // helper functions
1550 //
1551 //////////////////////////////////////////////////////////////////////
1552
1553
1554 // anonymous namespace for helper functions
1555 namespace {
1556
1557
1558 // parameters and data for [Cra] algorithm
1559 const cln::cl_N lambda = cln::cl_N("319/320");
1560 int L1;
1561 int L2;
1562 std::vector<std::vector<cln::cl_N> > f_kj;
1563 std::vector<cln::cl_N> crB;
1564 std::vector<std::vector<cln::cl_N> > crG;
1565 std::vector<cln::cl_N> crX;
1566
1567
1568 void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
1569 {
1570         const int size = a.size();
1571         for (int n=0; n<size; n++) {
1572                 c[n] = 0;
1573                 for (int m=0; m<=n; m++) {
1574                         c[n] = c[n] + a[m]*b[n-m];
1575                 }
1576         }
1577 }
1578
1579
1580 // [Cra] section 4
1581 void initcX(const std::vector<int>& s)
1582 {
1583         const int k = s.size();
1584
1585         crX.clear();
1586         crG.clear();
1587         crB.clear();
1588
1589         for (int i=0; i<=L2; i++) {
1590                 crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
1591         }
1592
1593         int Sm = 0;
1594         int Smp1 = 0;
1595         for (int m=0; m<k-1; m++) {
1596                 std::vector<cln::cl_N> crGbuf;
1597                 Sm = Sm + s[m];
1598                 Smp1 = Sm + s[m+1];
1599                 for (int i=0; i<=L2; i++) {
1600                         crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
1601                 }
1602                 crG.push_back(crGbuf);
1603         }
1604
1605         crX = crB;
1606
1607         for (int m=0; m<k-1; m++) {
1608                 std::vector<cln::cl_N> Xbuf;
1609                 for (int i=0; i<=L2; i++) {
1610                         Xbuf.push_back(crX[i] * crG[m][i]);
1611                 }
1612                 halfcyclic_convolute(Xbuf, crB, crX);
1613         }
1614 }
1615
1616
1617 // [Cra] section 4
1618 cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
1619 {
1620         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1621         cln::cl_N factor = cln::expt(lambda, Sqk);
1622         cln::cl_N res = factor / Sqk * crX[0] * one;
1623         cln::cl_N resbuf;
1624         int N = 0;
1625         do {
1626                 resbuf = res;
1627                 factor = factor * lambda;
1628                 N++;
1629                 res = res + crX[N] * factor / (N+Sqk);
1630         } while ((res != resbuf) || cln::zerop(crX[N]));
1631         return res;
1632 }
1633
1634
1635 // [Cra] section 4
1636 void calc_f(int maxr)
1637 {
1638         f_kj.clear();
1639         f_kj.resize(L1);
1640
1641         cln::cl_N t0, t1, t2, t3, t4;
1642         int i, j, k;
1643         std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
1644         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1645
1646         t0 = cln::exp(-lambda);
1647         t2 = 1;
1648         for (k=1; k<=L1; k++) {
1649                 t1 = k * lambda;
1650                 t2 = t0 * t2;
1651                 for (j=1; j<=maxr; j++) {
1652                         t3 = 1;
1653                         t4 = 1;
1654                         for (i=2; i<=j; i++) {
1655                                 t4 = t4 * (j-i+1);
1656                                 t3 = t1 * t3 + t4;
1657                         }
1658                         (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one);
1659                 }
1660                 it++;
1661         }
1662 }
1663
1664
1665 // [Cra] (3.1)
1666 cln::cl_N crandall_Z(const std::vector<int>& s)
1667 {
1668         const int j = s.size();
1669
1670         if (j == 1) {
1671                 cln::cl_N t0;
1672                 cln::cl_N t0buf;
1673                 int q = 0;
1674                 do {
1675                         t0buf = t0;
1676                         q++;
1677                         t0 = t0 + f_kj[q+j-2][s[0]-1];
1678                 } while (t0 != t0buf);
1679
1680                 return t0 / cln::factorial(s[0]-1);
1681         }
1682
1683         std::vector<cln::cl_N> t(j);
1684
1685         cln::cl_N t0buf;
1686         int q = 0;
1687         do {
1688                 t0buf = t[0];
1689                 q++;
1690                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1691                 for (int k=j-2; k>=1; k--) {
1692                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1693                 }
1694                 t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
1695         } while (t[0] != t0buf);
1696
1697         return t[0] / cln::factorial(s[0]-1);
1698 }
1699
1700
1701 // [Cra] (2.4)
1702 cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
1703 {
1704         std::vector<int> r = s;
1705         const int j = r.size();
1706
1707         // decide on maximal size of f_kj for crandall_Z
1708         if (Digits < 50) {
1709                 L1 = 150;
1710         } else {
1711                 L1 = Digits * 3 + j*2;
1712         }
1713
1714         // decide on maximal size of crX for crandall_Y
1715         if (Digits < 38) {
1716                 L2 = 63;
1717         } else if (Digits < 86) {
1718                 L2 = 127;
1719         } else if (Digits < 192) {
1720                 L2 = 255;
1721         } else if (Digits < 394) {
1722                 L2 = 511;
1723         } else if (Digits < 808) {
1724                 L2 = 1023;
1725         } else {
1726                 L2 = 2047;
1727         }
1728
1729         cln::cl_N res;
1730
1731         int maxr = 0;
1732         int S = 0;
1733         for (int i=0; i<j; i++) {
1734                 S += r[i];
1735                 if (r[i] > maxr) {
1736                         maxr = r[i];
1737                 }
1738         }
1739
1740         calc_f(maxr);
1741
1742         const cln::cl_N r0factorial = cln::factorial(r[0]-1);
1743
1744         std::vector<int> rz;
1745         int skp1buf;
1746         int Srun = S;
1747         for (int k=r.size()-1; k>0; k--) {
1748
1749                 rz.insert(rz.begin(), r.back());
1750                 skp1buf = rz.front();
1751                 Srun -= skp1buf;
1752                 r.pop_back();
1753
1754                 initcX(r);
1755
1756                 for (int q=0; q<skp1buf; q++) {
1757
1758                         cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
1759                         cln::cl_N pp2 = crandall_Z(rz);
1760
1761                         rz.front()--;
1762
1763                         if (q & 1) {
1764                                 res = res - pp1 * pp2 / cln::factorial(q);
1765                         } else {
1766                                 res = res + pp1 * pp2 / cln::factorial(q);
1767                         }
1768                 }
1769                 rz.front() = skp1buf;
1770         }
1771         rz.insert(rz.begin(), r.back());
1772
1773         initcX(rz);
1774
1775         res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
1776
1777         return res;
1778 }
1779
1780
1781 cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
1782 {
1783         const int j = r.size();
1784
1785         // buffer for subsums
1786         std::vector<cln::cl_N> t(j);
1787         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1788
1789         cln::cl_N t0buf;
1790         int q = 0;
1791         do {
1792                 t0buf = t[0];
1793                 q++;
1794                 t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
1795                 for (int k=j-2; k>=0; k--) {
1796                         t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
1797                 }
1798         } while (t[0] != t0buf);
1799
1800         return t[0];
1801 }
1802
1803
1804 // does Hoelder convolution. see [BBB] (7.0)
1805 cln::cl_N zeta_do_Hoelder_convolution(const std::vector<int>& m_, const std::vector<int>& s_)
1806 {
1807         // prepare parameters
1808         // holds Li arguments in [BBB] notation
1809         std::vector<int> s = s_;
1810         std::vector<int> m_p = m_;
1811         std::vector<int> m_q;
1812         // holds Li arguments in nested sums notation
1813         std::vector<cln::cl_N> s_p(s.size(), cln::cl_N(1));
1814         s_p[0] = s_p[0] * cln::cl_N("1/2");
1815         // convert notations
1816         int sig = 1;
1817         for (int i=0; i<s_.size(); i++) {
1818                 if (s_[i] < 0) {
1819                         sig = -sig;
1820                         s_p[i] = -s_p[i];
1821                 }
1822                 s[i] = sig * std::abs(s[i]);
1823         }
1824         std::vector<cln::cl_N> s_q;
1825         cln::cl_N signum = 1;
1826
1827         // first term
1828         cln::cl_N res = multipleLi_do_sum(m_p, s_p);
1829
1830         // middle terms
1831         do {
1832
1833                 // change parameters
1834                 if (s.front() > 0) {
1835                         if (m_p.front() == 1) {
1836                                 m_p.erase(m_p.begin());
1837                                 s_p.erase(s_p.begin());
1838                                 if (s_p.size() > 0) {
1839                                         s_p.front() = s_p.front() * cln::cl_N("1/2");
1840                                 }
1841                                 s.erase(s.begin());
1842                                 m_q.front()++;
1843                         } else {
1844                                 m_p.front()--;
1845                                 m_q.insert(m_q.begin(), 1);
1846                                 if (s_q.size() > 0) {
1847                                         s_q.front() = s_q.front() * 2;
1848                                 }
1849                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
1850                         }
1851                 } else {
1852                         if (m_p.front() == 1) {
1853                                 m_p.erase(m_p.begin());
1854                                 s_p.erase(s_p.begin());
1855                                 if (s_p.size() > 0) {
1856                                         s_p.front() = s_p.front() * cln::cl_N("1/2");
1857                                 }
1858                                 s.erase(s.begin());
1859                                 for (int i=0; i<s.size(); i++) {
1860                                         s_p[i] = -s_p[i];
1861                                         if (s[i] > 0) break;
1862
1863                                 }
1864                                 m_q.insert(m_q.begin(), 1);
1865                                 if (s_q.size() > 0) {
1866                                         s_q.front() = s_q.front() * 4;
1867                                 }
1868                                 s_q.insert(s_q.begin(), cln::cl_N("1/4"));
1869                                 signum = -signum;
1870                         } else {
1871                                 m_p.front()--;
1872                                 m_q.insert(m_q.begin(), 1);
1873                                 if (s_q.size() > 0) {
1874                                         s_q.front() = s_q.front() * 2;
1875                                 }
1876                                 s_q.insert(s_q.begin(), cln::cl_N("1/2"));
1877                         }
1878                 }
1879
1880                 // exiting the loop
1881                 if (m_p.size() == 0) break;
1882
1883                 res = res + signum * multipleLi_do_sum(m_p, s_p) * multipleLi_do_sum(m_q, s_q);
1884
1885         } while (true);
1886
1887         // last term
1888         res = res + signum * multipleLi_do_sum(m_q, s_q);
1889
1890         return res;
1891 }
1892
1893
1894 } // end of anonymous namespace
1895
1896
1897 //////////////////////////////////////////////////////////////////////
1898 //
1899 // Multiple zeta values  zeta(x)
1900 //
1901 // GiNaC function
1902 //
1903 //////////////////////////////////////////////////////////////////////
1904
1905
1906 static ex zeta1_evalf(const ex& x)
1907 {
1908         if (is_exactly_a<lst>(x) && (x.nops()>1)) {
1909
1910                 // multiple zeta value
1911                 const int count = x.nops();
1912                 const lst& xlst = ex_to<lst>(x);
1913                 std::vector<int> r(count);
1914
1915                 // check parameters and convert them
1916                 lst::const_iterator it1 = xlst.begin();
1917                 std::vector<int>::iterator it2 = r.begin();
1918                 do {
1919                         if (!(*it1).info(info_flags::posint)) {
1920                                 return zeta(x).hold();
1921                         }
1922                         *it2 = ex_to<numeric>(*it1).to_int();
1923                         it1++;
1924                         it2++;
1925                 } while (it2 != r.end());
1926
1927                 // check for divergence
1928                 if (r[0] == 1) {
1929                         return zeta(x).hold();
1930                 }
1931
1932                 // decide on summation algorithm
1933                 // this is still a bit clumsy
1934                 int limit = (Digits>17) ? 10 : 6;
1935                 if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) {
1936                         return numeric(zeta_do_sum_Crandall(r));
1937                 } else {
1938                         return numeric(zeta_do_sum_simple(r));
1939                 }
1940         }
1941
1942         // single zeta value
1943         if (is_exactly_a<numeric>(x) && (x != 1)) {
1944                 try {
1945                         return zeta(ex_to<numeric>(x));
1946                 } catch (const dunno &e) { }
1947         }
1948
1949         return zeta(x).hold();
1950 }
1951
1952
1953 static ex zeta1_eval(const ex& x)
1954 {
1955         if (is_exactly_a<lst>(x)) {
1956                 if (x.nops() == 1) {
1957                         return zeta(x.op(0));
1958                 }
1959                 return zeta(x).hold();
1960         }
1961
1962         if (x.info(info_flags::numeric)) {
1963                 const numeric& y = ex_to<numeric>(x);
1964                 // trap integer arguments:
1965                 if (y.is_integer()) {
1966                         if (y.is_zero()) {
1967                                 return _ex_1_2;
1968                         }
1969                         if (y.is_equal(_num1)) {
1970                                 return zeta(x).hold();
1971                         }
1972                         if (y.info(info_flags::posint)) {
1973                                 if (y.info(info_flags::odd)) {
1974                                         return zeta(x).hold();
1975                                 } else {
1976                                         return abs(bernoulli(y)) * pow(Pi, y) * pow(_num2, y-_num1) / factorial(y);
1977                                 }
1978                         } else {
1979                                 if (y.info(info_flags::odd)) {
1980                                         return -bernoulli(_num1-y) / (_num1-y);
1981                                 } else {
1982                                         return _ex0;
1983                                 }
1984                         }
1985                 }
1986                 // zeta(float)
1987                 if (y.info(info_flags::numeric) && !y.info(info_flags::crational))
1988                         return zeta1_evalf(x);
1989         }
1990         return zeta(x).hold();
1991 }
1992
1993
1994 static ex zeta1_deriv(const ex& x, unsigned deriv_param)
1995 {
1996         GINAC_ASSERT(deriv_param==0);
1997
1998         if (is_exactly_a<lst>(x)) {
1999                 return _ex0;
2000         } else {
2001                 return zeta(_ex1, x);
2002         }
2003 }
2004
2005
2006 unsigned zeta1_SERIAL::serial =
2007                         function::register_new(function_options("zeta").
2008                                                 eval_func(zeta1_eval).
2009                                                 evalf_func(zeta1_evalf).
2010                                                 do_not_evalf_params().
2011                                                 derivative_func(zeta1_deriv).
2012                                                 latex_name("\\zeta").
2013                                                 overloaded(2));
2014
2015
2016 //////////////////////////////////////////////////////////////////////
2017 //
2018 // Alternating Euler sum  zeta(x,s)
2019 //
2020 // GiNaC function
2021 //
2022 //////////////////////////////////////////////////////////////////////
2023
2024
2025 static ex zeta2_evalf(const ex& x, const ex& s)
2026 {
2027         if (is_exactly_a<lst>(x)) {
2028
2029                 // alternating Euler sum
2030                 const int count = x.nops();
2031                 const lst& xlst = ex_to<lst>(x);
2032                 const lst& slst = ex_to<lst>(s);
2033                 std::vector<int> xi(count);
2034                 std::vector<int> si(count);
2035
2036                 // check parameters and convert them
2037                 lst::const_iterator it_xread = xlst.begin();
2038                 lst::const_iterator it_sread = slst.begin();
2039                 std::vector<int>::iterator it_xwrite = xi.begin();
2040                 std::vector<int>::iterator it_swrite = si.begin();
2041                 do {
2042                         if (!(*it_xread).info(info_flags::posint)) {
2043                                 return zeta(x, s).hold();
2044                         }
2045                         *it_xwrite = ex_to<numeric>(*it_xread).to_int();
2046                         if (*it_sread > 0) {
2047                                 *it_swrite = 1;
2048                         } else {
2049                                 *it_swrite = -1;
2050                         }
2051                         it_xread++;
2052                         it_sread++;
2053                         it_xwrite++;
2054                         it_swrite++;
2055                 } while (it_xwrite != xi.end());
2056
2057                 // check for divergence
2058                 if ((xi[0] == 1) && (si[0] == 1)) {
2059                         return zeta(x, s).hold();
2060                 }
2061
2062                 // use Hoelder convolution
2063                 return numeric(zeta_do_Hoelder_convolution(xi, si));
2064         }
2065
2066         return zeta(x, s).hold();
2067 }
2068
2069
2070 static ex zeta2_eval(const ex& x, const ex& s)
2071 {
2072         if (is_exactly_a<lst>(s)) {
2073                 const lst& l = ex_to<lst>(s);
2074                 lst::const_iterator it = l.begin();
2075                 while (it != l.end()) {
2076                         if ((*it).info(info_flags::negative)) {
2077                                 return zeta(x, s).hold();
2078                         }
2079                         it++;
2080                 }
2081                 return zeta(x);
2082         } else {
2083                 if (s.info(info_flags::positive)) {
2084                         return zeta(x);
2085                 }
2086         }
2087
2088         return zeta(x, s).hold();
2089 }
2090
2091
2092 static ex zeta2_deriv(const ex& x, const ex& s, unsigned deriv_param)
2093 {
2094         GINAC_ASSERT(deriv_param==0);
2095
2096         if (is_exactly_a<lst>(x)) {
2097                 return _ex0;
2098         } else {
2099                 if ((is_exactly_a<lst>(s) && (s.op(0) > 0)) || (s > 0)) {
2100                         return zeta(_ex1, x);
2101                 }
2102                 return _ex0;
2103         }
2104 }
2105
2106
2107 unsigned zeta2_SERIAL::serial =
2108                         function::register_new(function_options("zeta").
2109                                                 eval_func(zeta2_eval).
2110                                                 evalf_func(zeta2_evalf).
2111                                                 do_not_evalf_params().
2112                                                 derivative_func(zeta2_deriv).
2113                                                 latex_name("\\zeta").
2114                                                 overloaded(2));
2115
2116
2117 } // namespace GiNaC
2118