ginac/inifcns_nstdsums.cpp

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