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 publication:
  13  *      Nielsen's Generalized Polylogarithms, K.S.Kolbig, SIAM J.Math.Anal. 17 (1986), pp. 1232-1258.
  14  *      This document will be referenced as [Kol] throughout this source code.
  15  *    - Classical polylogarithms (Li) and nielsen's generalized polylogarithms (S) can be numerically
  16  *      evaluated in the whole complex plane except for S(n,p,-1) when p is not unit (no formula yet
  17  *      to tackle these points). And of course, there is still room for speed optimizations ;-).
  18  *    - The remaining functions can only be numerically evaluated with arguments lying in the unit sphere
  19  *      at the moment. Sorry. The evaluation especially for mZeta is very slow ... better not use it
  20  *      right now.
  21  *    - The functions have no series expansion. To do it, you have to convert these functions
  22  *      into the appropriate objects from the nestedsums library, do the expansion and convert the
  23  *      result back.
  24  *    - Numerical testing of this implementation has been performed by doing a comparison of results
  25  *      between this software and the commercial M.......... 4.1.
  26  *
  27  */
  28
  29 /*
  30  *  GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
  31  *
  32  *  This program is free software; you can redistribute it and/or modify
  33  *  it under the terms of the GNU General Public License as published by
  34  *  the Free Software Foundation; either version 2 of the License, or
  35  *  (at your option) any later version.
  36  *
  37  *  This program is distributed in the hope that it will be useful,
  38  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  39  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  40  *  GNU General Public License for more details.
  41  *
  42  *  You should have received a copy of the GNU General Public License
  43  *  along with this program; if not, write to the Free Software
  44  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  45  */
  46
  47 #include <stdexcept>
  48 #include <vector>
  49 #include <cln/cln.h>
  50
  51 #include "inifcns.h"
  52 #include "lst.h"
  53 #include "numeric.h"
  54 #include "operators.h"
  55 #include "relational.h"
  56 #include "pseries.h"
  57
  58
  59 namespace GiNaC {
  60
  61
  62 //////////////////////
  63 // helper functions //
  64 //////////////////////
  65
  66
  67 // helper function for classical polylog Li
  68 static cln::cl_N Li_series(int n, const cln::cl_N& x, const cln::float_format_t& prec)
  69 {
  70         // Note: argument must be in the unit circle
  71         cln::cl_N aug, acc;
  72         cln::cl_N num = cln::complex(cln::cl_float(1, prec), 0);
  73         cln::cl_N den = 0;
  74         int i = 1;
  75         do {
  76                 num = num * x;
  77                 cln::cl_R ii = i;
  78                 den = cln::expt(ii, n);
  79                 i++;
  80                 aug = num / den;
  81                 acc = acc + aug;
  82         } while (acc != acc+aug);
  83         return acc;
  84 }
  85
  86
  87 // helper function for classical polylog Li
  88 static cln::cl_N Li_projection(int n, const cln::cl_N& x, const cln::float_format_t& prec)
  89 {
  90         return Li_series(n, x, prec);
  91 }
  92
  93
  94 // helper function for classical polylog Li
  95 static numeric Li_num(int n, const numeric& x)
  96 {
  97         if (n == 1) {
  98                 // just a log
  99                 return -cln::log(1-x.to_cl_N());
 100         }
 101         if (x.is_zero()) {
 102                 return 0;
 103         }
 104         if (x == 1) {
 105                 // [Kol] (2.22)
 106                 return cln::zeta(n);
 107         }
 108         else if (x == -1) {
 109                 // [Kol] (2.22)
 110                 return -(1-cln::expt(cln::cl_I(2),1-n)) * cln::zeta(n);
 111         }
 112
 113         // what is the desired float format?
 114         // first guess: default format
 115         cln::float_format_t prec = cln::default_float_format;
 116         const cln::cl_N value = x.to_cl_N();
 117         // second guess: the argument's format
 118         if (!x.real().is_rational())
 119                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
 120         else if (!x.imag().is_rational())
 121                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 122
 123         // [Kol] (5.15)
 124         if (cln::abs(value) > 1) {
 125                 cln::cl_N result = -cln::expt(cln::log(-value),n) / cln::factorial(n);
 126                 // check if argument is complex. if it is real, the new polylog has to be conjugated.
 127                 if (cln::zerop(cln::imagpart(value))) {
 128                         if (n & 1) {
 129                                 result = result + conjugate(Li_projection(n, cln::recip(value), prec));
 130                         }
 131                         else {
 132                                 result = result - conjugate(Li_projection(n, cln::recip(value), prec));
 133                         }
 134                 }
 135                 else {
 136                         if (n & 1) {
 137                                 result = result + Li_projection(n, cln::recip(value), prec);
 138                         }
 139                         else {
 140                                 result = result - Li_projection(n, cln::recip(value), prec);
 141                         }
 142                 }
 143                 cln::cl_N add;
 144                 for (int j=0; j<n-1; j++) {
 145                         add = add + (1+cln::expt(cln::cl_I(-1),n-j)) * (1-cln::expt(cln::cl_I(2),1-n+j))
 146                                         * Li_num(n-j,1).to_cl_N() * cln::expt(cln::log(-value),j) / cln::factorial(j);
 147                 }
 148                 result = result - add;
 149                 return result;
 150         }
 151         else {
 152                 return Li_projection(n, value, prec);
 153         }
 154 }
 155
 156
 157 // helper function for S(n,p,x)
 158 static cln::cl_N numeric_nielsen(int n, int step)
 159 {
 160         if (step) {
 161                 cln::cl_N res;
 162                 for (int i=1; i<n; i++) {
 163                         res = res + numeric_nielsen(i, step-1) / cln::cl_I(i);
 164                 }
 165                 return res;
 166         }
 167         else {
 168                 return 1;
 169         }
 170 }
 171
 172
 173 // forward declaration needed by function C below
 174 static numeric S_num(int n, int p, const numeric& x);
 175
 176
 177 // helper function for S(n,p,x)
 178 // [Kol] (7.2)
 179 static cln::cl_N C(int n, int p)
 180 {
 181         cln::cl_N result;
 182
 183         for (int k=0; k<p; k++) {
 184                 for (int j=0; j<=(n+k-1)/2; j++) {
 185                         if (k == 0) {
 186                                 if (n & 1) {
 187                                         if (j & 1) {
 188                                                 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 189                                         }
 190                                         else {
 191                                                 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1).to_cl_N() / cln::factorial(2*j);
 192                                         }
 193                                 }
 194                         }
 195                         else {
 196                                 if (k & 1) {
 197                                         if (j & 1) {
 198                                                 result = result + cln::factorial(n+k-1)
 199                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 200                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 201                                         }
 202                                         else {
 203                                                 result = result - cln::factorial(n+k-1)
 204                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 205                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 206                                         }
 207                                 }
 208                                 else {
 209                                         if (j & 1) {
 210                                                 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()
 211                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 212                                         }
 213                                         else {
 214                                                 result = result + cln::factorial(n+k-1)
 215                                                         * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1).to_cl_N()
 216                                                         / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
 217                                         }
 218                                 }
 219                         }
 220                 }
 221         }
 222         int np = n+p;
 223         if ((np-1) & 1) {
 224                 if (((np)/2+n) & 1) {
 225                         result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 226                 }
 227                 else {
 228                         result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
 229                 }
 230         }
 231
 232         return result;
 233 }
 234
 235
 236 // helper function for S(n,p,x)
 237 // [Kol] remark to (9.1)
 238 static cln::cl_N a_k(int k)
 239 {
 240         cln::cl_N result;
 241
 242         if (k == 0) {
 243                 return 1;
 244         }
 245
 246         result = result;
 247         for (int m=2; m<=k; m++) {
 248                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
 249         }
 250
 251         return -result / k;
 252 }
 253
 254
 255 // helper function for S(n,p,x)
 256 // [Kol] remark to (9.1)
 257 static cln::cl_N b_k(int k)
 258 {
 259         cln::cl_N result;
 260
 261         if (k == 0) {
 262                 return 1;
 263         }
 264
 265         result = result;
 266         for (int m=2; m<=k; m++) {
 267                 result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
 268         }
 269
 270         return result / k;
 271 }
 272
 273
 274 // helper function for S(n,p,x)
 275 static cln::cl_N S_series(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 276 {
 277         n++;
 278         int i = p;
 279         p--;
 280         cln::cl_N aug, acc;
 281         cln::cl_N num = cln::expt(x,p);
 282         cln::cl_N converter = cln::complex(cln::cl_float(1, prec), 0);
 283         cln::cl_N den = 0;
 284         do {
 285                 num = num * x;
 286                 den = cln::expt(cln::cl_I(i), n);
 287                 aug = num / den * numeric_nielsen(i, p);
 288                 i++;
 289                 acc = acc + aug;
 290         } while (acc != acc+aug);
 291
 292         return acc;
 293 }
 294
 295
 296 // helper function for S(n,p,x)
 297 static cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
 298 {
 299         // [Kol] (5.3)
 300         if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
 301
 302                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
 303                         * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
 304
 305                 for (int s=0; s<n; s++) {
 306                         cln::cl_N res2;
 307                         for (int r=0; r<p; r++) {
 308                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
 309                                         * S_series(p-r,n-s,1-x,prec) / cln::factorial(r);
 310                         }
 311                         result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
 312                 }
 313
 314                 return result;
 315         }
 316
 317         return S_series(n, p, x, prec);
 318 }
 319
 320
 321 // helper function for S(n,p,x)
 322 static numeric S_num(int n, int p, const numeric& x)
 323 {
 324         if (x == 1) {
 325                 if (n == 1) {
 326                     // [Kol] (2.22) with (2.21)
 327                         return cln::zeta(p+1);
 328                 }
 329
 330                 if (p == 1) {
 331                     // [Kol] (2.22)
 332                         return cln::zeta(n+1);
 333                 }
 334
 335                 // [Kol] (9.1)
 336                 cln::cl_N result;
 337                 for (int nu=0; nu<n; nu++) {
 338                         for (int rho=0; rho<=p; rho++) {
 339                                 result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
 340                                         * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
 341                         }
 342                 }
 343                 result = result * cln::expt(cln::cl_I(-1),n+p-1);
 344
 345                 return result;
 346         }
 347         else if (x == -1) {
 348                 // [Kol] (2.22)
 349                 if (p == 1) {
 350                         return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
 351                 }
 352                 throw std::runtime_error("don't know how to evaluate this function!");
 353         }
 354
 355         // what is the desired float format?
 356         // first guess: default format
 357         cln::float_format_t prec = cln::default_float_format;
 358         const cln::cl_N value = x.to_cl_N();
 359         // second guess: the argument's format
 360         if (!x.real().is_rational())
 361                 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
 362         else if (!x.imag().is_rational())
 363                 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
 364
 365
 366         // [Kol] (5.3)
 367         if (cln::realpart(value) < -0.5) {
 368
 369                 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
 370                         * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
 371
 372                 for (int s=0; s<n; s++) {
 373                         cln::cl_N res2;
 374                         for (int r=0; r<p; r++) {
 375                                 res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
 376                                         * S_num(p-r,n-s,1-value).to_cl_N() / cln::factorial(r);
 377                         }
 378                         result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1).to_cl_N() - res2) / cln::factorial(s);
 379                 }
 380
 381                 return result;
 382
 383         }
 384         // [Kol] (5.12)
 385         else if (cln::abs(value) > 1) {
 386
 387                 cln::cl_N result;
 388
 389                 for (int s=0; s<p; s++) {
 390                         for (int r=0; r<=s; r++) {
 391                                 result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
 392                                         / cln::factorial(r) / cln::factorial(s-r) / cln::factorial(n-1)
 393                                         * S_num(n+s-r,p-s,cln::recip(value)).to_cl_N();
 394                         }
 395                 }
 396                 result = result * cln::expt(cln::cl_I(-1),n);
 397
 398                 cln::cl_N res2;
 399                 for (int r=0; r<n; r++) {
 400                         res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
 401                 }
 402                 res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
 403
 404                 result = result + cln::expt(cln::cl_I(-1),p) * res2;
 405
 406                 return result;
 407         }
 408         else {
 409                 return S_projection(n, p, value, prec);
 410         }
 411 }
 412
 413
 414 // helper function for multiple polylogarithm
 415 static cln::cl_N numeric_zsum(int n, std::vector<cln::cl_N>& x, std::vector<cln::cl_N>& m)
 416 {
 417         cln::cl_N res;
 418         if (x.empty()) {
 419                 return 1;
 420         }
 421         for (int i=1; i<n; i++) {
 422                 std::vector<cln::cl_N>::iterator be;
 423                 std::vector<cln::cl_N>::iterator en;
 424                 be = x.begin();
 425                 be++;
 426                 en = x.end();
 427                 std::vector<cln::cl_N> xbuf(be, en);
 428                 be = m.begin();
 429                 be++;
 430                 en = m.end();
 431                 std::vector<cln::cl_N> mbuf(be, en);
 432                 res = res + cln::expt(x[0],i) / cln::expt(i,m[0]) * numeric_zsum(i, xbuf, mbuf);
 433         }
 434         return res;
 435 }
 436
 437
 438 // helper function for harmonic polylogarithm
 439 static cln::cl_N numeric_harmonic(int n, std::vector<cln::cl_N>& m)
 440 {
 441         cln::cl_N res;
 442         if (m.empty()) {
 443                 return 1;
 444         }
 445         for (int i=1; i<n; i++) {
 446                 std::vector<cln::cl_N>::iterator be;
 447                 std::vector<cln::cl_N>::iterator en;
 448                 be = m.begin();
 449                 be++;
 450                 en = m.end();
 451                 std::vector<cln::cl_N> mbuf(be, en);
 452                 res = res + cln::recip(cln::expt(i,m[0])) * numeric_harmonic(i, mbuf);
 453         }
 454         return res;
 455 }
 456
 457
 458 /////////////////////////////
 459 // end of helper functions //
 460 /////////////////////////////
 461
 462
 463 // Polylogarithm and multiple polylogarithm
 464
 465 static ex Li_eval(const ex& x1, const ex& x2)
 466 {
 467         if (x2.is_zero()) {
 468                 return 0;
 469         }
 470         else {
 471                 return Li(x1,x2).hold();
 472         }
 473 }
 474
 475 static ex Li_evalf(const ex& x1, const ex& x2)
 476 {
 477         // classical polylogs
 478         if (is_a<numeric>(x1) && is_a<numeric>(x2)) {
 479                 return Li_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2));
 480         }
 481         // multiple polylogs
 482         else if (is_a<lst>(x1) && is_a<lst>(x2)) {
 483                 for (int i=0; i<x1.nops(); i++) {
 484                         if (!is_a<numeric>(x1.op(i)))
 485                                 return Li(x1,x2).hold();
 486                         if (!is_a<numeric>(x2.op(i)))
 487                                 return Li(x1,x2).hold();
 488                         if (x2.op(i) >= 1)
 489                                 return Li(x1,x2).hold();
 490                 }
 491
 492                 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
 493                 cln::cl_N x_1 = ex_to<numeric>(x2.op(x2.nops()-1)).to_cl_N();
 494                 std::vector<cln::cl_N> x;
 495                 std::vector<cln::cl_N> m;
 496                 const int nops = ex_to<numeric>(x1.nops()).to_int();
 497                 for (int i=nops-2; i>=0; i--) {
 498                         m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
 499                         x.push_back(ex_to<numeric>(x2.op(i)).to_cl_N());
 500                 }
 501
 502                 cln::cl_N res;
 503                 cln::cl_N resbuf;
 504                 for (int i=nops; true; i++) {
 505                         resbuf = res;
 506                         res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_zsum(i, x, m);
 507                         if (cln::zerop(res-resbuf))
 508                                 break;
 509                 }
 510
 511                 return numeric(res);
 512
 513         }
 514
 515         return Li(x1,x2).hold();
 516 }
 517
 518 static ex Li_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
 519 {
 520         epvector seq;
 521         seq.push_back(expair(Li(x1,x2), 0));
 522         return pseries(rel,seq);
 523 }
 524
 525 REGISTER_FUNCTION(Li, eval_func(Li_eval).evalf_func(Li_evalf).do_not_evalf_params().series_func(Li_series));
 526
 527
 528 // Nielsen's generalized polylogarithm
 529
 530 static ex S_eval(const ex& x1, const ex& x2, const ex& x3)
 531 {
 532         if (x2 == 1) {
 533                 return Li(x1+1,x3);
 534         }
 535         return S(x1,x2,x3).hold();
 536 }
 537
 538 static ex S_evalf(const ex& x1, const ex& x2, const ex& x3)
 539 {
 540         if (is_a<numeric>(x1) && is_a<numeric>(x2) && is_a<numeric>(x3)) {
 541                 if ((x3 == -1) && (x2 != 1)) {
 542                         // no formula to evaluate this ... sorry
 543                         return S(x1,x2,x3).hold();
 544                 }
 545                 return S_num(ex_to<numeric>(x1).to_int(), ex_to<numeric>(x2).to_int(), ex_to<numeric>(x3));
 546         }
 547         return S(x1,x2,x3).hold();
 548 }
 549
 550 static ex S_series(const ex& x1, const ex& x2, const ex& x3, const relational& rel, int order, unsigned options)
 551 {
 552         epvector seq;
 553         seq.push_back(expair(S(x1,x2,x3), 0));
 554         return pseries(rel,seq);
 555 }
 556
 557 REGISTER_FUNCTION(S, eval_func(S_eval).evalf_func(S_evalf).do_not_evalf_params().series_func(S_series));
 558
 559
 560 // Harmonic polylogarithm
 561
 562 static ex H_eval(const ex& x1, const ex& x2)
 563 {
 564         return H(x1,x2).hold();
 565 }
 566
 567 static ex H_evalf(const ex& x1, const ex& x2)
 568 {
 569         if (is_a<lst>(x1) && is_a<numeric>(x2)) {
 570                 for (int i=0; i<x1.nops(); i++) {
 571                         if (!is_a<numeric>(x1.op(i)))
 572                                 return H(x1,x2).hold();
 573                 }
 574                 if (x2 >= 1) {
 575                         return H(x1,x2).hold();
 576                 }
 577
 578                 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
 579                 cln::cl_N x_1 = ex_to<numeric>(x2).to_cl_N();
 580                 std::vector<cln::cl_N> m;
 581                 const int nops = ex_to<numeric>(x1.nops()).to_int();
 582                 for (int i=nops-2; i>=0; i--) {
 583                         m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
 584                 }
 585
 586                 cln::cl_N res;
 587                 cln::cl_N resbuf;
 588                 for (int i=nops; true; i++) {
 589                         resbuf = res;
 590                         res = res + cln::expt(x_1,i) / cln::expt(i,m_1) * numeric_harmonic(i, m);
 591                         if (cln::zerop(res-resbuf))
 592                                 break;
 593                 }
 594
 595                 return numeric(res);
 596
 597         }
 598
 599         return H(x1,x2).hold();
 600 }
 601
 602 static ex H_series(const ex& x1, const ex& x2, const relational& rel, int order, unsigned options)
 603 {
 604         epvector seq;
 605         seq.push_back(expair(H(x1,x2), 0));
 606         return pseries(rel,seq);
 607 }
 608
 609 REGISTER_FUNCTION(H, eval_func(H_eval).evalf_func(H_evalf).do_not_evalf_params().series_func(H_series));
 610
 611
 612 // Multiple zeta value
 613
 614 static ex mZeta_eval(const ex& x1)
 615 {
 616         return mZeta(x1).hold();
 617 }
 618
 619 static ex mZeta_evalf(const ex& x1)
 620 {
 621         if (is_a<lst>(x1)) {
 622                 for (int i=0; i<x1.nops(); i++) {
 623                         if (!is_a<numeric>(x1.op(i)))
 624                                 return mZeta(x1).hold();
 625                 }
 626
 627                 cln::cl_N m_1 = ex_to<numeric>(x1.op(x1.nops()-1)).to_cl_N();
 628                 std::vector<cln::cl_N> m;
 629                 const int nops = ex_to<numeric>(x1.nops()).to_int();
 630                 for (int i=nops-2; i>=0; i--) {
 631                         m.push_back(ex_to<numeric>(x1.op(i)).to_cl_N());
 632                 }
 633
 634                 cln::float_format_t prec = cln::default_float_format;
 635                 cln::cl_N res = cln::complex(cln::cl_float(0, prec), 0);
 636                 cln::cl_N resbuf;
 637                 for (int i=nops; true; i++) {
 638                         // to infinity and beyond ... timewise
 639                         resbuf = res;
 640                         res = res + cln::recip(cln::expt(i,m_1)) * numeric_harmonic(i, m);
 641                         if (cln::zerop(res-resbuf))
 642                                 break;
 643                 }
 644
 645                 return numeric(res);
 646
 647         }
 648
 649         return mZeta(x1).hold();
 650 }
 651
 652 static ex mZeta_series(const ex& x1, const relational& rel, int order, unsigned options)
 653 {
 654         epvector seq;
 655         seq.push_back(expair(mZeta(x1), 0));
 656         return pseries(rel,seq);
 657 }
 658
 659 REGISTER_FUNCTION(mZeta, eval_func(mZeta_eval).evalf_func(mZeta_evalf).do_not_evalf_params().series_func(mZeta_series));
 660
 661
 662 } // namespace GiNaC
 663