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

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.