ginac/inifcns_nstdsums.cpp

   1 /** @file inifcns_nstdsums.cpp
   2  *
   3  *  Implementation of some special functions that have a representation as nested sums.
   4  *  The functions are:
   5  *    classical polylogarithm              Li(n,x)
   6  *    multiple polylogarithm               Li(lst(n_1,...,n_k),lst(x_1,...,x_k))
   7  *    nielsen's generalized polylogarithm  S(n,p,x)
   8  *    harmonic polylogarithm               H(lst(m_1,...,m_k),x)
   9  *    multiple zeta value                  mZeta(lst(m_1,...,m_k))
  10  *
  11  *  Some remarks:
  12  *    - All formulae used can be looked up in the following publications:
  13  *      [Kol] Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
  14  *      [Cra] Fast Evaluation of Multiple Zeta Sums, R.E.Crandall, Math.Comp. 67 (1998), pp. 1163-1172.
  15  *    - Classical polylogarithms (Li) and nielsen's generalized polylogarithms (S) can be numerically
  16  *      evaluated in the whole complex plane.
  17  *    - The calculation of classical polylogarithms is speed up by using Euler-Maclaurin summation (EuMac).
  18  *    - The remaining functions can only be numerically evaluated with arguments lying in the unit sphere
  19  *      at the moment.
  20  *    - The functions have no series expansion. To do it, you have to convert these functions
  21  *      into the appropriate objects from the nestedsums library, do the expansion and convert the
  22  *      result back.
  23  *    - Numerical testing of this implementation has been performed by doing a comparison of results
  24  *      between this software and the commercial M.......... 4.1.
  25  *
  26  */
  27
  28 /*
  29  *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
  30  *
  31  *  This program is free software; you can redistribute it and/or modify
  32  *  it under the terms of the GNU General Public License as published by
  33  *  the Free Software Foundation; either version 2 of the License, or
  34  *  (at your option) any later version.
  35  *
  36  *  This program is distributed in the hope that it will be useful,
  37  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  38  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  39  *  GNU General Public License for more details.
  40  *
  41  *  You should have received a copy of the GNU General Public License
  42  *  along with this program; if not, write to the Free Software
  43  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  44  */
  45
  46 #include <stdexcept>
  47 #include <vector>
  48 #include <cln/cln.h>
  49
  50 #include "inifcns.h"
  51 #include "lst.h"
  52 #include "numeric.h"
  53 #include "operators.h"
  54 #include "relational.h"
  55 #include "pseries.h"
  56
  57
  58 namespace GiNaC {
  59
  60
  61 //////////////////////////////////////////////////////////////////////
  62 //
  63 // Classical polylogarithm  Li
  64 //
  65 // helper functions
  66 //
  67 //////////////////////////////////////////////////////////////////////
  68
  69
  70 namespace {
  71 // lookup table for Euler-Maclaurin optimization
  72 // see fill_Xn()
  73 std::vector<std::vector<cln::cl_N> > Xn;
  74 int xnsize = 0;
  75 }
  76
  77
  78 // This function calculates the X_n. The X_n are needed for the Euler-Maclaurin summation (EMS) of
  79 // classical polylogarithms.
  80 // With EMS the polylogs can be calculated as follows:
  81 //   Li_p (x)  =  \sum_{n=0}^\infty X_{p-2}(n) u^{n+1}/(n+1)! with  u = -log(1-x)
  82 //   X_0(n) = B_n (Bernoulli numbers)
  83 //   X_p(n) = \sum_{k=0}^n binomial(n,k) B_{n-k} / (k+1) * X_{p-1}(k)
  84 // The calculation of Xn depends on X0 and X{n-1}.
  85 // X_0 is special, it holds only the non-zero Bernoulli numbers with index 2 or greater.
  86 // This results in a slightly more complicated algorithm for the X_n.
  87 // The first index in Xn corresponds to the index of the polylog minus 2.
  88 // The second index in Xn corresponds to the index from the EMS.
  89 static void fill_Xn(int n)
  90 {
  91         // rule of thumb. needs to be improved. TODO
  92         const int initsize = Digits * 3 / 2;
  93
  94         if (n>1) {
  95                 // calculate X_2 and higher (corresponding to Li_4 and higher)
  96                 std::vector<cln::cl_N> buf(initsize);
  97                 std::vector<cln::cl_N>::iterator it = buf.begin();
  98                 cln::cl_N result;
  99                 *it = -(cln::expt(cln::cl_I(2),n+1) - 1) / cln::expt(cln::cl_I(2),n+1); // i == 1
 100                 it++;
 101                 for (int i=2; i<=initsize; i++) {
 102                         if (i&1) {
 103                                 result = 0; // k == 0
 104                         } else {
 105                                 result = Xn[0][i/2-1]; // k == 0
 106                         }
 107                         for (int k=1; k<i-1; k++) {
 108                                 if ( !(((i-k) & 1) && ((i-k) > 1)) ) {
 109                                         result = result + cln::binomial(i,k) * Xn[0][(i-k)/2-1] * Xn[n-1][k-1] / (k+1);
 110                                 }
 111                         }
 112                         result = result - cln::binomial(i,i-1) * Xn[n-1][i-2] / 2 / i; // k == i-1
 113                         result = result + Xn[n-1][i-1] / (i+1); // k == i
 114
 115                         *it = result;
 116                         it++;
 117                 }
 118                 Xn.push_back(buf);
 119         } else if (n==1) {
 120                 // special case to handle the X_0 correct
 121                 std::vector<cln::cl_N> buf(initsize);
 122                 std::vector<cln::cl_N>::iterator it = buf.begin();
 123                 cln::cl_N result;
 124                 *it = cln::cl_I(-3)/cln::cl_I(4); // i == 1
 125                 it++;
 126                 *it = cln::cl_I(17)/cln::cl_I(36); // i == 2
 127                 it++;
 128                 for (int i=3; i<=initsize; i++) {
 129                         if (i & 1) {
 130                                 result = -Xn[0][(i-3)/2]/2;
 131                                 *it = (cln::binomial(i,1)/cln::cl_I(2) + cln::binomial(i,i-1)/cln::cl_I(i))*result;
 132                                 it++;
 133                         } else {
 134                                 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
 135                                 for (int k=1; k<i/2; k++) {
 136                                         result = result + cln::binomial(i,k*2) * Xn[0][k-1] * Xn[0][i/2-k-1] / (k*2+1);
 137                                 }
 138                                 *it = result;
 139                                 it++;
 140                         }
 141                 }
 142                 Xn.push_back(buf);
 143         } else {
 144                 // calculate X_0
 145                 std::vector<cln::cl_N> buf(initsize/2);
 146                 std::vector<cln::cl_N>::iterator it = buf.begin();
 147                 for (int i=1; i<=initsize/2; i++) {
 148                         *it = bernoulli(i*2).to_cl_N();
 149                         it++;
 150                 }
 151                 Xn.push_back(buf);
 152         }
 153
 154         xnsize++;
 155 }
 156
 157
 158 // calculates Li(2,x) without EuMac
 159 static cln::cl_N Li2_do_sum(const cln::cl_N& x)
 160 {
 161         cln::cl_N res = x;
 162         cln::cl_N resbuf;
 163         cln::cl_N num = x;
 164         cln::cl_I den = 1; // n^2 = 1
 165         unsigned i = 3;
 166         do {
 167                 resbuf = res;
 168                 num = num * x;
 169                 den = den + i;  // n^2 = 4, 9, 16, ...
 170                 i += 2;
 171                 res = res + num / den;
 172         } while (res != resbuf);
 173         return res;
 174 }
 175
 176
 177 // calculates Li(2,x) with EuMac
 178 static cln::cl_N Li2_do_sum_EuMac(const cln::cl_N& x)
 179 {
 180         std::vector<cln::cl_N>::const_iterator it = Xn[0].begin();
 181         cln::cl_N u = -cln::log(1-x);
 182         cln::cl_N factor = u;
 183         cln::cl_N res = u - u*u/4;
 184         cln::cl_N resbuf;
 185         unsigned i = 1;
 186         do {
 187                 resbuf = res;
 188                 factor = factor * u*u / (2*i * (2*i+1));
 189                 res = res + (*it) * factor;
 190                 it++; // should we check it? or rely on initsize? ...
 191                 i++;
 192         } while (res != resbuf);
 193         return res;
 194 }
 195
 196
 197 // calculates Li(n,x), n>2 without EuMac
 198 static cln::cl_N Lin_do_sum(int n, const cln::cl_N& x)
 199 {
 200         cln::cl_N factor = x;
 201         cln::cl_N res = x;
 202         cln::cl_N resbuf;
 203         int i=2;
 204         do {
 205                 resbuf = res;
 206                 factor = factor * x;
 207                 res = res + factor / cln::expt(cln::cl_I(i),n);
 208                 i++;
 209         } while (res != resbuf);
 210         return res;
 211 }
 212
 213
 214 // calculates Li(n,x), n>2 with EuMac
 215 static cln::cl_N Lin_do_sum_EuMac(int n, const cln::cl_N& x)
 216 {
 217         std::vector<cln::cl_N>::const_iterator it = Xn[n-2].begin();
 218         cln::cl_N u = -cln::log(1-x);
 219         cln::cl_N factor = u;
 220         cln::cl_N res = u;
 221         cln::cl_N resbuf;
 222         unsigned i=2;
 223         do {
 224                 resbuf = res;
 225                 factor = factor * u / i;
 226                 res = res + (*it) * factor;
 227                 it++; // should we check it? or rely on initsize? ...
 228                 i++;
 229         } while (res != resbuf);
 230         return res;
 231 }
 232
 233
 234 // forward declaration needed by function Li_projection and C below
 235 static numeric S_num(int n, int p, const numeric& x);
 236
 237
 238 // helper function for classical polylog Li
 239 static cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
 240 {
 241         // treat n=2 as special case
 242         if (n == 2) {
 243                 // check if precalculated X0 exists
 244                 if (xnsize == 0) {
 245                         fill_Xn(0);
 246                 }
 247
 248                 if (cln::realpart(x) < 0.5) {
 249                         // choose the faster algorithm
 250                         // the switching point was empirically determined. the optimal point
 251                         // depends on hardware, Digits, ... so an approx value is okay.
 252                         // it solves also the problem with precision due to the u=-log(1-x) transformation
 253                         if (cln::abs(cln::realpart(x)) < 0.25) {
 254
 255                                 return Li2_do_sum(x);
 256                         } else {
 257                                 return Li2_do_sum_EuMac(x);
 258                         }
 259                 } else {
 260                         // choose the faster algorithm
 261                         if (cln::abs(cln::realpart(x)) > 0.75) {
 262                                 return -Li2_do_sum(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
 263                         } else {
 264                                 return -Li2_do_sum_EuMac(1-x) - cln::log(x) * cln::log(1-x) + cln::zeta(2);
 265                         }
 266                 }
 267         } else {
 268                 // check if precalculated Xn exist
 269                 if (n > xnsize+1) {
 270                         for (int i=xnsize; i<n-1; i++) {
 271                                 fill_Xn(i);
 272                         }
 273                 }
 274
 275                 if (cln::realpart(x) < 0.5) {
 276                         // choose the faster algorithm
 277                         // with n>=12 the "normal" summation always wins against EuMac
 278                         if ((cln::abs(cln::realpart(x)) < 0.3) || (n >= 12)) {
 279                                 return Lin_do_sum(n, x);
 280                         } else {
 281                                 return Lin_do_sum_EuMac(n, x);
 282                         }
 283                 } else {
 284                         cln::cl_N result = -cln::expt(cln::log(x), n-1) * cln::log(1-x) / cln::factorial(n-1);
 285                         for (int j=0; j<n-1; j++) {
 286                                 result = result + (S_num(n-j-1, 1, 1).to_cl_N() - S_num(1, n-j-1, 1-x).to_cl_N())
 287                                         * cln::expt(cln::log(x), j) / cln::factorial(j);
 288                         }
 289                         return result;
 290                 }
 291         }
 292 }
 293
 294
 295 // helper function for classical polylog Li
 296 static numeric Li_num(int n, const numeric& x)
 297 {
 298         if (n == 1) {
 299                 // just a log
 300                 return -cln::log(1-x.to_cl_N());
 301         }
 302         if (x.is_zero()) {
 303                 return 0;
 304         }
 305         if (x == 1) {
 306                 // [Kol] (2.22)
 307                 return cln::zeta(n);
 308         }
 309         else if (x == -1) {
 310                 // [Kol] (2.22)
 311                 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
 312         }
 313
 314         // what is the desired float format?
 315         // first guess: default format
 316         cln::float_format_t prec = cln::default_float_format;
 317         const cln::cl_N value = x.to_cl_N();
 318         // second guess: the argument's format
 319         if (!x.real().is_rational())
 320                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
 321         else if (!x.imag().is_rational())
 322                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 323
 324         // [Kol] (5.15)
 325         if (cln::abs(value) > 1) {
 326                 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
 327                 // check if argument is complex. if it is real, the new polylog has to be conjugated.
 328                 if (cln::zerop(cln::imagpart(value))) {
 329                         if (n & 1) {
 330                                 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
 331                         }
 332                         else {
 333                                 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
 334                         }
 335                 }
 336                 else {
 337                         if (n & 1) {
 338                                 result = result + Li_projection(n, cln::recip(value), prec);
 339                         }
 340                         else {
 341                                 result = result - Li_projection(n, cln::recip(value), prec);
 342                         }
 343                 }
 344                 cln::cl_N add;
 345                 for (int j=0; j<n-1; j++) {
 346                         add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
 347                                         * Li_num(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
 348                 }
 349                 result = result - add;
 350                 return result;
 351         }
 352         else {
 353                 return Li_projection(n, value, prec);
 354         }
 355 }
 356
 357
 358 //////////////////////////////////////////////////////////////////////
 359 //
 360 // Multiple polylogarithm  Li
 361 //
 362 // helper function
 363 //
 364 //////////////////////////////////////////////////////////////////////
 365
 366
 367 static cln::cl_N multipleLi_do_sum(const std::vector<int>& s, const std::vector<cln::cl_N>& x)
 368 {
 369         const int j = s.size();
 370
 371         std::vector<cln::cl_N> t(j);
 372         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
 373
 374         cln::cl_N t0buf;
 375         int q = 0;
 376         do {
 377                 t0buf = t[0];
 378                 q++;
 379                 t[j-1] = t[j-1] + cln::expt(x[j-1], q) / cln::expt(cln::cl_I(q),s[j-1]) * one;
 380                 for (int k=j-2; k>=0; k--) {
 381                         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]);
 382                 }
 383         } while (t[0] != t0buf);
 384
 385         return t[0];
 386 }
 387
 388
 389 //////////////////////////////////////////////////////////////////////
 390 //
 391 // Classical polylogarithm and multiple polylogarithm  Li
 392 //
 393 // GiNaC function
 394 //
 395 //////////////////////////////////////////////////////////////////////
 396
 397
 398 static ex Li_eval(const ex& x1, const ex& x2)
 399 {
 400         if (x2.is_zero()) {
 401                 return 0;
 402         }
 403         else {
 404                 if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational)))
 405                         return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
 406                 if (is_a<lst>(x2)) {
 407                         for (int i=0; i<x2.nops(); i++) {
 408                                 if (!is_a<numeric>(x2.op(i))) {
 409                                         return Li(x1,x2).hold();
 410                                 }
 411                         }
 412                         return Li(x1,x2).evalf();
 413                 }
 414                 return Li(x1,x2).hold();
 415         }
 416 }
 417
 418
 419 static ex Li_evalf(const ex& x1, const ex& x2)
 420 {
 421         // classical polylogs
 422         if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
 423                 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
 424         }
 425         // multiple polylogs
 426         else if (is_a<lst>(x1) && is_a<lst>(x2)) {
 427                 for (int i=0; i<x1.nops(); i++) {
 428                         if (!x1.op(i).info(info_flags::posint)) {
 429                                 return Li(x1,x2).hold();
 430                         }
 431                         if (!is_a<numeric>(x2.op(i))) {
 432                                 return Li(x1,x2).hold();
 433                         }
 434                         if (x2.op(i) >= 1) {
 435                                 return Li(x1,x2).hold();
 436                         }
 437                 }
 438
 439                 std::vector<int> m;
 440                 std::vector<cln::cl_N> x;
 441                 for (int i=ex_to<numeric>(x1.nops()).to_int()-1; i>=0; i--) {
 442                         m.push_back(ex_to<numeric>(x1.op(i)).to_int());
 443                         x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
 444                 }
 445
 446                 return numeric(multipleLi_do_sum(m, x));
 447         }
 448
 449         return Li(x1,x2).hold();
 450 }
 451
 452
 453 static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
 454 {
 455         epvector seq;
 456         seq.push_back(expair(Li(x1,x2), 0));
 457         return pseries(rel,seq);
 458 }
 459
 460
 461 static ex Li_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
 462 {
 463         GINAC_ASSERT(deriv_param < 2);
 464         if (deriv_param == 0) {
 465                 return 0;
 466         }
 467         if (x1 > 0) {
 468                 return Li(x1-1, x2) / x2;
 469         } else {
 470                 return 1/(1-x2);
 471         }
 472 }
 473
 474
 475 REGISTER_FUNCTION(Li,
 476                 eval_func(Li_eval).
 477                 evalf_func(Li_evalf).
 478                 do_not_evalf_params().
 479                 series_func(Li_series).
 480                 derivative_func(Li_deriv));
 481
 482
 483 //////////////////////////////////////////////////////////////////////
 484 //
 485 // Nielsen's generalized polylogarithm  S
 486 //
 487 // helper functions
 488 //
 489 //////////////////////////////////////////////////////////////////////
 490
 491
 492 namespace {
 493 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
 494 // see fill_Yn()
 495 std::vector<std::vector<cln::cl_N> > Yn;
 496 int ynsize = 0; // number of Yn[]
 497 int ynlength = 100; // initial length of all Yn[i]
 498 }
 499
 500
 501 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
 502 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
 503 // representing S_{n,p}(x).
 504 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
 505 // equivalent Z-sum.
 506 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
 507 // representing S_{n,p}(x).
 508 // The calculation of Y_n uses the values from Y_{n-1}.
 509 static void fill_Yn(int n, const cln::float_format_t& prec)
 510 {
 511         const int initsize = ynlength;
 512         //const int initsize = initsize_Yn;
 513         cln::cl_N one = cln::cl_float(1, prec);
 514
 515         if (n) {
 516                 std::vector<cln::cl_N> buf(initsize);
 517                 std::vector<cln::cl_N>::iterator it = buf.begin();
 518                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
 519                 *it = (*itprev) / cln::cl_N(n+1) * one;
 520                 it++;
 521                 itprev++;
 522                 // sums with an index smaller than the depth are zero and need not to be calculated.
 523                 // calculation starts with depth, which is n+2)
 524                 for (int i=n+2; i<=initsize+n; i++) {
 525                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
 526                         it++;
 527                         itprev++;
 528                 }
 529                 Yn.push_back(buf);
 530         } else {
 531                 std::vector<cln::cl_N> buf(initsize);
 532                 std::vector<cln::cl_N>::iterator it = buf.begin();
 533                 *it = 1 * one;
 534                 it++;
 535                 for (int i=2; i<=initsize; i++) {
 536                         *it = *(it-1) + 1 / cln::cl_N(i) * one;
 537                         it++;
 538                 }
 539                 Yn.push_back(buf);
 540         }
 541         ynsize++;
 542 }
 543
 544
 545 // make Yn longer ...
 546 static void make_Yn_longer(int newsize, const cln::float_format_t& prec)
 547 {
 548
 549         cln::cl_N one = cln::cl_float(1, prec);
 550
 551         Yn[0].resize(newsize);
 552         std::vector<cln::cl_N>::iterator it = Yn[0].begin();
 553         it += ynlength;
 554         for (int i=ynlength+1; i<=newsize; i++) {
 555                 *it = *(it-1) + 1 / cln::cl_N(i) * one;
 556                 it++;
 557         }
 558
 559         for (int n=1; n<ynsize; n++) {
 560                 Yn[n].resize(newsize);
 561                 std::vector<cln::cl_N>::iterator it = Yn[n].begin();
 562                 std::vector<cln::cl_N>::iterator itprev = Yn[n-1].begin();
 563                 it += ynlength;
 564                 itprev += ynlength;
 565                 for (int i=ynlength+n+1; i<=newsize+n; i++) {
 566                         *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
 567                         it++;
 568                         itprev++;
 569                 }
 570         }
 571
 572         ynlength = newsize;
 573 }
 574
 575
 576 // helper function for S(n,p,x)
 577 // [Kol] (7.2)
 578 static cln::cl_N C(int n, int p)
 579 {
 580         cln::cl_N result;
 581
 582         for (int k=0; k<p; k++) {
 583                 for (int j=0; j<=(n+k-1)/2; j++) {
 584                         if (k == 0) {
 585                                 if (n & 1) {
 586                                         if (j & 1) {
 587                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 588                                         }
 589                                         else {
 590                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 591                                         }
 592                                 }
 593                         }
 594                         else {
 595                                 if (k & 1) {
 596                                         if (j & 1) {
 597                                                 result = result + cln::factorial(n+k-1)
 598                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 599                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 600                                         }
 601                                         else {
 602                                                 result = result - cln::factorial(n+k-1)
 603                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 604                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 605                                         }
 606                                 }
 607                                 else {
 608                                         if (j & 1) {
 609                                                 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()
 610                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 611                                         }
 612                                         else {
 613                                                 result = result + cln::factorial(n+k-1)
 614                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 615                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 616                                         }
 617                                 }
 618                         }
 619                 }
 620         }
 621         int np = n+p;
 622         if ((np-1) & 1) {
 623                 if (((np)/2+n) & 1) {
 624                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 625                 }
 626                 else {
 627                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 628                 }
 629         }
 630
 631         return result;
 632 }
 633
 634
 635 // helper function for S(n,p,x)
 636 // [Kol] remark to (9.1)
 637 static cln::cl_N a_k(int k)
 638 {
 639         cln::cl_N result;
 640
 641         if (k == 0) {
 642                 return 1;
 643         }
 644
 645         result = result;
 646         for (int m=2; m<=k; m++) {
 647                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
 648         }
 649
 650         return -result / k;
 651 }
 652
 653
 654 // helper function for S(n,p,x)
 655 // [Kol] remark to (9.1)
 656 static cln::cl_N b_k(int k)
 657 {
 658         cln::cl_N result;
 659
 660         if (k == 0) {
 661                 return 1;
 662         }
 663
 664         result = result;
 665         for (int m=2; m<=k; m++) {
 666                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
 667         }
 668
 669         return result / k;
 670 }
 671
 672
 673 // helper function for S(n,p,x)
 674 static cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 675 {
 676         if (p==1) {
 677                 return Li_projection(n+1, x, prec);
 678         }
 679
 680         // check if precalculated values are sufficient
 681         if (p > ynsize+1) {
 682                 for (int i=ynsize; i<p-1; i++) {
 683                         fill_Yn(i, prec);
 684                 }
 685         }
 686
 687         // should be done otherwise
 688         cln::cl_N xf = x * cln::cl_float(1, prec);
 689
 690         cln::cl_N res;
 691         cln::cl_N resbuf;
 692         cln::cl_N factor = cln::expt(xf, p);
 693         int i = p;
 694         do {
 695                 resbuf = res;
 696                 if (i-p >= ynlength) {
 697                         // make Yn longer
 698                         make_Yn_longer(ynlength*2, prec);
 699                 }
 700                 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? ...
 701                 //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
 702                 factor = factor * xf;
 703                 i++;
 704         } while (res != resbuf);
 705
 706         return res;
 707 }
 708
 709
 710 // helper function for S(n,p,x)
 711 static cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 712 {
 713         // [Kol] (5.3)
 714         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
 715
 716                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
 717                         * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
 718
 719                 for (int s=0; s<n; s++) {
 720                         cln::cl_N res2;
 721                         for (int r=0; r<p; r++) {
 722                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
 723                                         * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
 724                         }
 725                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
 726                 }
 727
 728                 return result;
 729         }
 730
 731         return S_do_sum(n, p, x, prec);
 732 }
 733
 734
 735 // helper function for S(n,p,x)
 736 static numeric S_num(int n, int p, const numeric& x)
 737 {
 738         if (x == 1) {
 739                 if (n == 1) {
 740                     // [Kol] (2.22) with (2.21)
 741                         return cln::zeta(p+1);
 742                 }
 743
 744                 if (p == 1) {
 745                     // [Kol] (2.22)
 746                         return cln::zeta(n+1);
 747                 }
 748
 749                 // [Kol] (9.1)
 750                 cln::cl_N result;
 751                 for (int nu=0; nu<n; nu++) {
 752                         for (int rho=0; rho<=p; rho++) {
 753                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
 754                                         * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
 755                         }
 756                 }
 757                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
 758
 759                 return result;
 760         }
 761         else if (x == -1) {
 762                 // [Kol] (2.22)
 763                 if (p == 1) {
 764                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
 765                 }
 766 //              throw std::runtime_error("don't know how to evaluate this function!");
 767         }
 768
 769         // what is the desired float format?
 770         // first guess: default format
 771         cln::float_format_t prec = cln::default_float_format;
 772         const cln::cl_N value = x.to_cl_N();
 773         // second guess: the argument's format
 774         if (!x.real().is_rational())
 775                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
 776         else if (!x.imag().is_rational())
 777                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 778
 779
 780         // [Kol] (5.3)
 781         if (cln::realpart(value) < -0.5) {
 782
 783                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
 784                         * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
 785
 786                 for (int s=0; s<n; s++) {
 787                         cln::cl_N res2;
 788                         for (int r=0; r<p; r++) {
 789                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
 790                                         * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
 791                         }
 792                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
 793                 }
 794
 795                 return result;
 796
 797         }
 798         // [Kol] (5.12)
 799         if (cln::abs(value) > 1) {
 800
 801                 cln::cl_N result;
 802
 803                 for (int s=0; s<p; s++) {
 804                         for (int r=0; r<=s; r++) {
 805                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
 806                                         / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
 807                                         * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
 808                         }
 809                 }
 810                 result = result * cln::expt(cln::cl_I(-1),n);
 811
 812                 cln::cl_N res2;
 813                 for (int r=0; r<n; r++) {
 814                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
 815                 }
 816                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
 817
 818                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
 819
 820                 return result;
 821         }
 822         else {
 823                 return S_projection(n, p, value, prec);
 824         }
 825 }
 826
 827
 828 //////////////////////////////////////////////////////////////////////
 829 //
 830 // Nielsen's generalized polylogarithm  S
 831 //
 832 // GiNaC function
 833 //
 834 //////////////////////////////////////////////////////////////////////
 835
 836
 837 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
 838 {
 839         if (x2 == 1) {
 840                 return Li(x1+1,x3);
 841         }
 842         if (x3.info(info_flags::numeric) && (!x3.info(info_flags::crational)) &&
 843                         x1.info(info_flags::posint) && x2.info(info_flags::posint)) {
 844                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
 845         }
 846         return S(x1,x2,x3).hold();
 847 }
 848
 849
 850 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
 851 {
 852         if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
 853                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
 854         }
 855         return S(x1,x2,x3).hold();
 856 }
 857
 858
 859 static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options)
 860 {
 861         epvector seq;
 862         seq.push_back(expair(S(x1,x2,x3), 0));
 863         return pseries(rel,seq);
 864 }
 865
 866
 867 static ex S_deriv(const ex& x1, const ex& x2, const ex& x3, unsigned deriv_param)
 868 {
 869         GINAC_ASSERT(deriv_param < 3);
 870         if (deriv_param < 2) {
 871                 return 0;
 872         }
 873         if (x1 > 0) {
 874                 return S(x1-1, x2, x3) / x3;
 875         } else {
 876                 return S(x1, x2-1, x3) / (1-x3);
 877         }
 878 }
 879
 880
 881 REGISTER_FUNCTION(S,
 882                 eval_func(S_eval).
 883                 evalf_func(S_evalf).
 884                 do_not_evalf_params().
 885                 series_func(S_series).
 886                 derivative_func(S_deriv));
 887
 888
 889 //////////////////////////////////////////////////////////////////////
 890 //
 891 // Harmonic polylogarithm  H
 892 //
 893 // helper function
 894 //
 895 //////////////////////////////////////////////////////////////////////
 896
 897
 898 static cln::cl_N H_do_sum(const std::vector<int>& s, const cln::cl_N& x)
 899 {
 900         const int j = s.size();
 901
 902         std::vector<cln::cl_N> t(j);
 903
 904         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
 905         cln::cl_N factor = cln::expt(x, j) * one;
 906         cln::cl_N t0buf;
 907         int q = 0;
 908         do {
 909                 t0buf = t[0];
 910                 q++;
 911                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
 912                 for (int k=j-2; k>=1; k--) {
 913                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
 914                 }
 915                 t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), s[0]);
 916                 factor = factor * x;
 917         } while (t[0] != t0buf);
 918
 919         return t[0];
 920 }
 921
 922
 923 //////////////////////////////////////////////////////////////////////
 924 //
 925 // Harmonic polylogarithm  H
 926 //
 927 // GiNaC function
 928 //
 929 //////////////////////////////////////////////////////////////////////
 930
 931
 932 static ex H_eval(const ex& x1, const ex& x2)
 933 {
 934         if (x2.info(info_flags::numeric) && (!x2.info(info_flags::crational))) {
 935                 return H(x1,x2).evalf();
 936         }
 937         return H(x1,x2).hold();
 938 }
 939
 940
 941 static ex H_evalf(const ex& x1, const ex& x2)
 942 {
 943         if (is_a<lst>(x1) && is_a<numeric>(x2)) {
 944                 for (int i=0; i<x1.nops(); i++) {
 945                         if (!x1.op(i).info(info_flags::posint)) {
 946                                 return H(x1,x2).hold();
 947                         }
 948                 }
 949                 if (x1.nops() < 1) {
 950                         return 1;
 951                 }
 952                 cln::cl_N x = ex_to<numeric>(x2).to_cl_N();
 953                 if (x == 1) {
 954                         lst x1rev;
 955                         for (int i=x1.nops()-1; i>=0; i--) {
 956                                 x1rev.append(x1.op(i));
 957                         }
 958                         return mZeta(x1rev).evalf();
 959                 }
 960                 const int j = x1.nops();
 961                 if (x2 > 1 || j < 2) {
 962                         return H(x1,x2).hold();
 963                 }
 964
 965                 std::vector<int> r(j);
 966                 for (int i=0; i<j; i++) {
 967                         r[i] = ex_to<numeric>(x1.op(i)).to_int();
 968                 }
 969
 970                 return numeric(H_do_sum(r,x));
 971         }
 972
 973         return H(x1,x2).hold();
 974 }
 975
 976
 977 static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
 978 {
 979         epvector seq;
 980         seq.push_back(expair(H(x1,x2), 0));
 981         return pseries(rel,seq);
 982 }
 983
 984
 985 static ex H_deriv(const ex& x1, const ex& x2, unsigned deriv_param)
 986 {
 987         GINAC_ASSERT(deriv_param < 2);
 988         if (deriv_param == 0) {
 989                 return 0;
 990         }
 991         if (is_a<lst>(x1)) {
 992                 lst newparameter = ex_to<lst>(x1);
 993                 if (x1.op(0) == 1) {
 994                         newparameter.remove_first();
 995                         return 1/(1-x2) * H(newparameter, x2);
 996                 } else {
 997                         newparameter[0]--;
 998                         return H(newparameter, x2).hold() / x2;
 999                 }
1000         } else {
1001                 if (x1 == 1) {
1002                         return 1/(1-x2);
1003                 } else {
1004                         return H(x1-1, x2).hold() / x2;
1005                 }
1006         }
1007 }
1008
1009
1010 REGISTER_FUNCTION(H,
1011                 eval_func(H_eval).
1012                 evalf_func(H_evalf).
1013                 do_not_evalf_params().
1014                 series_func(H_series).
1015                 derivative_func(H_deriv));
1016
1017
1018 //////////////////////////////////////////////////////////////////////
1019 //
1020 // Multiple zeta values  mZeta
1021 //
1022 // helper functions
1023 //
1024 //////////////////////////////////////////////////////////////////////
1025
1026
1027 namespace {
1028 const cln::cl_N lambda = cln::cl_N("319/320");
1029 int L1;
1030 int L2;
1031 std::vector<std::vector<cln::cl_N> > f_kj;
1032 std::vector<cln::cl_N> crB;
1033 std::vector<std::vector<cln::cl_N> > crG;
1034 std::vector<cln::cl_N> crX;
1035 }
1036
1037
1038 static void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
1039 {
1040         const int size = a.size();
1041         for (int n=0; n<size; n++) {
1042                 c[n] = 0;
1043                 for (int m=0; m<=n; m++) {
1044                         c[n] = c[n] + a[m]*b[n-m];
1045                 }
1046         }
1047 }
1048
1049
1050 // [Cra] section 4
1051 static void initcX(const std::vector<int>& s)
1052 {
1053         const int k = s.size();
1054
1055         crX.clear();
1056         crG.clear();
1057         crB.clear();
1058
1059         for (int i=0; i<=L2; i++) {
1060                 crB.push_back(bernoulli(i).to_cl_N() / cln::factorial(i));
1061         }
1062
1063         int Sm = 0;
1064         int Smp1 = 0;
1065         for (int m=0; m<k-1; m++) {
1066                 std::vector<cln::cl_N> crGbuf;
1067                 Sm = Sm + s[m];
1068                 Smp1 = Sm + s[m+1];
1069                 for (int i=0; i<=L2; i++) {
1070                         crGbuf.push_back(cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2));
1071                 }
1072                 crG.push_back(crGbuf);
1073         }
1074
1075         crX = crB;
1076
1077         for (int m=0; m<k-1; m++) {
1078                 std::vector<cln::cl_N> Xbuf;
1079                 for (int i=0; i<=L2; i++) {
1080                         Xbuf.push_back(crX[i] * crG[m][i]);
1081                 }
1082                 halfcyclic_convolute(Xbuf, crB, crX);
1083         }
1084 }
1085
1086
1087 // [Cra] section 4
1088 static cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk)
1089 {
1090         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1091         cln::cl_N factor = cln::expt(lambda, Sqk);
1092         cln::cl_N res = factor / Sqk * crX[0] * one;
1093         cln::cl_N resbuf;
1094         int N = 0;
1095         do {
1096                 resbuf = res;
1097                 factor = factor * lambda;
1098                 N++;
1099                 res = res + crX[N] * factor / (N+Sqk);
1100         } while ((res != resbuf) || cln::zerop(crX[N]));
1101         return res;
1102 }
1103
1104
1105 // [Cra] section 4
1106 static void calc_f(int maxr)
1107 {
1108         f_kj.clear();
1109         f_kj.resize(L1);
1110
1111         cln::cl_N t0, t1, t2, t3, t4;
1112         int i, j, k;
1113         std::vector<std::vector<cln::cl_N> >::iterator it = f_kj.begin();
1114         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1115
1116         t0 = cln::exp(-lambda);
1117         t2 = 1;
1118         for (k=1; k<=L1; k++) {
1119                 t1 = k * lambda;
1120                 t2 = t0 * t2;
1121                 for (j=1; j<=maxr; j++) {
1122                         t3 = 1;
1123                         t4 = 1;
1124                         for (i=2; i<=j; i++) {
1125                                 t4 = t4 * (j-i+1);
1126                                 t3 = t1 * t3 + t4;
1127                         }
1128                         (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one);
1129                 }
1130                 it++;
1131         }
1132 }
1133
1134
1135 // [Cra] (3.1)
1136 static cln::cl_N crandall_Z(const std::vector<int>& s)
1137 {
1138         const int j = s.size();
1139
1140         if (j == 1) {
1141                 cln::cl_N t0;
1142                 cln::cl_N t0buf;
1143                 int q = 0;
1144                 do {
1145                         t0buf = t0;
1146                         q++;
1147                         t0 = t0 + f_kj[q+j-2][s[0]-1];
1148                 } while (t0 != t0buf);
1149
1150                 return t0 / cln::factorial(s[0]-1);
1151         }
1152
1153         std::vector<cln::cl_N> t(j);
1154
1155         cln::cl_N t0buf;
1156         int q = 0;
1157         do {
1158                 t0buf = t[0];
1159                 q++;
1160                 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
1161                 for (int k=j-2; k>=1; k--) {
1162                         t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
1163                 }
1164                 t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
1165         } while (t[0] != t0buf);
1166
1167         return t[0] / cln::factorial(s[0]-1);
1168 }
1169
1170
1171 // [Cra] (2.4)
1172 static cln::cl_N mZeta_do_sum_Crandall(const std::vector<int>& s)
1173 {
1174         std::vector<int> r = s;
1175         const int j = r.size();
1176
1177         // decide on maximal size of f_kj for crandall_Z
1178         if (Digits < 50) {
1179                 L1 = 150;
1180         } else {
1181                 L1 = Digits * 3 + j*2;
1182         }
1183
1184         // decide on maximal size of crX for crandall_Y
1185         if (Digits < 38) {
1186                 L2 = 63;
1187         } else if (Digits < 86) {
1188                 L2 = 127;
1189         } else if (Digits < 192) {
1190                 L2 = 255;
1191         } else if (Digits < 394) {
1192                 L2 = 511;
1193         } else if (Digits < 808) {
1194                 L2 = 1023;
1195         } else {
1196                 L2 = 2047;
1197         }
1198
1199         cln::cl_N res;
1200
1201         int maxr = 0;
1202         int S = 0;
1203         for (int i=0; i<j; i++) {
1204                 S += r[i];
1205                 if (r[i] > maxr) {
1206                         maxr = r[i];
1207                 }
1208         }
1209
1210         calc_f(maxr);
1211
1212         const cln::cl_N r0factorial = cln::factorial(r[0]-1);
1213
1214         std::vector<int> rz;
1215         int skp1buf;
1216         int Srun = S;
1217         for (int k=r.size()-1; k>0; k--) {
1218
1219                 rz.insert(rz.begin(), r.back());
1220                 skp1buf = rz.front();
1221                 Srun -= skp1buf;
1222                 r.pop_back();
1223
1224                 initcX(r);
1225
1226                 for (int q=0; q<skp1buf; q++) {
1227
1228                         cln::cl_N pp1 = crandall_Y_loop(Srun+q-k);
1229                         cln::cl_N pp2 = crandall_Z(rz);
1230
1231                         rz.front()--;
1232
1233                         if (q & 1) {
1234                                 res = res - pp1 * pp2 / cln::factorial(q);
1235                         } else {
1236                                 res = res + pp1 * pp2 / cln::factorial(q);
1237                         }
1238                 }
1239                 rz.front() = skp1buf;
1240         }
1241         rz.insert(rz.begin(), r.back());
1242
1243         initcX(rz);
1244
1245         res = (res + crandall_Y_loop(S-j)) / r0factorial + crandall_Z(rz);
1246
1247         return res;
1248 }
1249
1250
1251 static cln::cl_N mZeta_do_sum_simple(const std::vector<int>& r)
1252 {
1253         const int j = r.size();
1254
1255         // buffer for subsums
1256         std::vector<cln::cl_N> t(j);
1257         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
1258
1259         cln::cl_N t0buf;
1260         int q = 0;
1261         do {
1262                 t0buf = t[0];
1263                 q++;
1264                 t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
1265                 for (int k=j-2; k>=0; k--) {
1266                         t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
1267                 }
1268         } while (t[0] != t0buf);
1269
1270         return t[0];
1271 }
1272
1273
1274 //////////////////////////////////////////////////////////////////////
1275 //
1276 // Multiple zeta values  mZeta
1277 //
1278 // GiNaC function
1279 //
1280 //////////////////////////////////////////////////////////////////////
1281
1282
1283 static ex mZeta_eval(const ex& x1)
1284 {
1285         return mZeta(x1).hold();
1286 }
1287
1288
1289 static ex mZeta_evalf(const ex& x1)
1290 {
1291         if (is_a<lst>(x1)) {
1292                 for (int i=0; i<x1.nops(); i++) {
1293                         if (!x1.op(i).info(info_flags::posint))
1294                                 return mZeta(x1).hold();
1295                 }
1296
1297                 const int j = x1.nops();
1298
1299                 std::vector<int> r(j);
1300                 for (int i=0; i<j; i++) {
1301                         r[j-1-i] = ex_to<numeric>(x1.op(i)).to_int();
1302                 }
1303
1304                 // check for divergence
1305                 if (r[0] == 1) {
1306                         return mZeta(x1).hold();
1307                 }
1308
1309                 // if only one argument, use cln::zeta
1310                 if (j == 1) {
1311                         return numeric(cln::zeta(r[0]));
1312                 }
1313
1314                 // decide on summation algorithm
1315                 // this is still a bit clumsy
1316                 int limit = (Digits>17) ? 10 : 6;
1317                 if (r[0]<limit || (j>3 && r[1]<limit/2)) {
1318                         return numeric(mZeta_do_sum_Crandall(r));
1319                 } else {
1320                         return numeric(mZeta_do_sum_simple(r));
1321                 }
1322         } else if (x1.info(info_flags::posint) && (x1 != 1)) {
1323                 return numeric(cln::zeta(ex_to<numeric>(x1).to_int()));
1324         }
1325
1326         return mZeta(x1).hold();
1327 }
1328
1329
1330 static ex mZeta_series(const ex& x1, const relational& rel, int order, unsigned options)
1331 {
1332         epvector seq;
1333         seq.push_back(expair(mZeta(x1), 0));
1334         return pseries(rel,seq);
1335 }
1336
1337
1338 static ex mZeta_deriv(const ex& x, unsigned deriv_param)
1339 {
1340         return 0;
1341 }
1342
1343
1344 REGISTER_FUNCTION(mZeta,
1345                 eval_func(mZeta_eval).
1346                 evalf_func(mZeta_evalf).
1347                 do_not_evalf_params().series_func(mZeta_series).
1348                 derivative_func(mZeta_deriv));
1349
1350
1351 } // namespace GiNaC
1352