GiNaC  1.8.3
inifcns_nstdsums.cpp
Go to the documentation of this file.
1 
49 /*
50  * GiNaC Copyright (C) 1999-2022 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  // 24.03.2021: this block can be outside the loop over r
1017  cln::cl_RA p(2);
1018  bool adjustp;
1019  do {
1020  adjustp = false;
1021  for (std::size_t i = 0; i < size; ++i) {
1022  // 24.03.2021: replaced (x[i] == cln::cl_RA(1)/p) by (cln::zerop(x[i] - cln::cl_RA(1)/p)
1023  // in the case where we compare a float with a rational, CLN behaves differently in the two versions
1024  if (cln::zerop(x[i] - cln::cl_RA(1)/p) ) {
1025  p = p/2 + cln::cl_RA(3)/2;
1026  adjustp = true;
1027  continue;
1028  }
1029  }
1030  } while (adjustp);
1031  cln::cl_RA q = p/(p-1);
1032 
1033  for (std::size_t r = 0; r <= size; ++r) {
1034  cln::cl_N buffer(1 & r ? -1 : 1);
1035  std::vector<cln::cl_N> qlstx;
1036  std::vector<int> qlsts;
1037  for (std::size_t j = r; j >= 1; --j) {
1038  qlstx.push_back(cln::cl_N(1) - x[j-1]);
1039  if (imagpart(x[j-1])==0 && realpart(x[j-1]) >= 1) {
1040  qlsts.push_back(1);
1041  } else {
1042  qlsts.push_back(-s[j-1]);
1043  }
1044  }
1045  if (qlstx.size() > 0) {
1046  buffer = buffer*G_numeric(qlstx, qlsts, 1/q);
1047  }
1048  std::vector<cln::cl_N> plstx;
1049  std::vector<int> plsts;
1050  for (std::size_t j = r+1; j <= size; ++j) {
1051  plstx.push_back(x[j-1]);
1052  plsts.push_back(s[j-1]);
1053  }
1054  if (plstx.size() > 0) {
1055  buffer = buffer*G_numeric(plstx, plsts, 1/p);
1056  }
1057  result = result + buffer;
1058  }
1059  return result;
1060 }
1061 
1062 class less_object_for_cl_N
1063 {
1064 public:
1065  bool operator() (const cln::cl_N & a, const cln::cl_N & b) const
1066  {
1067  // absolute value?
1068  if (abs(a) != abs(b))
1069  return (abs(a) < abs(b)) ? true : false;
1070 
1071  // complex phase?
1072  if (phase(a) != phase(b))
1073  return (phase(a) < phase(b)) ? true : false;
1074 
1075  // equal, therefore "less" is not true
1076  return false;
1077  }
1078 };
1079 
1080 
1081 // convergence transformation, used for numerical evaluation of G function.
1082 // the parameter x, s and y must only contain numerics
1083 static cln::cl_N
1084 G_do_trafo(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
1085  const cln::cl_N& y, bool flag_trailing_zeros_only)
1086 {
1087  // sort (|x|<->position) to determine indices
1088  typedef std::multimap<cln::cl_N, std::size_t, less_object_for_cl_N> sortmap_t;
1089  sortmap_t sortmap;
1090  std::size_t size = 0;
1091  for (std::size_t i = 0; i < x.size(); ++i) {
1092  if (!zerop(x[i])) {
1093  sortmap.insert(std::make_pair(x[i], i));
1094  ++size;
1095  }
1096  }
1097  // include upper limit (scale)
1098  sortmap.insert(std::make_pair(y, x.size()));
1099 
1100  // generate missing dummy-symbols
1101  int i = 1;
1102  // holding dummy-symbols for the G/Li transformations
1103  exvector gsyms;
1104  gsyms.push_back(symbol("GSYMS_ERROR"));
1105  cln::cl_N lastentry(0);
1106  for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
1107  if (it != sortmap.begin()) {
1108  if (it->second < x.size()) {
1109  if (x[it->second] == lastentry) {
1110  gsyms.push_back(gsyms.back());
1111  continue;
1112  }
1113  } else {
1114  if (y == lastentry) {
1115  gsyms.push_back(gsyms.back());
1116  continue;
1117  }
1118  }
1119  }
1120  std::ostringstream os;
1121  os << "a" << i;
1122  gsyms.push_back(symbol(os.str()));
1123  ++i;
1124  if (it->second < x.size()) {
1125  lastentry = x[it->second];
1126  } else {
1127  lastentry = y;
1128  }
1129  }
1130 
1131  // fill position data according to sorted indices and prepare substitution list
1132  Gparameter a(x.size());
1133  exmap subslst;
1134  std::size_t pos = 1;
1135  int scale = pos;
1136  for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
1137  if (it->second < x.size()) {
1138  if (s[it->second] > 0) {
1139  a[it->second] = pos;
1140  } else {
1141  a[it->second] = -int(pos);
1142  }
1143  subslst[gsyms[pos]] = numeric(x[it->second]);
1144  } else {
1145  scale = pos;
1146  subslst[gsyms[pos]] = numeric(y);
1147  }
1148  ++pos;
1149  }
1150 
1151  // do transformation
1152  Gparameter pendint;
1153  ex result = G_transform(pendint, a, scale, gsyms, flag_trailing_zeros_only);
1154  // replace dummy symbols with their values
1155  result = result.expand();
1156  result = result.subs(subslst).evalf();
1157  if (!is_a<numeric>(result))
1158  throw std::logic_error("G_do_trafo: G_transform returned non-numeric result");
1159 
1160  cln::cl_N ret = ex_to<numeric>(result).to_cl_N();
1161  return ret;
1162 }
1163 
1164 // handles the transformations and the numerical evaluation of G
1165 // the parameter x, s and y must only contain numerics
1166 static cln::cl_N
1167 G_numeric(const std::vector<cln::cl_N>& x, const std::vector<int>& s,
1168  const cln::cl_N& y)
1169 {
1170  // check for convergence and necessary accelerations
1171  bool need_trafo = false;
1172  bool need_hoelder = false;
1173  bool have_trailing_zero = false;
1174  std::size_t depth = 0;
1175  for (auto & xi : x) {
1176  if (!zerop(xi)) {
1177  ++depth;
1178  const cln::cl_N x_y = abs(xi) - y;
1179  if (instanceof(x_y, cln::cl_R_ring) &&
1180  realpart(x_y) < cln::least_negative_float(cln::float_format(Digits - 2)))
1181  need_trafo = true;
1182 
1183  if (abs(abs(xi/y) - 1) < 0.01)
1184  need_hoelder = true;
1185  }
1186  }
1187  if (zerop(x.back())) {
1188  have_trailing_zero = true;
1189  need_trafo = true;
1190  }
1191 
1192  if (depth == 1 && x.size() == 2 && !need_trafo)
1193  return - Li_projection(2, y/x[1], cln::float_format(Digits));
1194 
1195  // do acceleration transformation (hoelder convolution [BBB])
1196  if (need_hoelder && !have_trailing_zero)
1197  return G_do_hoelder(x, s, y);
1198 
1199  // convergence transformation
1200  if (need_trafo)
1201  return G_do_trafo(x, s, y, have_trailing_zero);
1202 
1203  // do summation
1204  std::vector<cln::cl_N> newx;
1205  newx.reserve(x.size());
1206  std::vector<int> m;
1207  m.reserve(x.size());
1208  int mcount = 1;
1209  int sign = 1;
1210  cln::cl_N factor = y;
1211  for (auto & xi : x) {
1212  if (zerop(xi)) {
1213  ++mcount;
1214  } else {
1215  newx.push_back(factor/xi);
1216  factor = xi;
1217  m.push_back(mcount);
1218  mcount = 1;
1219  sign = -sign;
1220  }
1221  }
1222 
1223  return sign*multipleLi_do_sum(m, newx);
1224 }
1225 
1226 
1227 ex mLi_numeric(const lst& m, const lst& x)
1228 {
1229  // let G_numeric do the transformation
1230  std::vector<cln::cl_N> newx;
1231  newx.reserve(x.nops());
1232  std::vector<int> s;
1233  s.reserve(x.nops());
1234  cln::cl_N factor(1);
1235  for (auto itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
1236  for (int i = 1; i < *itm; ++i) {
1237  newx.push_back(cln::cl_N(0));
1238  s.push_back(1);
1239  }
1240  const cln::cl_N xi = ex_to<numeric>(*itx).to_cl_N();
1241  factor = factor/xi;
1242  newx.push_back(factor);
1243  if ( !instanceof(factor, cln::cl_R_ring) && imagpart(factor) < 0 ) {
1244  s.push_back(-1);
1245  }
1246  else {
1247  s.push_back(1);
1248  }
1249  }
1250  return numeric(cln::cl_N(1 & m.nops() ? - 1 : 1)*G_numeric(newx, s, cln::cl_N(1)));
1251 }
1252 
1253 
1254 } // end of anonymous namespace
1255 
1256 
1258 //
1259 // Generalized multiple polylogarithm G(x, y) and G(x, s, y)
1260 //
1261 // GiNaC function
1262 //
1264 
1265 
1266 static ex G2_evalf(const ex& x_, const ex& y)
1267 {
1268  if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
1269  return G(x_, y).hold();
1270  }
1271  lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
1272  if (x.nops() == 0) {
1273  return _ex1;
1274  }
1275  if (x.op(0) == y) {
1276  return G(x_, y).hold();
1277  }
1278  std::vector<int> s;
1279  s.reserve(x.nops());
1280  bool all_zero = true;
1281  for (const auto & it : x) {
1282  if (!it.info(info_flags::numeric)) {
1283  return G(x_, y).hold();
1284  }
1285  if (it != _ex0) {
1286  all_zero = false;
1287  }
1288  if ( !ex_to<numeric>(it).is_real() && ex_to<numeric>(it).imag() < 0 ) {
1289  s.push_back(-1);
1290  }
1291  else {
1292  s.push_back(1);
1293  }
1294  }
1295  if (all_zero) {
1296  return pow(log(y), x.nops()) / factorial(x.nops());
1297  }
1298  std::vector<cln::cl_N> xv;
1299  xv.reserve(x.nops());
1300  for (const auto & it : x)
1301  xv.push_back(ex_to<numeric>(it).to_cl_N());
1302  cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
1303  return numeric(result);
1304 }
1305 
1306 
1307 static ex G2_eval(const ex& x_, const ex& y)
1308 {
1309  //TODO eval to MZV or H or S or Lin
1310 
1311  if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
1312  return G(x_, y).hold();
1313  }
1314  lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
1315  if (x.nops() == 0) {
1316  return _ex1;
1317  }
1318  if (x.op(0) == y) {
1319  return G(x_, y).hold();
1320  }
1321  std::vector<int> s;
1322  s.reserve(x.nops());
1323  bool all_zero = true;
1324  bool crational = true;
1325  for (const auto & it : x) {
1326  if (!it.info(info_flags::numeric)) {
1327  return G(x_, y).hold();
1328  }
1329  if (!it.info(info_flags::crational)) {
1330  crational = false;
1331  }
1332  if (it != _ex0) {
1333  all_zero = false;
1334  }
1335  if ( !ex_to<numeric>(it).is_real() && ex_to<numeric>(it).imag() < 0 ) {
1336  s.push_back(-1);
1337  }
1338  else {
1339  s.push_back(+1);
1340  }
1341  }
1342  if (all_zero) {
1343  return pow(log(y), x.nops()) / factorial(x.nops());
1344  }
1345  if (!y.info(info_flags::crational)) {
1346  crational = false;
1347  }
1348  if (crational) {
1349  return G(x_, y).hold();
1350  }
1351  std::vector<cln::cl_N> xv;
1352  xv.reserve(x.nops());
1353  for (const auto & it : x)
1354  xv.push_back(ex_to<numeric>(it).to_cl_N());
1355  cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
1356  return numeric(result);
1357 }
1358 
1359 
1360 // option do_not_evalf_params() removed.
1361 unsigned G2_SERIAL::serial = function::register_new(function_options("G", 2).
1362  evalf_func(G2_evalf).
1363  eval_func(G2_eval).
1364  overloaded(2));
1365 //TODO
1366 // derivative_func(G2_deriv).
1367 // print_func<print_latex>(G2_print_latex).
1368 
1369 
1370 static ex G3_evalf(const ex& x_, const ex& s_, const ex& y)
1371 {
1372  if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
1373  return G(x_, s_, y).hold();
1374  }
1375  lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
1376  lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst{s_};
1377  if (x.nops() != s.nops()) {
1378  return G(x_, s_, y).hold();
1379  }
1380  if (x.nops() == 0) {
1381  return _ex1;
1382  }
1383  if (x.op(0) == y) {
1384  return G(x_, s_, y).hold();
1385  }
1386  std::vector<int> sn;
1387  sn.reserve(s.nops());
1388  bool all_zero = true;
1389  for (auto itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
1390  if (!(*itx).info(info_flags::numeric)) {
1391  return G(x_, y).hold();
1392  }
1393  if (!(*its).info(info_flags::real)) {
1394  return G(x_, y).hold();
1395  }
1396  if (*itx != _ex0) {
1397  all_zero = false;
1398  }
1399  if ( ex_to<numeric>(*itx).is_real() ) {
1400  if ( ex_to<numeric>(*itx).is_positive() ) {
1401  if ( *its >= 0 ) {
1402  sn.push_back(1);
1403  }
1404  else {
1405  sn.push_back(-1);
1406  }
1407  } else {
1408  sn.push_back(1);
1409  }
1410  }
1411  else {
1412  if ( ex_to<numeric>(*itx).imag() > 0 ) {
1413  sn.push_back(1);
1414  }
1415  else {
1416  sn.push_back(-1);
1417  }
1418  }
1419  }
1420  if (all_zero) {
1421  return pow(log(y), x.nops()) / factorial(x.nops());
1422  }
1423  std::vector<cln::cl_N> xn;
1424  xn.reserve(x.nops());
1425  for (const auto & it : x)
1426  xn.push_back(ex_to<numeric>(it).to_cl_N());
1427  cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
1428  return numeric(result);
1429 }
1430 
1431 
1432 static ex G3_eval(const ex& x_, const ex& s_, const ex& y)
1433 {
1434  //TODO eval to MZV or H or S or Lin
1435 
1436  if ((!y.info(info_flags::numeric)) || (!y.info(info_flags::positive))) {
1437  return G(x_, s_, y).hold();
1438  }
1439  lst x = is_a<lst>(x_) ? ex_to<lst>(x_) : lst{x_};
1440  lst s = is_a<lst>(s_) ? ex_to<lst>(s_) : lst{s_};
1441  if (x.nops() != s.nops()) {
1442  return G(x_, s_, y).hold();
1443  }
1444  if (x.nops() == 0) {
1445  return _ex1;
1446  }
1447  if (x.op(0) == y) {
1448  return G(x_, s_, y).hold();
1449  }
1450  std::vector<int> sn;
1451  sn.reserve(s.nops());
1452  bool all_zero = true;
1453  bool crational = true;
1454  for (auto itx = x.begin(), its = s.begin(); itx != x.end(); ++itx, ++its) {
1455  if (!(*itx).info(info_flags::numeric)) {
1456  return G(x_, s_, y).hold();
1457  }
1458  if (!(*its).info(info_flags::real)) {
1459  return G(x_, s_, y).hold();
1460  }
1461  if (!(*itx).info(info_flags::crational)) {
1462  crational = false;
1463  }
1464  if (*itx != _ex0) {
1465  all_zero = false;
1466  }
1467  if ( ex_to<numeric>(*itx).is_real() ) {
1468  if ( ex_to<numeric>(*itx).is_positive() ) {
1469  if ( *its >= 0 ) {
1470  sn.push_back(1);
1471  }
1472  else {
1473  sn.push_back(-1);
1474  }
1475  } else {
1476  sn.push_back(1);
1477  }
1478  }
1479  else {
1480  if ( ex_to<numeric>(*itx).imag() > 0 ) {
1481  sn.push_back(1);
1482  }
1483  else {
1484  sn.push_back(-1);
1485  }
1486  }
1487  }
1488  if (all_zero) {
1489  return pow(log(y), x.nops()) / factorial(x.nops());
1490  }
1491  if (!y.info(info_flags::crational)) {
1492  crational = false;
1493  }
1494  if (crational) {
1495  return G(x_, s_, y).hold();
1496  }
1497  std::vector<cln::cl_N> xn;
1498  xn.reserve(x.nops());
1499  for (const auto & it : x)
1500  xn.push_back(ex_to<numeric>(it).to_cl_N());
1501  cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
1502  return numeric(result);
1503 }
1504 
1505 
1506 // option do_not_evalf_params() removed.
1507 // This is safe: in the code above it only matters if s_ > 0 or s_ < 0,
1508 // s_ is allowed to be of floating type.
1509 unsigned G3_SERIAL::serial = function::register_new(function_options("G", 3).
1510  evalf_func(G3_evalf).
1511  eval_func(G3_eval).
1512  overloaded(2));
1513 //TODO
1514 // derivative_func(G3_deriv).
1515 // print_func<print_latex>(G3_print_latex).
1516 
1517 
1519 //
1520 // Classical polylogarithm and multiple polylogarithm Li(m,x)
1521 //
1522 // GiNaC function
1523 //
1525 
1526 
1527 static ex Li_evalf(const ex& m_, const ex& x_)
1528 {
1529  // classical polylogs
1530  if (m_.info(info_flags::posint)) {
1531  if (x_.info(info_flags::numeric)) {
1532  int m__ = ex_to<numeric>(m_).to_int();
1533  const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
1534  const cln::cl_N result = Lin_numeric(m__, x__);
1535  return numeric(result);
1536  } else {
1537  // try to numerically evaluate second argument
1538  ex x_val = x_.evalf();
1539  if (x_val.info(info_flags::numeric)) {
1540  int m__ = ex_to<numeric>(m_).to_int();
1541  const cln::cl_N x__ = ex_to<numeric>(x_val).to_cl_N();
1542  const cln::cl_N result = Lin_numeric(m__, x__);
1543  return numeric(result);
1544  }
1545  }
1546  }
1547  // multiple polylogs
1548  if (is_a<lst>(m_) && is_a<lst>(x_)) {
1549 
1550  const lst& m = ex_to<lst>(m_);
1551  const lst& x = ex_to<lst>(x_);
1552  if (m.nops() != x.nops()) {
1553  return Li(m_,x_).hold();
1554  }
1555  if (x.nops() == 0) {
1556  return _ex1;
1557  }
1558  if ((m.op(0) == _ex1) && (x.op(0) == _ex1)) {
1559  return Li(m_,x_).hold();
1560  }
1561 
1562  for (auto itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
1563  if (!(*itm).info(info_flags::posint)) {
1564  return Li(m_, x_).hold();
1565  }
1566  if (!(*itx).info(info_flags::numeric)) {
1567  return Li(m_, x_).hold();
1568  }
1569  if (*itx == _ex0) {
1570  return _ex0;
1571  }
1572  }
1573 
1574  return mLi_numeric(m, x);
1575  }
1576 
1577  return Li(m_,x_).hold();
1578 }
1579 
1580 
1581 static ex Li_eval(const ex& m_, const ex& x_)
1582 {
1583  if (is_a<lst>(m_)) {
1584  if (is_a<lst>(x_)) {
1585  // multiple polylogs
1586  const lst& m = ex_to<lst>(m_);
1587  const lst& x = ex_to<lst>(x_);
1588  if (m.nops() != x.nops()) {
1589  return Li(m_,x_).hold();
1590  }
1591  if (x.nops() == 0) {
1592  return _ex1;
1593  }
1594  bool is_H = true;
1595  bool is_zeta = true;
1596  bool do_evalf = true;
1597  bool crational = true;
1598  for (auto itm = m.begin(), itx = x.begin(); itm != m.end(); ++itm, ++itx) {
1599  if (!(*itm).info(info_flags::posint)) {
1600  return Li(m_,x_).hold();
1601  }
1602  if ((*itx != _ex1) && (*itx != _ex_1)) {
1603  if (itx != x.begin()) {
1604  is_H = false;
1605  }
1606  is_zeta = false;
1607  }
1608  if (*itx == _ex0) {
1609  return _ex0;
1610  }
1611  if (!(*itx).info(info_flags::numeric)) {
1612  do_evalf = false;
1613  }
1614  if (!(*itx).info(info_flags::crational)) {
1615  crational = false;
1616  }
1617  }
1618  if (is_zeta) {
1619  lst newx;
1620  for (const auto & itx : x) {
1621  GINAC_ASSERT((itx == _ex1) || (itx == _ex_1));
1622  // XXX: 1 + 0.0*I is considered equal to 1. However
1623  // the former is a not automatically converted
1624  // to a real number. Do the conversion explicitly
1625  // to avoid the "numeric::operator>(): complex inequality"
1626  // exception (and similar problems).
1627  newx.append(itx != _ex_1 ? _ex1 : _ex_1);
1628  }
1629  return zeta(m_, newx);
1630  }
1631  if (is_H) {
1632  ex prefactor;
1633  lst newm = convert_parameter_Li_to_H(m, x, prefactor);
1634  return prefactor * H(newm, x[0]);
1635  }
1636  if (do_evalf && !crational) {
1637  return mLi_numeric(m,x);
1638  }
1639  }
1640  return Li(m_, x_).hold();
1641  } else if (is_a<lst>(x_)) {
1642  return Li(m_, x_).hold();
1643  }
1644 
1645  // classical polylogs
1646  if (x_ == _ex0) {
1647  return _ex0;
1648  }
1649  if (x_ == _ex1) {
1650  return zeta(m_);
1651  }
1652  if (x_ == _ex_1) {
1653  return (pow(2,1-m_)-1) * zeta(m_);
1654  }
1655  if (m_ == _ex1) {
1656  return -log(1-x_);
1657  }
1658  if (m_ == _ex2) {
1659  if (x_.is_equal(I)) {
1660  return power(Pi,_ex2)/_ex_48 + Catalan*I;
1661  }
1662  if (x_.is_equal(-I)) {
1663  return power(Pi,_ex2)/_ex_48 - Catalan*I;
1664  }
1665  }
1667  int m__ = ex_to<numeric>(m_).to_int();
1668  const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
1669  const cln::cl_N result = Lin_numeric(m__, x__);
1670  return numeric(result);
1671  }
1672 
1673  return Li(m_, x_).hold();
1674 }
1675 
1676 
1677 static ex Li_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
1678 {
1679  if (is_a<lst>(m) || is_a<lst>(x)) {
1680  // multiple polylog
1681  epvector seq { expair(Li(m, x), 0) };
1682  return pseries(rel, std::move(seq));
1683  }
1684 
1685  // classical polylog
1686  const ex x_pt = x.subs(rel, subs_options::no_pattern);
1687  if (m.info(info_flags::numeric) && x_pt.info(info_flags::numeric)) {
1688  // First special case: x==0 (derivatives have poles)
1689  if (x_pt.is_zero()) {
1690  const symbol s;
1691  ex ser;
1692  // manually construct the primitive expansion
1693  for (int i=1; i<order; ++i)
1694  ser += pow(s,i) / pow(numeric(i), m);
1695  // substitute the argument's series expansion
1696  ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
1697  // maybe that was terminating, so add a proper order term
1698  epvector nseq { expair(Order(_ex1), order) };
1699  ser += pseries(rel, std::move(nseq));
1700  // reexpanding it will collapse the series again
1701  return ser.series(rel, order);
1702  }
1703  // TODO special cases: x==1 (branch point) and x real, >=1 (branch cut)
1704  throw std::runtime_error("Li_series: don't know how to do the series expansion at this point!");
1705  }
1706  // all other cases should be safe, by now:
1707  throw do_taylor(); // caught by function::series()
1708 }
1709 
1710 
1711 static ex Li_deriv(const ex& m_, const ex& x_, unsigned deriv_param)
1712 {
1713  GINAC_ASSERT(deriv_param < 2);
1714  if (deriv_param == 0) {
1715  return _ex0;
1716  }
1717  if (m_.nops() > 1) {
1718  throw std::runtime_error("don't know how to derivate multiple polylogarithm!");
1719  }
1720  ex m;
1721  if (is_a<lst>(m_)) {
1722  m = m_.op(0);
1723  } else {
1724  m = m_;
1725  }
1726  ex x;
1727  if (is_a<lst>(x_)) {
1728  x = x_.op(0);
1729  } else {
1730  x = x_;
1731  }
1732  if (m > 0) {
1733  return Li(m-1, x) / x;
1734  } else {
1735  return 1/(1-x);
1736  }
1737 }
1738 
1739 
1740 static void Li_print_latex(const ex& m_, const ex& x_, const print_context& c)
1741 {
1742  lst m;
1743  if (is_a<lst>(m_)) {
1744  m = ex_to<lst>(m_);
1745  } else {
1746  m = lst{m_};
1747  }
1748  lst x;
1749  if (is_a<lst>(x_)) {
1750  x = ex_to<lst>(x_);
1751  } else {
1752  x = lst{x_};
1753  }
1754  c.s << "\\mathrm{Li}_{";
1755  auto itm = m.begin();
1756  (*itm).print(c);
1757  itm++;
1758  for (; itm != m.end(); itm++) {
1759  c.s << ",";
1760  (*itm).print(c);
1761  }
1762  c.s << "}(";
1763  auto itx = x.begin();
1764  (*itx).print(c);
1765  itx++;
1766  for (; itx != x.end(); itx++) {
1767  c.s << ",";
1768  (*itx).print(c);
1769  }
1770  c.s << ")";
1771 }
1772 
1773 
1775  evalf_func(Li_evalf).
1776  eval_func(Li_eval).
1777  series_func(Li_series).
1778  derivative_func(Li_deriv).
1779  print_func<print_latex>(Li_print_latex).
1780  do_not_evalf_params());
1781 
1782 
1784 //
1785 // Nielsen's generalized polylogarithm S(n,p,x)
1786 //
1787 // helper functions
1788 //
1790 
1791 
1792 // anonymous namespace for helper functions
1793 namespace {
1794 
1795 
1796 // lookup table for special Euler-Zagier-Sums (used for S_n,p(x))
1797 // see fill_Yn()
1798 std::vector<std::vector<cln::cl_N>> Yn;
1799 int ynsize = 0; // number of Yn[]
1800 int ynlength = 100; // initial length of all Yn[i]
1801 
1802 
1803 // This function calculates the Y_n. The Y_n are needed for the evaluation of S_{n,p}(x).
1804 // The Y_n are basically Euler-Zagier sums with all m_i=1. They are subsums in the Z-sum
1805 // representing S_{n,p}(x).
1806 // The first index in Y_n corresponds to the parameter p minus one, i.e. the depth of the
1807 // equivalent Z-sum.
1808 // The second index in Y_n corresponds to the running index of the outermost sum in the full Z-sum
1809 // representing S_{n,p}(x).
1810 // The calculation of Y_n uses the values from Y_{n-1}.
1811 void fill_Yn(int n, const cln::float_format_t& prec)
1812 {
1813  const int initsize = ynlength;
1814  //const int initsize = initsize_Yn;
1815  cln::cl_N one = cln::cl_float(1, prec);
1816 
1817  if (n) {
1818  std::vector<cln::cl_N> buf(initsize);
1819  auto it = buf.begin();
1820  auto itprev = Yn[n-1].begin();
1821  *it = (*itprev) / cln::cl_N(n+1) * one;
1822  it++;
1823  itprev++;
1824  // sums with an index smaller than the depth are zero and need not to be calculated.
1825  // calculation starts with depth, which is n+2)
1826  for (int i=n+2; i<=initsize+n; i++) {
1827  *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
1828  it++;
1829  itprev++;
1830  }
1831  Yn.push_back(buf);
1832  } else {
1833  std::vector<cln::cl_N> buf(initsize);
1834  auto it = buf.begin();
1835  *it = 1 * one;
1836  it++;
1837  for (int i=2; i<=initsize; i++) {
1838  *it = *(it-1) + 1 / cln::cl_N(i) * one;
1839  it++;
1840  }
1841  Yn.push_back(buf);
1842  }
1843  ynsize++;
1844 }
1845 
1846 
1847 // make Yn longer ...
1848 void make_Yn_longer(int newsize, const cln::float_format_t& prec)
1849 {
1850 
1851  cln::cl_N one = cln::cl_float(1, prec);
1852 
1853  Yn[0].resize(newsize);
1854  auto it = Yn[0].begin();
1855  it += ynlength;
1856  for (int i=ynlength+1; i<=newsize; i++) {
1857  *it = *(it-1) + 1 / cln::cl_N(i) * one;
1858  it++;
1859  }
1860 
1861  for (int n=1; n<ynsize; n++) {
1862  Yn[n].resize(newsize);
1863  auto it = Yn[n].begin();
1864  auto itprev = Yn[n-1].begin();
1865  it += ynlength;
1866  itprev += ynlength;
1867  for (int i=ynlength+n+1; i<=newsize+n; i++) {
1868  *it = *(it-1) + (*itprev) / cln::cl_N(i) * one;
1869  it++;
1870  itprev++;
1871  }
1872  }
1873 
1874  ynlength = newsize;
1875 }
1876 
1877 
1878 // helper function for S(n,p,x)
1879 // [Kol] (7.2)
1880 cln::cl_N C(int n, int p)
1881 {
1882  cln::cl_N result;
1883 
1884  for (int k=0; k<p; k++) {
1885  for (int j=0; j<=(n+k-1)/2; j++) {
1886  if (k == 0) {
1887  if (n & 1) {
1888  if (j & 1) {
1889  result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
1890  }
1891  else {
1892  result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(n-2*j,p,1) / cln::factorial(2*j);
1893  }
1894  }
1895  }
1896  else {
1897  if (k & 1) {
1898  if (j & 1) {
1899  result = result + cln::factorial(n+k-1)
1900  * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1901  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1902  }
1903  else {
1904  result = result - cln::factorial(n+k-1)
1905  * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1906  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1907  }
1908  }
1909  else {
1910  if (j & 1) {
1911  result = result - cln::factorial(n+k-1) * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1912  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1913  }
1914  else {
1915  result = result + cln::factorial(n+k-1)
1916  * cln::expt(cln::pi(),2*j) * S_num(n+k-2*j,p-k,1)
1917  / (cln::factorial(k) * cln::factorial(n-1) * cln::factorial(2*j));
1918  }
1919  }
1920  }
1921  }
1922  }
1923  int np = n+p;
1924  if ((np-1) & 1) {
1925  if (((np)/2+n) & 1) {
1926  result = -result - cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
1927  }
1928  else {
1929  result = -result + cln::expt(cln::pi(),np) / (np * cln::factorial(n-1) * cln::factorial(p));
1930  }
1931  }
1932 
1933  return result;
1934 }
1935 
1936 
1937 // helper function for S(n,p,x)
1938 // [Kol] remark to (9.1)
1939 cln::cl_N a_k(int k)
1940 {
1941  cln::cl_N result;
1942 
1943  if (k == 0) {
1944  return 1;
1945  }
1946 
1947  result = result;
1948  for (int m=2; m<=k; m++) {
1949  result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * a_k(k-m);
1950  }
1951 
1952  return -result / k;
1953 }
1954 
1955 
1956 // helper function for S(n,p,x)
1957 // [Kol] remark to (9.1)
1958 cln::cl_N b_k(int k)
1959 {
1960  cln::cl_N result;
1961 
1962  if (k == 0) {
1963  return 1;
1964  }
1965 
1966  result = result;
1967  for (int m=2; m<=k; m++) {
1968  result = result + cln::expt(cln::cl_N(-1),m) * cln::zeta(m) * b_k(k-m);
1969  }
1970 
1971  return result / k;
1972 }
1973 
1974 
1975 // helper function for S(n,p,x)
1976 cln::cl_N S_do_sum(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
1977 {
1978  static cln::float_format_t oldprec = cln::default_float_format;
1979 
1980  if (p==1) {
1981  return Li_projection(n+1, x, prec);
1982  }
1983 
1984  // precision has changed, we need to clear lookup table Yn
1985  if ( oldprec != prec ) {
1986  Yn.clear();
1987  ynsize = 0;
1988  ynlength = 100;
1989  oldprec = prec;
1990  }
1991 
1992  // check if precalculated values are sufficient
1993  if (p > ynsize+1) {
1994  for (int i=ynsize; i<p-1; i++) {
1995  fill_Yn(i, prec);
1996  }
1997  }
1998 
1999  // should be done otherwise
2000  cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
2001  cln::cl_N xf = x * one;
2002  //cln::cl_N xf = x * cln::cl_float(1, prec);
2003 
2004  cln::cl_N res;
2005  cln::cl_N resbuf;
2006  cln::cl_N factor = cln::expt(xf, p);
2007  int i = p;
2008  do {
2009  resbuf = res;
2010  if (i-p >= ynlength) {
2011  // make Yn longer
2012  make_Yn_longer(ynlength*2, prec);
2013  }
2014  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? ...
2015  //res = res + factor / cln::expt(cln::cl_I(i),n+1) * (*it); // should we check it? or rely on magic number? ...
2016  factor = factor * xf;
2017  i++;
2018  } while (res != resbuf);
2019 
2020  return res;
2021 }
2022 
2023 
2024 // helper function for S(n,p,x)
2025 cln::cl_N S_projection(int n, int p, const cln::cl_N& x, const cln::float_format_t& prec)
2026 {
2027  // [Kol] (5.3)
2028  if (cln::abs(cln::realpart(x)) > cln::cl_F("0.5")) {
2029 
2030  cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(x),n)
2031  * cln::expt(cln::log(1-x),p) / cln::factorial(n) / cln::factorial(p);
2032 
2033  for (int s=0; s<n; s++) {
2034  cln::cl_N res2;
2035  for (int r=0; r<p; r++) {
2036  res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-x),r)
2037  * S_do_sum(p-r,n-s,1-x,prec) / cln::factorial(r);
2038  }
2039  result = result + cln::expt(cln::log(x),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
2040  }
2041 
2042  return result;
2043  }
2044 
2045  return S_do_sum(n, p, x, prec);
2046 }
2047 
2048 
2049 // helper function for S(n,p,x)
2050 const cln::cl_N S_num(int n, int p, const cln::cl_N& x)
2051 {
2052  if (x == 1) {
2053  if (n == 1) {
2054  // [Kol] (2.22) with (2.21)
2055  return cln::zeta(p+1);
2056  }
2057 
2058  if (p == 1) {
2059  // [Kol] (2.22)
2060  return cln::zeta(n+1);
2061  }
2062 
2063  // [Kol] (9.1)
2064  cln::cl_N result;
2065  for (int nu=0; nu<n; nu++) {
2066  for (int rho=0; rho<=p; rho++) {
2067  result = result + b_k(n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
2068  * cln::factorial(nu+rho+1) / cln::factorial(rho) / cln::factorial(nu+1);
2069  }
2070  }
2071  result = result * cln::expt(cln::cl_I(-1),n+p-1);
2072 
2073  return result;
2074  }
2075  else if (x == -1) {
2076  // [Kol] (2.22)
2077  if (p == 1) {
2078  return -(1-cln::expt(cln::cl_I(2),-n)) * cln::zeta(n+1);
2079  }
2080 // throw std::runtime_error("don't know how to evaluate this function!");
2081  }
2082 
2083  // what is the desired float format?
2084  // first guess: default format
2085  cln::float_format_t prec = cln::default_float_format;
2086  const cln::cl_N value = x;
2087  // second guess: the argument's format
2088  if (!instanceof(realpart(value), cln::cl_RA_ring))
2089  prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(value)));
2090  else if (!instanceof(imagpart(value), cln::cl_RA_ring))
2091  prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(value)));
2092 
2093  // [Kol] (5.3)
2094  // 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
2095  // we don't care here about abs(value)<1 && real(value)>0.5, this will be taken care of in S_projection
2096  if ((cln::realpart(value) < -0.5) || (n == 0) || ((cln::abs(value) <= 1) && (cln::abs(value) > 0.95) && (cln::abs(1-value) > 1) )) {
2097 
2098  cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(cln::log(value),n)
2099  * cln::expt(cln::log(1-value),p) / cln::factorial(n) / cln::factorial(p);
2100 
2101  for (int s=0; s<n; s++) {
2102  cln::cl_N res2;
2103  for (int r=0; r<p; r++) {
2104  res2 = res2 + cln::expt(cln::cl_I(-1),r) * cln::expt(cln::log(1-value),r)
2105  * S_num(p-r,n-s,1-value) / cln::factorial(r);
2106  }
2107  result = result + cln::expt(cln::log(value),s) * (S_num(n-s,p,1) - res2) / cln::factorial(s);
2108  }
2109 
2110  return result;
2111 
2112  }
2113  // [Kol] (5.12)
2114  if (cln::abs(value) > 1) {
2115 
2116  cln::cl_N result;
2117 
2118  for (int s=0; s<p; s++) {
2119  for (int r=0; r<=s; r++) {
2120  result = result + cln::expt(cln::cl_I(-1),s) * cln::expt(cln::log(-value),r) * cln::factorial(n+s-r-1)
2122  * S_num(n+s-r,p-s,cln::recip(value));
2123  }
2124  }
2125  result = result * cln::expt(cln::cl_I(-1),n);
2126 
2127  cln::cl_N res2;
2128  for (int r=0; r<n; r++) {
2129  res2 = res2 + cln::expt(cln::log(-value),r) * C(n-r,p) / cln::factorial(r);
2130  }
2131  res2 = res2 + cln::expt(cln::log(-value),n+p) / cln::factorial(n+p);
2132 
2133  result = result + cln::expt(cln::cl_I(-1),p) * res2;
2134 
2135  return result;
2136  }
2137 
2138  if ((cln::abs(value) > 0.95) && (cln::abs(value-9.53) < 9.47)) {
2139  lst m;
2140  m.append(n+1);
2141  for (int s=0; s<p-1; s++)
2142  m.append(1);
2143 
2144  ex res = H(m,numeric(value)).evalf();
2145  return ex_to<numeric>(res).to_cl_N();
2146  }
2147  else {
2148  return S_projection(n, p, value, prec);
2149  }
2150 }
2151 
2152 
2153 } // end of anonymous namespace
2154 
2155 
2157 //
2158 // Nielsen's generalized polylogarithm S(n,p,x)
2159 //
2160 // GiNaC function
2161 //
2163 
2164 
2165 static ex S_evalf(const ex& n, const ex& p, const ex& x)
2166 {
2167  if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
2168  const int n_ = ex_to<numeric>(n).to_int();
2169  const int p_ = ex_to<numeric>(p).to_int();
2170  if (is_a<numeric>(x)) {
2171  const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
2172  const cln::cl_N result = S_num(n_, p_, x_);
2173  return numeric(result);
2174  } else {
2175  ex x_val = x.evalf();
2176  if (is_a<numeric>(x_val)) {
2177  const cln::cl_N x_val_ = ex_to<numeric>(x_val).to_cl_N();
2178  const cln::cl_N result = S_num(n_, p_, x_val_);
2179  return numeric(result);
2180  }
2181  }
2182  }
2183  return S(n, p, x).hold();
2184 }
2185 
2186 
2187 static ex S_eval(const ex& n, const ex& p, const ex& x)
2188 {
2189  if (n.info(info_flags::posint) && p.info(info_flags::posint)) {
2190  if (x == 0) {
2191  return _ex0;
2192  }
2193  if (x == 1) {
2194  lst m{n+1};
2195  for (int i=ex_to<numeric>(p).to_int()-1; i>0; i--) {
2196  m.append(1);
2197  }
2198  return zeta(m);
2199  }
2200  if (p == 1) {
2201  return Li(n+1, x);
2202  }
2204  int n_ = ex_to<numeric>(n).to_int();
2205  int p_ = ex_to<numeric>(p).to_int();
2206  const cln::cl_N x_ = ex_to<numeric>(x).to_cl_N();
2207  const cln::cl_N result = S_num(n_, p_, x_);
2208  return numeric(result);
2209  }
2210  }
2211  if (n.is_zero()) {
2212  // [Kol] (5.3)
2213  return pow(-log(1-x), p) / factorial(p);
2214  }
2215  return S(n, p, x).hold();
2216 }
2217 
2218 
2219 static ex S_series(const ex& n, const ex& p, const ex& x, const relational& rel, int order, unsigned options)
2220 {
2221  if (p == _ex1) {
2222  return Li(n+1, x).series(rel, order, options);
2223  }
2224 
2225  const ex x_pt = x.subs(rel, subs_options::no_pattern);
2227  // First special case: x==0 (derivatives have poles)
2228  if (x_pt.is_zero()) {
2229  const symbol s;
2230  ex ser;
2231  // manually construct the primitive expansion
2232  // subsum = Euler-Zagier-Sum is needed
2233  // dirty hack (slow ...) calculation of subsum:
2234  std::vector<ex> presubsum, subsum;
2235  subsum.push_back(0);
2236  for (int i=1; i<order-1; ++i) {
2237  subsum.push_back(subsum[i-1] + numeric(1, i));
2238  }
2239  for (int depth=2; depth<p; ++depth) {
2240  presubsum = subsum;
2241  for (int i=1; i<order-1; ++i) {
2242  subsum[i] = subsum[i-1] + numeric(1, i) * presubsum[i-1];
2243  }
2244  }
2245 
2246  for (int i=1; i<order; ++i) {
2247  ser += pow(s,i) / pow(numeric(i), n+1) * subsum[i-1];
2248  }
2249  // substitute the argument's series expansion
2250  ser = ser.subs(s==x.series(rel, order), subs_options::no_pattern);
2251  // maybe that was terminating, so add a proper order term
2252  epvector nseq { expair(Order(_ex1), order) };
2253  ser += pseries(rel, std::move(nseq));
2254  // reexpanding it will collapse the series again
2255  return ser.series(rel, order);
2256  }
2257  // TODO special cases: x==1 (branch point) and x real, >=1 (branch cut)
2258  throw std::runtime_error("S_series: don't know how to do the series expansion at this point!");
2259  }
2260  // all other cases should be safe, by now:
2261  throw do_taylor(); // caught by function::series()
2262 }
2263 
2264 
2265 static ex S_deriv(const ex& n, const ex& p, const ex& x, unsigned deriv_param)
2266 {
2267  GINAC_ASSERT(deriv_param < 3);
2268  if (deriv_param < 2) {
2269  return _ex0;
2270  }
2271  if (n > 0) {
2272  return S(n-1, p, x) / x;
2273  } else {
2274  return S(n, p-1, x) / (1-x);
2275  }
2276 }
2277 
2278 
2279 static void S_print_latex(const ex& n, const ex& p, const ex& x, const print_context& c)
2280 {
2281  c.s << "\\mathrm{S}_{";
2282  n.print(c);
2283  c.s << ",";
2284  p.print(c);
2285  c.s << "}(";
2286  x.print(c);
2287  c.s << ")";
2288 }
2289 
2290 
2292  evalf_func(S_evalf).
2293  eval_func(S_eval).
2294  series_func(S_series).
2295  derivative_func(S_deriv).
2296  print_func<print_latex>(S_print_latex).
2297  do_not_evalf_params());
2298 
2299 
2301 //
2302 // Harmonic polylogarithm H(m,x)
2303 //
2304 // helper functions
2305 //
2307 
2308 
2309 // anonymous namespace for helper functions
2310 namespace {
2311 
2312 
2313 // regulates the pole (used by 1/x-transformation)
2314 symbol H_polesign("IMSIGN");
2315 
2316 
2317 // convert parameters from H to Li representation
2318 // parameters are expected to be in expanded form, i.e. only 0, 1 and -1
2319 // returns true if some parameters are negative
2320 bool convert_parameter_H_to_Li(const lst& l, lst& m, lst& s, ex& pf)
2321 {
2322  // expand parameter list
2323  lst mexp;
2324  for (const auto & it : l) {
2325  if (it > 1) {
2326  for (ex count=it-1; count > 0; count--) {
2327  mexp.append(0);
2328  }
2329  mexp.append(1);
2330  } else if (it < -1) {
2331  for (ex count=it+1; count < 0; count++) {
2332  mexp.append(0);
2333  }
2334  mexp.append(-1);
2335  } else {
2336  mexp.append(it);
2337  }
2338  }
2339 
2340  ex signum = 1;
2341  pf = 1;
2342  bool has_negative_parameters = false;
2343  ex acc = 1;
2344  for (const auto & it : mexp) {
2345  if (it == 0) {
2346  acc++;
2347  continue;
2348  }
2349  if (it > 0) {
2350  m.append((it+acc-1) * signum);
2351  } else {
2352  m.append((it-acc+1) * signum);
2353  }
2354  acc = 1;
2355  signum = it;
2356  pf *= it;
2357  if (pf < 0) {
2358  has_negative_parameters = true;
2359  }
2360  }
2361  if (has_negative_parameters) {
2362  for (std::size_t i=0; i<m.nops(); i++) {
2363  if (m.op(i) < 0) {
2364  m.let_op(i) = -m.op(i);
2365  s.append(-1);
2366  } else {
2367  s.append(1);
2368  }
2369  }
2370  }
2371 
2372  return has_negative_parameters;
2373 }
2374 
2375 
2376 // recursivly transforms H to corresponding multiple polylogarithms
2377 struct map_trafo_H_convert_to_Li : public map_function
2378 {
2379  ex operator()(const ex& e) override
2380  {
2381  if (is_a<add>(e) || is_a<mul>(e)) {
2382  return e.map(*this);
2383  }
2384  if (is_a<function>(e)) {
2385  std::string name = ex_to<function>(e).get_name();
2386  if (name == "H") {
2387  lst parameter;
2388  if (is_a<lst>(e.op(0))) {
2389  parameter = ex_to<lst>(e.op(0));
2390  } else {
2391  parameter = lst{e.op(0)};
2392  }
2393  ex arg = e.op(1);
2394 
2395  lst m;
2396  lst s;
2397  ex pf;
2398  if (convert_parameter_H_to_Li(parameter, m, s, pf)) {
2399  s.let_op(0) = s.op(0) * arg;
2400  return pf * Li(m, s).hold();
2401  } else {
2402  for (std::size_t i=0; i<m.nops(); i++) {
2403  s.append(1);
2404  }
2405  s.let_op(0) = s.op(0) * arg;
2406  return Li(m, s).hold();
2407  }
2408  }
2409  }
2410  return e;
2411  }
2412 };
2413 
2414 
2415 // recursivly transforms H to corresponding zetas
2416 struct map_trafo_H_convert_to_zeta : public map_function
2417 {
2418  ex operator()(const ex& e) override
2419  {
2420  if (is_a<add>(e) || is_a<mul>(e)) {
2421  return e.map(*this);
2422  }
2423  if (is_a<function>(e)) {
2424  std::string name = ex_to<function>(e).get_name();
2425  if (name == "H") {
2426  lst parameter;
2427  if (is_a<lst>(e.op(0))) {
2428  parameter = ex_to<lst>(e.op(0));
2429  } else {
2430  parameter = lst{e.op(0)};
2431  }
2432 
2433  lst m;
2434  lst s;
2435  ex pf;
2436  if (convert_parameter_H_to_Li(parameter, m, s, pf)) {
2437  return pf * zeta(m, s);
2438  } else {
2439  return zeta(m);
2440  }
2441  }
2442  }
2443  return e;
2444  }
2445 };
2446 
2447 
2448 // remove trailing zeros from H-parameters
2449 struct map_trafo_H_reduce_trailing_zeros : public map_function
2450 {
2451  ex operator()(const ex& e) override
2452  {
2453  if (is_a<add>(e) || is_a<mul>(e)) {
2454  return e.map(*this);
2455  }
2456  if (is_a<function>(e)) {
2457  std::string name = ex_to<function>(e).get_name();
2458  if (name == "H") {
2459  lst parameter;
2460  if (is_a<lst>(e.op(0))) {
2461  parameter = ex_to<lst>(e.op(0));
2462  } else {
2463  parameter = lst{e.op(0)};
2464  }
2465  ex arg = e.op(1);
2466  if (parameter.op(parameter.nops()-1) == 0) {
2467 
2468  //
2469  if (parameter.nops() == 1) {
2470  return log(arg);
2471  }
2472 
2473  //
2474  auto it = parameter.begin();
2475  while ((it != parameter.end()) && (*it == 0)) {
2476  it++;
2477  }
2478  if (it == parameter.end()) {
2479  return pow(log(arg),parameter.nops()) / factorial(parameter.nops());
2480  }
2481 
2482  //
2483  parameter.remove_last();
2484  std::size_t lastentry = parameter.nops();
2485  while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
2486  lastentry--;
2487  }
2488 
2489  //
2490  ex result = log(arg) * H(parameter,arg).hold();
2491  ex acc = 0;
2492  for (ex i=0; i<lastentry; i++) {
2493  if (parameter[i] > 0) {
2494  parameter[i]++;
2495  result -= (acc + parameter[i]-1) * H(parameter, arg).hold();
2496  parameter[i]--;
2497  acc = 0;
2498  } else if (parameter[i] < 0) {
2499  parameter[i]--;
2500  result -= (acc + abs(parameter[i]+1)) * H(parameter, arg).hold();
2501  parameter[i]++;
2502  acc = 0;
2503  } else {
2504  acc++;
2505  }
2506  }
2507 
2508  if (lastentry < parameter.nops()) {
2509  result = result / (parameter.nops()-lastentry+1);
2510  return result.map(*this);
2511  } else {
2512  return result;
2513  }
2514  }
2515  }
2516  }
2517  return e;
2518  }
2519 };
2520 
2521 
2522 // returns an expression with zeta functions corresponding to the parameter list for H
2523 ex convert_H_to_zeta(const lst& m)
2524 {
2525  symbol xtemp("xtemp");
2526  map_trafo_H_reduce_trailing_zeros filter;
2527  map_trafo_H_convert_to_zeta filter2;
2528  return filter2(filter(H(m, xtemp).hold())).subs(xtemp == 1);
2529 }
2530 
2531 
2532 // convert signs form Li to H representation
2533 lst convert_parameter_Li_to_H(const lst& m, const lst& x, ex& pf)
2534 {
2535  lst res;
2536  auto itm = m.begin();
2537  auto itx = ++x.begin();
2538  int signum = 1;
2539  pf = _ex1;
2540  res.append(*itm);
2541  itm++;
2542  while (itx != x.end()) {
2543  GINAC_ASSERT((*itx == _ex1) || (*itx == _ex_1));
2544  // XXX: 1 + 0.0*I is considered equal to 1. However the former
2545  // is not automatically converted to a real number.
2546  // Do the conversion explicitly to avoid the
2547  // "numeric::operator>(): complex inequality" exception.
2548  signum *= (*itx != _ex_1) ? 1 : -1;
2549  pf *= signum;
2550  res.append((*itm) * signum);
2551  itm++;
2552  itx++;
2553  }
2554  return res;
2555 }
2556 
2557 
2558 // multiplies an one-dimensional H with another H
2559 // [ReV] (18)
2560 ex trafo_H_mult(const ex& h1, const ex& h2)
2561 {
2562  ex res;
2563  ex hshort;
2564  lst hlong;
2565  ex h1nops = h1.op(0).nops();
2566  ex h2nops = h2.op(0).nops();
2567  if (h1nops > 1) {
2568  hshort = h2.op(0).op(0);
2569  hlong = ex_to<lst>(h1.op(0));
2570  } else {
2571  hshort = h1.op(0).op(0);
2572  if (h2nops > 1) {
2573  hlong = ex_to<lst>(h2.op(0));
2574  } else {
2575  hlong = lst{h2.op(0).op(0)};
2576  }
2577  }
2578  for (std::size_t i=0; i<=hlong.nops(); i++) {
2579  lst newparameter;
2580  std::size_t j=0;
2581  for (; j<i; j++) {
2582  newparameter.append(hlong[j]);
2583  }
2584  newparameter.append(hshort);
2585  for (; j<hlong.nops(); j++) {
2586  newparameter.append(hlong[j]);
2587  }
2588  res += H(newparameter, h1.op(1)).hold();
2589  }
2590  return res;
2591 }
2592 
2593 
2594 // applies trafo_H_mult recursively on expressions
2595 struct map_trafo_H_mult : public map_function
2596 {
2597  ex operator()(const ex& e) override
2598  {
2599  if (is_a<add>(e)) {
2600  return e.map(*this);
2601  }
2602 
2603  if (is_a<mul>(e)) {
2604 
2605  ex result = 1;
2606  ex firstH;
2607  lst Hlst;
2608  for (std::size_t pos=0; pos<e.nops(); pos++) {
2609  if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
2610  std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
2611  if (name == "H") {
2612  for (ex i=0; i<e.op(pos).op(1); i++) {
2613  Hlst.append(e.op(pos).op(0));
2614  }
2615  continue;
2616  }
2617  } else if (is_a<function>(e.op(pos))) {
2618  std::string name = ex_to<function>(e.op(pos)).get_name();
2619  if (name == "H") {
2620  if (e.op(pos).op(0).nops() > 1) {
2621  firstH = e.op(pos);
2622  } else {
2623  Hlst.append(e.op(pos));
2624  }
2625  continue;
2626  }
2627  }
2628  result *= e.op(pos);
2629  }
2630  if (firstH == 0) {
2631  if (Hlst.nops() > 0) {
2632  firstH = Hlst[Hlst.nops()-1];
2633  Hlst.remove_last();
2634  } else {
2635  return e;
2636  }
2637  }
2638 
2639  if (Hlst.nops() > 0) {
2640  ex buffer = trafo_H_mult(firstH, Hlst.op(0));
2641  result *= buffer;
2642  for (std::size_t i=1; i<Hlst.nops(); i++) {
2643  result *= Hlst.op(i);
2644  }
2645  result = result.expand();
2646  map_trafo_H_mult recursion;
2647  return recursion(result);
2648  } else {
2649  return e;
2650  }
2651 
2652  }
2653  return e;
2654  }
2655 };
2656 
2657 
2658 // do integration [ReV] (55)
2659 // put parameter 0 in front of existing parameters
2660 ex trafo_H_1tx_prepend_zero(const ex& e, const ex& arg)
2661 {
2662  ex h;
2663  std::string name;
2664  if (is_a<function>(e)) {
2665  name = ex_to<function>(e).get_name();
2666  }
2667  if (name == "H") {
2668  h = e;
2669  } else {
2670  for (std::size_t i=0; i<e.nops(); i++) {
2671  if (is_a<function>(e.op(i))) {
2672  std::string name = ex_to<function>(e.op(i)).get_name();
2673  if (name == "H") {
2674  h = e.op(i);
2675  }
2676  }
2677  }
2678  }
2679  if (h != 0) {
2680  lst newparameter = ex_to<lst>(h.op(0));
2681  newparameter.prepend(0);
2682  ex addzeta = convert_H_to_zeta(newparameter);
2683  return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
2684  } else {
2685  return e * (-H(lst{ex(0)},1/arg).hold());
2686  }
2687 }
2688 
2689 
2690 // do integration [ReV] (49)
2691 // put parameter 1 in front of existing parameters
2692 ex trafo_H_prepend_one(const ex& e, const ex& arg)
2693 {
2694  ex h;
2695  std::string name;
2696  if (is_a<function>(e)) {
2697  name = ex_to<function>(e).get_name();
2698  }
2699  if (name == "H") {
2700  h = e;
2701  } else {
2702  for (std::size_t i=0; i<e.nops(); i++) {
2703  if (is_a<function>(e.op(i))) {
2704  std::string name = ex_to<function>(e.op(i)).get_name();
2705  if (name == "H") {
2706  h = e.op(i);
2707  }
2708  }
2709  }
2710  }
2711  if (h != 0) {
2712  lst newparameter = ex_to<lst>(h.op(0));
2713  newparameter.prepend(1);
2714  return e.subs(h == H(newparameter, h.op(1)).hold());
2715  } else {
2716  return e * H(lst{ex(1)},1-arg).hold();
2717  }
2718 }
2719 
2720 
2721 // do integration [ReV] (55)
2722 // put parameter -1 in front of existing parameters
2723 ex trafo_H_1tx_prepend_minusone(const ex& e, const ex& arg)
2724 {
2725  ex h;
2726  std::string name;
2727  if (is_a<function>(e)) {
2728  name = ex_to<function>(e).get_name();
2729  }
2730  if (name == "H") {
2731  h = e;
2732  } else {
2733  for (std::size_t i=0; i<e.nops(); i++) {
2734  if (is_a<function>(e.op(i))) {
2735  std::string name = ex_to<function>(e.op(i)).get_name();
2736  if (name == "H") {
2737  h = e.op(i);
2738  }
2739  }
2740  }
2741  }
2742  if (h != 0) {
2743  lst newparameter = ex_to<lst>(h.op(0));
2744  newparameter.prepend(-1);
2745  ex addzeta = convert_H_to_zeta(newparameter);
2746  return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
2747  } else {
2748  ex addzeta = convert_H_to_zeta(lst{ex(-1)});
2749  return (e * (addzeta - H(lst{ex(-1)},1/arg).hold())).expand();
2750  }
2751 }
2752 
2753 
2754 // do integration [ReV] (55)
2755 // put parameter -1 in front of existing parameters
2756 ex trafo_H_1mxt1px_prepend_minusone(const ex& e, const ex& arg)
2757 {
2758  ex h;
2759  std::string name;
2760  if (is_a<function>(e)) {
2761  name = ex_to<function>(e).get_name();
2762  }
2763  if (name == "H") {
2764  h = e;
2765  } else {
2766  for (std::size_t i = 0; i < e.nops(); i++) {
2767  if (is_a<function>(e.op(i))) {
2768  std::string name = ex_to<function>(e.op(i)).get_name();
2769  if (name == "H") {
2770  h = e.op(i);
2771  }
2772  }
2773  }
2774  }
2775  if (h != 0) {
2776  lst newparameter = ex_to<lst>(h.op(0));
2777  newparameter.prepend(-1);
2778  return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
2779  } else {
2780  return (e * H(lst{ex(-1)},(1-arg)/(1+arg)).hold()).expand();
2781  }
2782 }
2783 
2784 
2785 // do integration [ReV] (55)
2786 // put parameter 1 in front of existing parameters
2787 ex trafo_H_1mxt1px_prepend_one(const ex& e, const ex& arg)
2788 {
2789  ex h;
2790  std::string name;
2791  if (is_a<function>(e)) {
2792  name = ex_to<function>(e).get_name();
2793  }
2794  if (name == "H") {
2795  h = e;
2796  } else {
2797  for (std::size_t i = 0; i < e.nops(); i++) {
2798  if (is_a<function>(e.op(i))) {
2799  std::string name = ex_to<function>(e.op(i)).get_name();
2800  if (name == "H") {
2801  h = e.op(i);
2802  }
2803  }
2804  }
2805  }
2806  if (h != 0) {
2807  lst newparameter = ex_to<lst>(h.op(0));
2808  newparameter.prepend(1);
2809  return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
2810  } else {
2811  return (e * H(lst{ex(1)},(1-arg)/(1+arg)).hold()).expand();
2812  }
2813 }
2814 
2815 
2816 // do x -> 1-x transformation
2817 struct map_trafo_H_1mx : public map_function
2818 {
2819  ex operator()(const ex& e) override
2820  {
2821  if (is_a<add>(e) || is_a<mul>(e)) {
2822  return e.map(*this);
2823  }
2824 
2825  if (is_a<function>(e)) {
2826  std::string name = ex_to<function>(e).get_name();
2827  if (name == "H") {
2828 
2829  lst parameter = ex_to<lst>(e.op(0));
2830  ex arg = e.op(1);
2831 
2832  // special cases if all parameters are either 0, 1 or -1
2833  bool allthesame = true;
2834  if (parameter.op(0) == 0) {
2835  for (std::size_t i = 1; i < parameter.nops(); i++) {
2836  if (parameter.op(i) != 0) {
2837  allthesame = false;
2838  break;
2839  }
2840  }
2841  if (allthesame) {
2842  lst newparameter;
2843  for (int i=parameter.nops(); i>0; i--) {
2844  newparameter.append(1);
2845  }
2846  return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
2847  }
2848  } else if (parameter.op(0) == -1) {
2849  throw std::runtime_error("map_trafo_H_1mx: cannot handle weights equal -1!");
2850  } else {
2851  for (std::size_t i = 1; i < parameter.nops(); i++) {
2852  if (parameter.op(i) != 1) {
2853  allthesame = false;
2854  break;
2855  }
2856  }
2857  if (allthesame) {
2858  lst newparameter;
2859  for (int i=parameter.nops(); i>0; i--) {
2860  newparameter.append(0);
2861  }
2862  return pow(-1, parameter.nops()) * H(newparameter, 1-arg).hold();
2863  }
2864  }
2865 
2866  lst newparameter = parameter;
2867  newparameter.remove_first();
2868 
2869  if (parameter.op(0) == 0) {
2870 
2871  // leading zero
2872  ex res = convert_H_to_zeta(parameter);
2873  //ex res = convert_from_RV(parameter, 1).subs(H(wild(1),wild(2))==zeta(wild(1)));
2874  map_trafo_H_1mx recursion;
2875  ex buffer = recursion(H(newparameter, arg).hold());
2876  if (is_a<add>(buffer)) {
2877  for (std::size_t i = 0; i < buffer.nops(); i++) {
2878  res -= trafo_H_prepend_one(buffer.op(i), arg);
2879  }
2880  } else {
2881  res -= trafo_H_prepend_one(buffer, arg);
2882  }
2883  return res;
2884 
2885  } else {
2886 
2887  // leading one
2888  map_trafo_H_1mx recursion;
2889  map_trafo_H_mult unify;
2890  ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
2891  std::size_t firstzero = 0;
2892  while (parameter.op(firstzero) == 1) {
2893  firstzero++;
2894  }
2895  for (std::size_t i = firstzero-1; i < parameter.nops()-1; i++) {
2896  lst newparameter;
2897  std::size_t j=0;
2898  for (; j<=i; j++) {
2899  newparameter.append(parameter[j+1]);
2900  }
2901  newparameter.append(1);
2902  for (; j<parameter.nops()-1; j++) {
2903  newparameter.append(parameter[j+1]);
2904  }
2905  res -= H(newparameter, arg).hold();
2906  }
2907  res = recursion(res).expand() / firstzero;
2908  return unify(res);
2909  }
2910  }
2911  }
2912  return e;
2913  }
2914 };
2915 
2916 
2917 // do x -> 1/x transformation
2918 struct map_trafo_H_1overx : public map_function
2919 {
2920  ex operator()(const ex& e) override
2921  {
2922  if (is_a<add>(e) || is_a<mul>(e)) {
2923  return e.map(*this);
2924  }
2925 
2926  if (is_a<function>(e)) {
2927  std::string name = ex_to<function>(e).get_name();
2928  if (name == "H") {
2929 
2930  lst parameter = ex_to<lst>(e.op(0));
2931  ex arg = e.op(1);
2932 
2933  // special cases if all parameters are either 0, 1 or -1
2934  bool allthesame = true;
2935  if (parameter.op(0) == 0) {
2936  for (std::size_t i = 1; i < parameter.nops(); i++) {
2937  if (parameter.op(i) != 0) {
2938  allthesame = false;
2939  break;
2940  }
2941  }
2942  if (allthesame) {
2943  return pow(-1, parameter.nops()) * H(parameter, 1/arg).hold();
2944  }
2945  } else if (parameter.op(0) == -1) {
2946  for (std::size_t i = 1; i < parameter.nops(); i++) {
2947  if (parameter.op(i) != -1) {
2948  allthesame = false;
2949  break;
2950  }
2951  }
2952  if (allthesame) {
2953  map_trafo_H_mult unify;
2954  return unify((pow(H(lst{ex(-1)},1/arg).hold() - H(lst{ex(0)},1/arg).hold(), parameter.nops())
2955  / factorial(parameter.nops())).expand());
2956  }
2957  } else {
2958  for (std::size_t i = 1; i < parameter.nops(); i++) {
2959  if (parameter.op(i) != 1) {
2960  allthesame = false;
2961  break;
2962  }
2963  }
2964  if (allthesame) {
2965  map_trafo_H_mult unify;
2966  return unify((pow(H(lst{ex(1)},1/arg).hold() + H(lst{ex(0)},1/arg).hold() + H_polesign, parameter.nops())
2967  / factorial(parameter.nops())).expand());
2968  }
2969  }
2970 
2971  lst newparameter = parameter;
2972  newparameter.remove_first();
2973 
2974  if (parameter.op(0) == 0) {
2975 
2976  // leading zero
2977  ex res = convert_H_to_zeta(parameter);
2978  map_trafo_H_1overx recursion;
2979  ex buffer = recursion(H(newparameter, arg).hold());
2980  if (is_a<add>(buffer)) {
2981  for (std::size_t i = 0; i < buffer.nops(); i++) {
2982  res += trafo_H_1tx_prepend_zero(buffer.op(i), arg);
2983  }
2984  } else {
2985  res += trafo_H_1tx_prepend_zero(buffer, arg);
2986  }
2987  return res;
2988 
2989  } else if (parameter.op(0) == -1) {
2990 
2991  // leading negative one
2992  ex res = convert_H_to_zeta(parameter);
2993  map_trafo_H_1overx recursion;
2994  ex buffer = recursion(H(newparameter, arg).hold());
2995  if (is_a<add>(buffer)) {
2996  for (std::size_t i = 0; i < buffer.nops(); i++) {
2997  res += trafo_H_1tx_prepend_zero(buffer.op(i), arg) - trafo_H_1tx_prepend_minusone(buffer.op(i), arg);
2998  }
2999  } else {
3000  res += trafo_H_1tx_prepend_zero(buffer, arg) - trafo_H_1tx_prepend_minusone(buffer, arg);
3001  }
3002  return res;
3003 
3004  } else {
3005 
3006  // leading one
3007  map_trafo_H_1overx recursion;
3008  map_trafo_H_mult unify;
3009  ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
3010  std::size_t firstzero = 0;
3011  while (parameter.op(firstzero) == 1) {
3012  firstzero++;
3013  }
3014  for (std::size_t i = firstzero-1; i < parameter.nops() - 1; i++) {
3015  lst newparameter;
3016  std::size_t j = 0;
3017  for (; j<=i; j++) {
3018  newparameter.append(parameter[j+1]);
3019  }
3020  newparameter.append(1);
3021  for (; j<parameter.nops()-1; j++) {
3022  newparameter.append(parameter[j+1]);
3023  }
3024  res -= H(newparameter, arg).hold();
3025  }
3026  res = recursion(res).expand() / firstzero;
3027  return unify(res);
3028 
3029  }
3030 
3031  }
3032  }
3033  return e;
3034  }
3035 };
3036 
3037 
3038 // do x -> (1-x)/(1+x) transformation
3039 struct map_trafo_H_1mxt1px : public map_function
3040 {
3041  ex operator()(const ex& e) override
3042  {
3043  if (is_a<add>(e) || is_a<mul>(e)) {
3044  return e.map(*this);
3045  }
3046 
3047  if (is_a<function>(e)) {
3048  std::string name = ex_to<function>(e).get_name();
3049  if (name == "H") {
3050 
3051  lst parameter = ex_to<lst>(e.op(0));
3052  ex arg = e.op(1);
3053 
3054  // special cases if all parameters are either 0, 1 or -1
3055  bool allthesame = true;
3056  if (parameter.op(0) == 0) {
3057  for (std::size_t i = 1; i < parameter.nops(); i++) {
3058  if (parameter.op(i) != 0) {
3059  allthesame = false;
3060  break;
3061  }
3062  }
3063  if (allthesame) {
3064  map_trafo_H_mult unify;
3065  return unify((pow(-H(lst{ex(1)},(1-arg)/(1+arg)).hold() - H(lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.nops())
3066  / factorial(parameter.nops())).expand());
3067  }
3068  } else if (parameter.op(0) == -1) {
3069  for (std::size_t i = 1; i < parameter.nops(); i++) {
3070  if (parameter.op(i) != -1) {
3071  allthesame = false;
3072  break;
3073  }
3074  }
3075  if (allthesame) {
3076  map_trafo_H_mult unify;
3077  return unify((pow(log(2) - H(lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.nops())
3078  / factorial(parameter.nops())).expand());
3079  }
3080  } else {
3081  for (std::size_t i = 1; i < parameter.nops(); i++) {
3082  if (parameter.op(i) != 1) {
3083  allthesame = false;
3084  break;
3085  }
3086  }
3087  if (allthesame) {
3088  map_trafo_H_mult unify;
3089  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())
3090  / factorial(parameter.nops())).expand());
3091  }
3092  }
3093 
3094  lst newparameter = parameter;
3095  newparameter.remove_first();
3096 
3097  if (parameter.op(0) == 0) {
3098 
3099  // leading zero
3100  ex res = convert_H_to_zeta(parameter);
3101  map_trafo_H_1mxt1px recursion;
3102  ex buffer = recursion(H(newparameter, arg).hold());
3103  if (is_a<add>(buffer)) {
3104  for (std::size_t i = 0; i < buffer.nops(); i++) {
3105  res -= trafo_H_1mxt1px_prepend_one(buffer.op(i), arg) + trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
3106  }
3107  } else {
3108  res -= trafo_H_1mxt1px_prepend_one(buffer, arg) + trafo_H_1mxt1px_prepend_minusone(buffer, arg);
3109  }
3110  return res;
3111 
3112  } else if (parameter.op(0) == -1) {
3113 
3114  // leading negative one
3115  ex res = convert_H_to_zeta(parameter);
3116  map_trafo_H_1mxt1px recursion;
3117  ex buffer = recursion(H(newparameter, arg).hold());
3118  if (is_a<add>(buffer)) {
3119  for (std::size_t i = 0; i < buffer.nops(); i++) {
3120  res -= trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
3121  }
3122  } else {
3123  res -= trafo_H_1mxt1px_prepend_minusone(buffer, arg);
3124  }
3125  return res;
3126 
3127  } else {
3128 
3129  // leading one
3130  map_trafo_H_1mxt1px recursion;
3131  map_trafo_H_mult unify;
3132  ex res = H(lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
3133  std::size_t firstzero = 0;
3134  while (parameter.op(firstzero) == 1) {
3135  firstzero++;
3136  }
3137  for (std::size_t i = firstzero - 1; i < parameter.nops() - 1; i++) {
3138  lst newparameter;
3139  std::size_t j=0;
3140  for (; j<=i; j++) {
3141  newparameter.append(parameter[j+1]);
3142  }
3143  newparameter.append(1);
3144  for (; j<parameter.nops()-1; j++) {
3145  newparameter.append(parameter[j+1]);
3146  }
3147  res -= H(newparameter, arg).hold();
3148  }
3149  res = recursion(res).expand() / firstzero;
3150  return unify(res);
3151 
3152  }
3153 
3154  }
3155  }
3156  return e;
3157  }
3158 };
3159 
3160 
3161 // do the actual summation.
3162 cln::cl_N H_do_sum(const std::vector<int>& m, const cln::cl_N& x)
3163 {
3164  const int j = m.size();
3165 
3166  std::vector<cln::cl_N> t(j);
3167 
3168  cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
3169  cln::cl_N factor = cln::expt(x, j) * one;
3170  cln::cl_N t0buf;
3171  int q = 0;
3172  do {
3173  t0buf = t[0];
3174  q++;
3175  t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),m[j-1]);
3176  for (int k=j-2; k>=1; k--) {
3177  t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), m[k]);
3178  }
3179  t[0] = t[0] + t[1] * factor / cln::expt(cln::cl_I(q+j-1), m[0]);
3180  factor = factor * x;
3181  } while (t[0] != t0buf);
3182 
3183  return t[0];
3184 }
3185 
3186 
3187 } // end of anonymous namespace
3188 
3189 
3191 //
3192 // Harmonic polylogarithm H(m,x)
3193 //
3194 // GiNaC function
3195 //
3197 
3198 
3199 static ex H_evalf(const ex& x1, const ex& x2)
3200 {
3201  if (is_a<lst>(x1)) {
3202 
3203  cln::cl_N x;
3204  if (is_a<numeric>(x2)) {
3205  x = ex_to<numeric>(x2).to_cl_N();
3206  } else {
3207  ex x2_val = x2.evalf();
3208  if (is_a<numeric>(x2_val)) {
3209  x = ex_to<numeric>(x2_val).to_cl_N();
3210  }
3211  }
3212 
3213  for (std::size_t i = 0; i < x1.nops(); i++) {
3214  if (!x1.op(i).info(info_flags::integer)) {
3215  return H(x1, x2).hold();
3216  }
3217  }
3218  if (x1.nops() < 1) {
3219  return H(x1, x2).hold();
3220  }
3221 
3222  const lst& morg = ex_to<lst>(x1);
3223  // remove trailing zeros ...
3224  if (*(--morg.end()) == 0) {
3225  symbol xtemp("xtemp");
3226  map_trafo_H_reduce_trailing_zeros filter;
3227  return filter(H(x1, xtemp).hold()).subs(xtemp==x2).evalf();
3228  }
3229  // ... and expand parameter notation
3230  lst m;
3231  for (const auto & it : morg) {
3232  if (it > 1) {
3233  for (ex count=it-1; count > 0; count--) {
3234  m.append(0);
3235  }
3236  m.append(1);
3237  } else if (it <= -1) {
3238  for (ex count=it+1; count < 0; count++) {
3239  m.append(0);
3240  }
3241  m.append(-1);
3242  } else {
3243  m.append(it);
3244  }
3245  }
3246 
3247  // do summation
3248  if (cln::abs(x) < 0.95) {
3249  lst m_lst;
3250  lst s_lst;
3251  ex pf;
3252  if (convert_parameter_H_to_Li(m, m_lst, s_lst, pf)) {
3253  // negative parameters -> s_lst is filled
3254  std::vector<int> m_int;
3255  std::vector<cln::cl_N> x_cln;
3256  for (auto it_int = m_lst.begin(), it_cln = s_lst.begin();
3257  it_int != m_lst.end(); it_int++, it_cln++) {
3258  m_int.push_back(ex_to<numeric>(*it_int).to_int());
3259  x_cln.push_back(ex_to<numeric>(*it_cln).to_cl_N());
3260  }
3261  x_cln.front() = x_cln.front() * x;
3262  return pf * numeric(multipleLi_do_sum(m_int, x_cln));
3263  } else {
3264  // only positive parameters
3265  //TODO
3266  if (m_lst.nops() == 1) {
3267  return Li(m_lst.op(0), x2).evalf();
3268  }
3269  std::vector<int> m_int;
3270  for (const auto & it : m_lst) {
3271  m_int.push_back(ex_to<numeric>(it).to_int());
3272  }
3273  return numeric(H_do_sum(m_int, x));
3274  }
3275  }
3276 
3277  symbol xtemp("xtemp");
3278  ex res = 1;
3279 
3280  // ensure that the realpart of the argument is positive
3281  if (cln::realpart(x) < 0) {
3282  x = -x;
3283  for (std::size_t i = 0; i < m.nops(); i++) {
3284  if (m.op(i) != 0) {
3285  m.let_op(i) = -m.op(i);
3286  res *= -1;
3287  }
3288  }
3289  }
3290 
3291  // x -> 1/x
3292  if (cln::abs(x) >= 2.0) {
3293  map_trafo_H_1overx trafo;
3294  res *= trafo(H(m, xtemp).hold());
3295  if (cln::imagpart(x) <= 0) {
3296  res = res.subs(H_polesign == -I*Pi);
3297  } else {
3298  res = res.subs(H_polesign == I*Pi);
3299  }
3300  return res.subs(xtemp == numeric(x)).evalf();
3301  }
3302 
3303  // check for letters (-1)
3304  bool has_minus_one = false;
3305  for (const auto & it : m) {
3306  if (it == -1)
3307  has_minus_one = true;
3308  }
3309 
3310  // check transformations for 0.95 <= |x| < 2.0
3311 
3312  // |(1-x)/(1+x)| < 0.9 -> circular area with center=9.53+0i and radius=9.47
3313  if (cln::abs(x-9.53) <= 9.47) {
3314  // x -> (1-x)/(1+x)
3315  map_trafo_H_1mxt1px trafo;
3316  res *= trafo(H(m, xtemp).hold());
3317  } else {
3318  // x -> 1-x
3319  if (has_minus_one) {
3320  map_trafo_H_convert_to_Li filter;
3321  // 09.06.2021: bug fix: don't forget a possible minus sign from the case realpart(x) < 0
3322  res *= filter(H(m, numeric(x)).hold()).evalf();
3323  return res;
3324  }
3325  map_trafo_H_1mx trafo;
3326  res *= trafo(H(m, xtemp).hold());
3327  }
3328 
3329  return res.subs(xtemp == numeric(x)).evalf();
3330  }
3331 
3332  return H(x1,x2).hold();
3333 }
3334 
3335 
3336 static ex H_eval(const ex& m_, const ex& x)
3337 {
3338  lst m;
3339  if (is_a<lst>(m_)) {
3340  m = ex_to<lst>(m_);
3341  } else {
3342  m = lst{m_};
3343  }
3344  if (m.nops() == 0) {
3345  return _ex1;
3346  }
3347  ex pos1;
3348  ex pos2;
3349  ex n;
3350  ex p;
3351  int step = 0;
3352  if (*m.begin() > _ex1) {
3353  step++;
3354  pos1 = _ex0;
3355  pos2 = _ex1;
3356  n = *m.begin()-1;
3357  p = _ex1;
3358  } else if (*m.begin() < _ex_1) {
3359  step++;
3360  pos1 = _ex0;
3361  pos2 = _ex_1;
3362  n = -*m.begin()-1;
3363  p = _ex1;
3364  } else if (*m.begin() == _ex0) {
3365  pos1 = _ex0;
3366  n = _ex1;
3367  } else {
3368  pos1 = *m.begin();
3369  p = _ex1;
3370  }
3371  for (auto it = ++m.begin(); it != m.end(); it++) {
3372  if (it->info(info_flags::integer)) {
3373  if (step == 0) {
3374  if (*it > _ex1) {
3375  if (pos1 == _ex0) {
3376  step = 1;
3377  pos2 = _ex1;
3378  n += *it-1;
3379  p = _ex1;
3380  } else {
3381  step = 2;
3382  }
3383  } else if (*it < _ex_1) {
3384  if (pos1 == _ex0) {
3385  step = 1;
3386  pos2 = _ex_1;
3387  n += -*it-1;
3388  p = _ex1;
3389  } else {
3390  step = 2;
3391  }
3392  } else {
3393  if (*it != pos1) {
3394  step = 1;
3395  pos2 = *it;
3396  }
3397  if (*it == _ex0) {
3398  n++;
3399  } else {
3400  p++;
3401  }
3402  }
3403  } else if (step == 1) {
3404  if (*it != pos2) {
3405  step = 2;
3406  } else {
3407  if (*it == _ex0) {
3408  n++;
3409  } else {
3410  p++;
3411  }
3412  }
3413  }
3414  } else {
3415  // if some m_i is not an integer
3416  return H(m_, x).hold();
3417  }
3418  }
3419  if ((x == _ex1) && (*(--m.end()) != _ex0)) {
3420  return convert_H_to_zeta(m);
3421  }
3422  if (step == 0) {
3423  if (pos1 == _ex0) {
3424  // all zero
3425  if (x == _ex0) {
3426  return H(m_, x).hold();
3427  }
3428  return pow(log(x), m.nops()) / factorial(m.nops());
3429  } else {
3430  // all (minus) one
3431  return pow(-pos1*log(1-pos1*x), m.nops()) / factorial(m.nops());
3432  }
3433  } else if ((step == 1) && (pos1 == _ex0)){
3434  // convertible to S
3435  if (pos2 == _ex1) {
3436  return S(n, p, x);
3437  } else {
3438  return pow(-1, p) * S(n, p, -x);
3439  }
3440  }
3441  if (x == _ex0) {
3442  return _ex0;
3443  }
3445  return H(m_, x).evalf();
3446  }
3447  return H(m_, x).hold();
3448 }
3449 
3450 
3451 static ex H_series(const ex& m, const ex& x, const relational& rel, int order, unsigned options)
3452 {
3453  epvector seq { expair(H(m, x), 0) };
3454  return pseries(rel, std::move(seq));
3455 }
3456 
3457 
3458 static ex H_deriv(const ex& m_, const ex& x, unsigned deriv_param)
3459 {
3460  GINAC_ASSERT(deriv_param < 2);
3461  if (deriv_param == 0) {
3462  return _ex0;
3463  }
3464  lst m;
3465  if (is_a<lst>(m_)) {
3466  m = ex_to<lst>(m_);
3467  } else {
3468  m = lst{m_};
3469  }
3470  ex mb = *m.begin();
3471  if (mb > _ex1) {
3472  m[0]--;
3473  return H(m, x) / x;
3474  }
3475  if (mb < _ex_1) {
3476  m[0]++;
3477  return H(m, x) / x;
3478  }
3479  m.remove_first();
3480  if (mb == _ex1) {
3481  return 1/(1-x) * H(m, x);
3482  } else if (mb == _ex_1) {
3483  return 1/(1+x) * H(m, x);
3484  } else {
3485  return H(m, x) / x;
3486  }
3487 }
3488 
3489 
3490 static void H_print_latex(const ex& m_, const ex& x, const print_context& c)
3491 {
3492  lst m;
3493  if (is_a<lst>(m_)) {
3494  m = ex_to<lst>(m_);
3495  } else {
3496  m = lst{m_};
3497  }
3498  c.s << "\\mathrm{H}_{";
3499  auto itm = m.begin();
3500  (*itm).print(c);
3501  itm++;
3502  for (; itm != m.end(); itm++) {
3503  c.s << ",";
3504  (*itm).print(c);
3505  }
3506  c.s << "}(";
3507  x.print(c);
3508  c.s << ")";
3509 }
3510 
3511 
3513  evalf_func(H_evalf).
3514  eval_func(H_eval).
3515  series_func(H_series).
3516  derivative_func(H_deriv).
3517  print_func<print_latex>(H_print_latex).
3518  do_not_evalf_params());
3519 
3520 
3521 // takes a parameter list for H and returns an expression with corresponding multiple polylogarithms
3522 ex convert_H_to_Li(const ex& m, const ex& x)
3523 {
3524  map_trafo_H_reduce_trailing_zeros filter;
3525  map_trafo_H_convert_to_Li filter2;
3526  if (is_a<lst>(m)) {
3527  return filter2(filter(H(m, x).hold()));
3528  } else {
3529  return filter2(filter(H(lst{m}, x).hold()));
3530  }
3531 }
3532 
3533 
3535 //
3536 // Multiple zeta values zeta(x) and zeta(x,s)
3537 //
3538 // helper functions
3539 //
3541 
3542 
3543 // anonymous namespace for helper functions
3544 namespace {
3545 
3546 
3547 // parameters and data for [Cra] algorithm
3548 const cln::cl_N lambda = cln::cl_N("319/320");
3549 
3550 void halfcyclic_convolute(const std::vector<cln::cl_N>& a, const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>& c)
3551 {
3552  const int size = a.size();
3553  for (int n=0; n<size; n++) {
3554  c[n] = 0;
3555  for (int m=0; m<=n; m++) {
3556  c[n] = c[n] + a[m]*b[n-m];
3557  }
3558  }
3559 }
3560 
3561 
3562 // [Cra] section 4
3563 static void initcX(std::vector<cln::cl_N>& crX,
3564  const std::vector<int>& s,
3565  const int L2)
3566 {
3567  std::vector<cln::cl_N> crB(L2 + 1);
3568  for (int i=0; i<=L2; i++)
3569  crB[i] = bernoulli(i).to_cl_N() / cln::factorial(i);
3570 
3571  int Sm = 0;
3572  int Smp1 = 0;
3573  std::vector<std::vector<cln::cl_N>> crG(s.size() - 1, std::vector<cln::cl_N>(L2 + 1));
3574  for (int m=0; m < (int)s.size() - 1; m++) {
3575  Sm += s[m];
3576  Smp1 = Sm + s[m+1];
3577  for (int i = 0; i <= L2; i++)
3578  crG[m][i] = cln::factorial(i + Sm - m - 2) / cln::factorial(i + Smp1 - m - 2);
3579  }
3580 
3581  crX = crB;
3582 
3583  for (std::size_t m = 0; m < s.size() - 1; m++) {
3584  std::vector<cln::cl_N> Xbuf(L2 + 1);
3585  for (int i = 0; i <= L2; i++)
3586  Xbuf[i] = crX[i] * crG[m][i];
3587 
3588  halfcyclic_convolute(Xbuf, crB, crX);
3589  }
3590 }
3591 
3592 
3593 // [Cra] section 4
3594 static cln::cl_N crandall_Y_loop(const cln::cl_N& Sqk,
3595  const std::vector<cln::cl_N>& crX)
3596 {
3597  cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
3598  cln::cl_N factor = cln::expt(lambda, Sqk);
3599  cln::cl_N res = factor / Sqk * crX[0] * one;
3600  cln::cl_N resbuf;
3601  int N = 0;
3602  do {
3603  resbuf = res;
3604  factor = factor * lambda;
3605  N++;
3606  res = res + crX[N] * factor / (N+Sqk);
3607  } while (((res != resbuf) || cln::zerop(crX[N])) && (N+1 < crX.size()));
3608  return res;
3609 }
3610 
3611 
3612 // [Cra] section 4
3613 static void calc_f(std::vector<std::vector<cln::cl_N>>& f_kj,
3614  const int maxr, const int L1)
3615 {
3616  cln::cl_N t0, t1, t2, t3, t4;
3617  int i, j, k;
3618  auto it = f_kj.begin();
3619  cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
3620 
3621  t0 = cln::exp(-lambda);
3622  t2 = 1;
3623  for (k=1; k<=L1; k++) {
3624  t1 = k * lambda;
3625  t2 = t0 * t2;
3626  for (j=1; j<=maxr; j++) {
3627  t3 = 1;
3628  t4 = 1;
3629  for (i=2; i<=j; i++) {
3630  t4 = t4 * (j-i+1);
3631  t3 = t1 * t3 + t4;
3632  }
3633  (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(k),-j) * one);
3634  }
3635  it++;
3636  }
3637 }
3638 
3639 
3640 // [Cra] (3.1)
3641 static cln::cl_N crandall_Z(const std::vector<int>& s,
3642  const std::vector<std::vector<cln::cl_N>>& f_kj)
3643 {
3644  const int j = s.size();
3645 
3646  if (j == 1) {
3647  cln::cl_N t0;
3648  cln::cl_N t0buf;
3649  int q = 0;
3650  do {
3651  t0buf = t0;
3652  q++;
3653  t0 = t0 + f_kj[q+j-2][s[0]-1];
3654  } while ((t0 != t0buf) && (q+j-1 < f_kj.size()));
3655 
3656  return t0 / cln::factorial(s[0]-1);
3657  }
3658 
3659  std::vector<cln::cl_N> t(j);
3660 
3661  cln::cl_N t0buf;
3662  int q = 0;
3663  do {
3664  t0buf = t[0];
3665  q++;
3666  t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
3667  for (int k=j-2; k>=1; k--) {
3668  t[k] = t[k] + t[k+1] / cln::expt(cln::cl_I(q+j-1-k), s[k]);
3669  }
3670  t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
3671  } while ((t[0] != t0buf) && (q+j-1 < f_kj.size()));
3672 
3673  return t[0] / cln::factorial(s[0]-1);
3674 }
3675 
3676 
3677 // [Cra] (2.4)
3678 cln::cl_N zeta_do_sum_Crandall(const std::vector<int>& s)
3679 {
3680  std::vector<int> r = s;
3681  const int j = r.size();
3682 
3683  std::size_t L1;
3684 
3685  // decide on maximal size of f_kj for crandall_Z
3686  if (Digits < 50) {
3687  L1 = 150;
3688  } else {
3689  L1 = Digits * 3 + j*2;
3690  }
3691 
3692  std::size_t L2;
3693  // decide on maximal size of crX for crandall_Y
3694  if (Digits < 38) {
3695  L2 = 63;
3696  } else if (Digits < 86) {
3697  L2 = 127;
3698  } else if (Digits < 192) {
3699  L2 = 255;
3700  } else if (Digits < 394) {
3701  L2 = 511;
3702  } else if (Digits < 808) {
3703  L2 = 1023;
3704  } else if (Digits < 1636) {
3705  L2 = 2047;
3706  } else {
3707  // [Cra] section 6, log10(lambda/2/Pi) approx -0.79 for lambda=319/320, add some extra digits
3708  L2 = std::pow(2, ceil( std::log2((long(Digits))/0.79 + 40 )) ) - 1;
3709  }
3710 
3711  cln::cl_N res;
3712 
3713  int maxr = 0;
3714  int S = 0;
3715  for (int i=0; i<j; i++) {
3716  S += r[i];
3717  if (r[i] > maxr) {
3718  maxr = r[i];
3719  }
3720  }
3721 
3722  std::vector<std::vector<cln::cl_N>> f_kj(L1);
3723  calc_f(f_kj, maxr, L1);
3724 
3725  const cln::cl_N r0factorial = cln::factorial(r[0]-1);
3726 
3727  std::vector<int> rz;
3728  int skp1buf;
3729  int Srun = S;
3730  for (int k=r.size()-1; k>0; k--) {
3731 
3732  rz.insert(rz.begin(), r.back());
3733  skp1buf = rz.front();
3734  Srun -= skp1buf;
3735  r.pop_back();
3736 
3737  std::vector<cln::cl_N> crX;
3738  initcX(crX, r, L2);
3739 
3740  for (int q=0; q<skp1buf; q++) {
3741 
3742  cln::cl_N pp1 = crandall_Y_loop(Srun+q-k, crX);
3743  cln::cl_N pp2 = crandall_Z(rz, f_kj);
3744 
3745  rz.front()--;
3746 
3747  if (q & 1) {
3748  res = res - pp1 * pp2 / cln::factorial(q);
3749  } else {
3750  res = res + pp1 * pp2 / cln::factorial(q);
3751  }
3752  }
3753  rz.front() = skp1buf;
3754  }
3755  rz.insert(rz.begin(), r.back());
3756 
3757  std::vector<cln::cl_N> crX;
3758  initcX(crX, rz, L2);
3759 
3760  res = (res + crandall_Y_loop(S-j, crX)) / r0factorial
3761  + crandall_Z(rz, f_kj);
3762 
3763  return res;
3764 }
3765 
3766 
3767 cln::cl_N zeta_do_sum_simple(const std::vector<int>& r)
3768 {
3769  const int j = r.size();
3770 
3771  // buffer for subsums
3772  std::vector<cln::cl_N> t(j);
3773  cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
3774 
3775  cln::cl_N t0buf;
3776  int q = 0;
3777  do {
3778  t0buf = t[0];
3779  q++;
3780  t[j-1] = t[j-1] + one / cln::expt(cln::cl_I(q),r[j-1]);
3781  for (int k=j-2; k>=0; k--) {
3782  t[k] = t[k] + one * t[k+1] / cln::expt(cln::cl_I(q+j-1-k), r[k]);
3783  }
3784  } while (t[0] != t0buf);
3785 
3786  return t[0];
3787 }
3788 
3789 
3790 // does Hoelder convolution. see [BBB] (7.0)
3791 cln::cl_N zeta_do_Hoelder_convolution(const std::vector<int>& m_, const std::vector<int>& s_)
3792 {
3793  // prepare parameters
3794  // holds Li arguments in [BBB] notation
3795  std::vector<int> s = s_;
3796  std::vector<int> m_p = m_;
3797  std::vector<int> m_q;
3798  // holds Li arguments in nested sums notation
3799  std::vector<cln::cl_N> s_p(s.size(), cln::cl_N(1));
3800  s_p[0] = s_p[0] * cln::cl_N("1/2");
3801  // convert notations
3802  int sig = 1;
3803  for (std::size_t i = 0; i < s_.size(); i++) {
3804  if (s_[i] < 0) {
3805  sig = -sig;
3806  s_p[i] = -s_p[i];
3807  }
3808  s[i] = sig * std::abs(s[i]);
3809  }
3810  std::vector<cln::cl_N> s_q;
3811  cln::cl_N signum = 1;
3812 
3813  // first term
3814  cln::cl_N res = multipleLi_do_sum(m_p, s_p);
3815 
3816  // middle terms
3817  do {
3818 
3819  // change parameters
3820  if (s.front() > 0) {
3821  if (m_p.front() == 1) {
3822  m_p.erase(m_p.begin());
3823  s_p.erase(s_p.begin());
3824  if (s_p.size() > 0) {
3825  s_p.front() = s_p.front() * cln::cl_N("1/2");
3826  }
3827  s.erase(s.begin());
3828  m_q.front()++;
3829  } else {
3830  m_p.front()--;
3831  m_q.insert(m_q.begin(), 1);
3832  if (s_q.size() > 0) {
3833  s_q.front() = s_q.front() * 2;
3834  }
3835  s_q.insert(s_q.begin(), cln::cl_N("1/2"));
3836  }
3837  } else {
3838  if (m_p.front() == 1) {
3839  m_p.erase(m_p.begin());
3840  cln::cl_N spbuf = s_p.front();
3841  s_p.erase(s_p.begin());
3842  if (s_p.size() > 0) {
3843  s_p.front() = s_p.front() * spbuf;
3844  }
3845  s.erase(s.begin());
3846  m_q.insert(m_q.begin(), 1);
3847  if (s_q.size() > 0) {
3848  s_q.front() = s_q.front() * 4;
3849  }
3850  s_q.insert(s_q.begin(), cln::cl_N("1/4"));
3851  signum = -signum;
3852  } else {
3853  m_p.front()--;
3854  m_q.insert(m_q.begin(), 1);
3855  if (s_q.size() > 0) {
3856  s_q.front() = s_q.front() * 2;
3857  }
3858  s_q.insert(s_q.begin(), cln::cl_N("1/2"));
3859  }
3860  }
3861 
3862  // exiting the loop
3863  if (m_p.size() == 0) break;
3864 
3865  res = res + signum * multipleLi_do_sum(m_p, s_p) * multipleLi_do_sum(m_q, s_q);
3866 
3867  } while (true);
3868 
3869  // last term
3870  res = res + signum * multipleLi_do_sum(m_q, s_q);
3871 
3872  return res;
3873 }
3874 
3875 
3876 } // end of anonymous namespace
3877 
3878 
3880 //
3881 // Multiple zeta values zeta(x)
3882 //
3883 // GiNaC function
3884 //
3886 
3887 
3888 static ex zeta1_evalf(const ex& x)
3889 {
3890  if (is_exactly_a<lst>(x) && (x.nops()>1)) {
3891 
3892  // multiple zeta value
3893  const int count = x.nops();
3894  const lst& xlst = ex_to<lst>(x);
3895  std::vector<int> r(count);
3896  std::vector<int> si(count);
3897 
3898  // check parameters and convert them
3899  auto it1 = xlst.begin();
3900  auto it2 = r.begin();
3901  auto it_swrite = si.begin();
3902  do {
3903  if (!(*it1).info(info_flags::posint)) {
3904  return zeta(x).hold();
3905  }
3906  *it2 = ex_to<numeric>(*it1).to_int();
3907  *it_swrite = 1;
3908  it1++;
3909  it2++;
3910  it_swrite++;
3911  } while (it2 != r.end());
3912 
3913  // check for divergence
3914  if (r[0] == 1) {
3915  return zeta(x).hold();
3916  }
3917 
3918  // use Hoelder convolution if Digits is large
3919  if (Digits>50)
3920  return numeric(zeta_do_Hoelder_convolution(r, si));
3921 
3922  // decide on summation algorithm
3923  // this is still a bit clumsy
3924  int limit = (Digits>17) ? 10 : 6;
3925  if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) {
3926  return numeric(zeta_do_sum_Crandall(r));
3927  } else {
3928  return numeric(zeta_do_sum_simple(r));
3929  }
3930  }
3931 
3932  // single zeta value
3933  if (is_exactly_a<numeric>(x) && (x != 1)) {
3934  try {
3935  return zeta(ex_to<numeric>(x));
3936  } catch (const dunno &e) { }
3937  }
3938 
3939  return zeta(x).hold();
3940 }
3941 
3942 
3943 static ex zeta1_eval(const ex& m)
3944 {
3945  if (is_exactly_a<lst>(m)) {
3946  if (m.nops() == 1) {
3947  return zeta(m.op(0));
3948  }
3949  return zeta(m).hold();
3950  }
3951 
3952  if (m.info(info_flags::numeric)) {
3953  const numeric& y = ex_to<numeric>(m);
3954  // trap integer arguments:
3955  if (y.is_integer()) {
3956  if (y.is_zero()) {
3957  return _ex_1_2;
3958  }
3959  if (y.is_equal(*_num1_p)) {
3960  return zeta(m).hold();
3961  }
3962  if (y.info(info_flags::posint)) {
3963  if (y.info(info_flags::odd)) {
3964  return zeta(m).hold();
3965  } else {
3966  return abs(bernoulli(y)) * pow(Pi, y) * pow(*_num2_p, y-(*_num1_p)) / factorial(y);
3967  }
3968  } else {
3969  if (y.info(info_flags::odd)) {
3970  return -bernoulli((*_num1_p)-y) / ((*_num1_p)-y);
3971  } else {
3972  return _ex0;
3973  }
3974  }
3975  }
3976  // zeta(float)
3978  return zeta1_evalf(m);
3979  }
3980  }
3981  return zeta(m).hold();
3982 }
3983 
3984 
3985 static ex zeta1_deriv(const ex& m, unsigned deriv_param)
3986 {
3987  GINAC_ASSERT(deriv_param==0);
3988 
3989  if (is_exactly_a<lst>(m)) {
3990  return _ex0;
3991  } else {
3992  return zetaderiv(_ex1, m);
3993  }
3994 }
3995 
3996 
3997 static void zeta1_print_latex(const ex& m_, const print_context& c)
3998 {
3999  c.s << "\\zeta(";
4000  if (is_a<lst>(m_)) {
4001  const lst& m = ex_to<lst>(m_);
4002  auto it = m.begin();
4003  (*it).print(c);
4004  it++;
4005  for (; it != m.end(); it++) {
4006  c.s << ",";
4007  (*it).print(c);
4008  }
4009  } else {
4010  m_.print(c);
4011  }
4012  c.s << ")";
4013 }
4014 
4015 
4016 unsigned zeta1_SERIAL::serial = function::register_new(function_options("zeta", 1).
4017  evalf_func(zeta1_evalf).
4018  eval_func(zeta1_eval).
4019  derivative_func(zeta1_deriv).
4020  print_func<print_latex>(zeta1_print_latex).
4021  do_not_evalf_params().
4022  overloaded(2));
4023 
4024 
4026 //
4027 // Alternating Euler sum zeta(x,s)
4028 //
4029 // GiNaC function
4030 //
4032 
4033 
4034 static ex zeta2_evalf(const ex& x, const ex& s)
4035 {
4036  if (is_exactly_a<lst>(x)) {
4037 
4038  // alternating Euler sum
4039  const int count = x.nops();
4040  const lst& xlst = ex_to<lst>(x);
4041  const lst& slst = ex_to<lst>(s);
4042  std::vector<int> xi(count);
4043  std::vector<int> si(count);
4044 
4045  // check parameters and convert them
4046  auto it_xread = xlst.begin();
4047  auto it_sread = slst.begin();
4048  auto it_xwrite = xi.begin();
4049  auto it_swrite = si.begin();
4050  do {
4051  if (!(*it_xread).info(info_flags::posint)) {
4052  return zeta(x, s).hold();
4053  }
4054  *it_xwrite = ex_to<numeric>(*it_xread).to_int();
4055  if (*it_sread > 0) {
4056  *it_swrite = 1;
4057  } else {
4058  *it_swrite = -1;
4059  }
4060  it_xread++;
4061  it_sread++;
4062  it_xwrite++;
4063  it_swrite++;
4064  } while (it_xwrite != xi.end());
4065 
4066  // check for divergence
4067  if ((xi[0] == 1) && (si[0] == 1)) {
4068  return zeta(x, s).hold();
4069  }
4070 
4071  // use Hoelder convolution
4072  return numeric(zeta_do_Hoelder_convolution(xi, si));
4073  }
4074 
4075  // x and s are not lists: convert to lists
4076  return zeta(lst{x}, lst{s}).evalf();
4077 }
4078 
4079 
4080 static ex zeta2_eval(const ex& m, const ex& s_)
4081 {
4082  if (is_exactly_a<lst>(s_)) {
4083  const lst& s = ex_to<lst>(s_);
4084  for (const auto & it : s) {
4085  if (it.info(info_flags::positive)) {
4086  continue;
4087  }
4088  return zeta(m, s_).hold();
4089  }
4090  return zeta(m);
4091  } else if (s_.info(info_flags::positive)) {
4092  return zeta(m);
4093  }
4094 
4095  return zeta(m, s_).hold();
4096 }
4097 
4098 
4099 static ex zeta2_deriv(const ex& m, const ex& s, unsigned deriv_param)
4100 {
4101  GINAC_ASSERT(deriv_param==0);
4102 
4103  if (is_exactly_a<lst>(m)) {
4104  return _ex0;
4105  } else {
4106  if ((is_exactly_a<lst>(s) && s.op(0).info(info_flags::positive)) || s.info(info_flags::positive)) {
4107  return zetaderiv(_ex1, m);
4108  }
4109  return _ex0;
4110  }
4111 }
4112 
4113 
4114 static void zeta2_print_latex(const ex& m_, const ex& s_, const print_context& c)
4115 {
4116  lst m;
4117  if (is_a<lst>(m_)) {
4118  m = ex_to<lst>(m_);
4119  } else {
4120  m = lst{m_};
4121  }
4122  lst s;
4123  if (is_a<lst>(s_)) {
4124  s = ex_to<lst>(s_);
4125  } else {
4126  s = lst{s_};
4127  }
4128  c.s << "\\zeta(";
4129  auto itm = m.begin();
4130  auto its = s.begin();
4131  if (*its < 0) {
4132  c.s << "\\overline{";
4133  (*itm).print(c);
4134  c.s << "}";
4135  } else {
4136  (*itm).print(c);
4137  }
4138  its++;
4139  itm++;
4140  for (; itm != m.end(); itm++, its++) {
4141  c.s << ",";
4142  if (*its < 0) {
4143  c.s << "\\overline{";
4144  (*itm).print(c);
4145  c.s << "}";
4146  } else {
4147  (*itm).print(c);
4148  }
4149  }
4150  c.s << ")";
4151 }
4152 
4153 
4154 unsigned zeta2_SERIAL::serial = function::register_new(function_options("zeta", 2).
4155  evalf_func(zeta2_evalf).
4156  eval_func(zeta2_eval).
4157  derivative_func(zeta2_deriv).
4158  print_func<print_latex>(zeta2_print_latex).
4159  do_not_evalf_params().
4160  overloaded(2));
4161 
4162 
4163 } // namespace GiNaC
4164 
Interface to GiNaC's sums of expressions.
#define GINAC_ASSERT(X)
Assertion macro for checking invariances.
Definition: assertion.h:33
static unsigned serial
Definition: inifcns.h:128
static unsigned serial
Definition: inifcns.h:134
virtual size_t nops() const
Number of operands/members.
Definition: basic.cpp:229
const basic & hold() const
Stop further evaluation.
Definition: basic.cpp:887
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
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
const_iterator end() const
Definition: container.h:240
const_iterator begin() const
Definition: container.h:239
size_t nops() const override
Number of operands/members.
Definition: container.h:118
ex op(size_t i) const override
Return operand/member at position i.
Definition: container.h:295
container & append(const ex &b)
Add element at back.
Definition: container.h:391
Exception class thrown by classes which provide their own series expansion to signal that ordinary Ta...
Definition: function.h:668
Exception class thrown by functions to signal unimplemented functionality so the expression may just ...
Definition: utils.h:37
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
ex map(map_function &f) const
Definition: ex.h:162
const_iterator begin() const noexcept
Definition: ex.h:647
bool is_equal(const ex &other) const
Definition: ex.h:345
ex & let_op(size_t i)
Return modifiable operand/member at position i.
Definition: ex.cpp:206
ex evalf() const
Definition: ex.h:121
const_iterator end() const noexcept
Definition: ex.h:652
size_t nops() const
Definition: ex.h:135
ex series(const ex &r, int order, unsigned options=0) const
Compute the truncated series expansion of an expression.
Definition: pseries.cpp:1259
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
bool info(unsigned inf) const
Definition: ex.h:132
bool is_zero() const
Definition: ex.h:213
void print(const print_context &c, unsigned level=0) const
Print expression to stream.
Definition: ex.cpp:56
ex op(size_t i) const
Definition: ex.h:136
A pair of expressions.
Definition: expair.h:38
static unsigned register_new(function_options const &opt)
Definition: function.cpp:2243
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
bool info(unsigned inf) const override
Information about the object.
Definition: numeric.cpp:684
bool is_integer() const
True if object is a non-complex integer.
Definition: numeric.cpp:1154
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
bool is_equal(const numeric &other) const
Definition: numeric.cpp:1122
int to_int() const
Converts numeric types to machine's int.
Definition: numeric.cpp:1303
bool is_zero() const
True if object is zero.
Definition: numeric.cpp:1129
This class holds a two-component object, a basis and and exponent representing exponentiation.
Definition: power.h:39
Base class for print_contexts.
Definition: print.h:103
This class holds a extended truncated power series (positive and negative integer powers).
Definition: pseries.h:36
This class holds a relation consisting of two expressions and a logical relation between them.
Definition: relational.h:35
@ no_pattern
disable pattern matching
Definition: flags.h:51
Basic CAS symbol.
Definition: symbol.h:39
static unsigned serial
Definition: inifcns.h:109
static unsigned serial
Definition: inifcns.h:115
Interface to GiNaC's constant types and some special constants.
static const bool value
Definition: factor.cpp:231
unsigned options
Definition: factor.cpp:2480
vector< int > k
Definition: factor.cpp:1466
ex x
Definition: factor.cpp:1641
size_t n
Definition: factor.cpp:1463
size_t c
Definition: factor.cpp:770
size_t r
Definition: factor.cpp:770
umodpoly one
Definition: factor.cpp:1462
mvec m
Definition: factor.cpp:771
size_t last
Definition: factor.cpp:1465
Interface to GiNaC's initially known functions.
int order
Definition of GiNaC's lst.
Interface to GiNaC's products of expressions.
Definition: add.cpp:38
static ex G3_eval(const ex &x_, const ex &s_, const ex &y)
const numeric I
Imaginary unit.
Definition: numeric.cpp:1433
static ex zeta2_evalf(const ex &x, const ex &s)
static ex Li_evalf(const ex &m_, const ex &x_)
const ex _ex2
Definition: utils.cpp:389
const numeric pow(const numeric &x, const numeric &y)
Definition: numeric.h:251
const ex _ex_1_2
Definition: utils.cpp:356
static ex H_deriv(const ex &m_, const ex &x, unsigned deriv_param)
static ex S_eval(const ex &n, const ex &p, const ex &x)
container< std::list > lst
Definition: lst.h:32
static ex Li_deriv(const ex &m_, const ex &x_, unsigned deriv_param)
std::map< ex, ex, ex_is_less > exmap
Definition: basic.h:50
const numeric bernoulli(const numeric &nn)
Bernoulli number.
Definition: numeric.cpp:2166
std::vector< expair > epvector
expair-vector
Definition: expairseq.h:33
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
function zeta(const T1 &p1)
Definition: inifcns.h:111
const ex _ex1
Definition: utils.cpp:385
static void S_print_latex(const ex &n, const ex &p, const ex &x, const print_context &c)
static ex S_evalf(const ex &n, const ex &p, const ex &x)
ex conjugate(const ex &thisex)
Definition: ex.h:718
static ex zeta2_eval(const ex &m, const ex &s_)
const numeric imag(const numeric &x)
Definition: numeric.h:314
const numeric binomial(const numeric &n, const numeric &k)
The Binomial coefficients.
Definition: numeric.cpp:2143
const numeric * _num2_p
Definition: utils.cpp:388
const numeric exp(const numeric &x)
Exponential function.
Definition: numeric.cpp:1439
static ex S_deriv(const ex &n, const ex &p, const ex &x, unsigned deriv_param)
static ex Li_series(const ex &m, const ex &x, const relational &rel, int order, unsigned options)
const numeric factorial(const numeric &n)
Factorial combinatorial function.
Definition: numeric.cpp:2113
const ex _ex_1
Definition: utils.cpp:352
unsigned log2(unsigned n)
Integer binary logarithm.
Definition: utils.cpp:48
static void H_print_latex(const ex &m_, const ex &x, const print_context &c)
const constant Pi("Pi", PiEvalf, "\\pi", domain::positive)
Pi.
Definition: constant.h:82
static ex Li_eval(const ex &m_, const ex &x_)
static ex H_series(const ex &m, const ex &x, const relational &rel, int order, unsigned options)
static void Li_print_latex(const ex &m_, const ex &x_, const print_context &c)
static ex zeta1_eval(const ex &m)
static ex G2_eval(const ex &x_, const ex &y)
const numeric log(const numeric &x)
Natural logarithm.
Definition: numeric.cpp:1450
ex evalf(const ex &thisex)
Definition: ex.h:769
_numeric_digits Digits
Accuracy in decimal digits.
Definition: numeric.cpp:2586
bool is_real(const numeric &x)
Definition: numeric.h:293
ex op(const ex &thisex, size_t i)
Definition: ex.h:811
const numeric * _num1_p
Definition: utils.cpp:384
static void zeta1_print_latex(const ex &m_, const print_context &c)
const constant Catalan("Catalan", CatalanEvalf, "G", domain::positive)
Catalan's constant.
Definition: constant.h:83
static ex zeta2_deriv(const ex &m, const ex &s, unsigned deriv_param)
static void zeta2_print_latex(const ex &m_, const ex &s_, const print_context &c)
static ex G2_evalf(const ex &x_, const ex &y)
static ex S_series(const ex &n, const ex &p, const ex &x, const relational &rel, int order, unsigned options)
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
Definition: factor.cpp:2581
int to_int(const numeric &x)
Definition: numeric.h:302
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...
static ex H_evalf(const ex &x1, const ex &x2)
const ex _ex_48
Definition: utils.cpp:280
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"))
static ex G3_evalf(const ex &x_, const ex &s_, const ex &y)
const ex _ex0
Definition: utils.cpp:369
static ex zeta1_evalf(const ex &x)
std::vector< ex > exvector
Definition: basic.h:46
static ex zeta1_deriv(const ex &m, unsigned deriv_param)
static ex H_eval(const ex &m_, const ex &x)
numeric step(const numeric &x)
Definition: numeric.h:257
function G(const T1 &x, const T2 &y)
Definition: inifcns.h:130
ex expand(const ex &thisex, unsigned options=0)
Definition: ex.h:715
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to class for extended truncated power series.
Interface to relations between expressions.
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...
Interface to GiNaC's wildcard objects.

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