ginac/inifcns_nstdsums.cpp

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