106 std::vector<std::vector<cln::cl_N>> Xn;
108 const int xninitsizestep = 26;
109 int xninitsize = xninitsizestep;
127 std::vector<cln::cl_N> buf(xninitsize);
128 auto it = buf.begin();
130 *it = -(cln::expt(cln::cl_I(2),
n+1) - 1) / cln::expt(cln::cl_I(2),
n+1);
132 for (
int i=2; i<=xninitsize; i++) {
136 result = Xn[0][i/2-1];
138 for (
int k=1;
k<i-1;
k++) {
139 if ( !(((i-
k) & 1) && ((i-
k) > 1)) ) {
144 result = result + Xn[
n-1][i-1] / (i+1);
152 std::vector<cln::cl_N> buf(xninitsize);
153 auto it = buf.begin();
155 *it = cln::cl_I(-3)/cln::cl_I(4);
157 *it = cln::cl_I(17)/cln::cl_I(36);
159 for (
int i=3; i<=xninitsize; i++) {
161 result = -Xn[0][(i-3)/2]/2;
165 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
166 for (
int k=1;
k<i/2;
k++) {
176 std::vector<cln::cl_N> buf(xninitsize/2);
177 auto it = buf.begin();
178 for (
int i=1; i<=xninitsize/2; i++) {
191 const int pos0 = xninitsize / 2;
193 for (
int i=1; i<=xninitsizestep/2; ++i) {
194 Xn[0].push_back(
bernoulli((i+pos0)*2).to_cl_N());
197 int xend = xninitsize + xninitsizestep;
200 for (
int i=xninitsize+1; i<=xend; ++i) {
202 result = -Xn[0][(i-3)/2]/2;
205 result = Xn[0][i/2-1] + Xn[0][i/2-1]/(i+1);
206 for (
int k=1;
k<i/2;
k++) {
209 Xn[1].push_back(result);
213 for (
size_t n=2;
n<Xn.size(); ++
n) {
214 for (
int i=xninitsize+1; i<=xend; ++i) {
218 result = Xn[0][i/2-1];
220 for (
int k=1;
k<i-1; ++
k) {
221 if ( !(((i-
k) & 1) && ((i-
k) > 1)) ) {
226 result = result + Xn[
n-1][i-1] / (i+1);
227 Xn[
n].push_back(result);
231 xninitsize += xninitsizestep;
236 cln::cl_N Li2_do_sum(
const cln::cl_N&
x)
240 cln::cl_N num =
x * cln::cl_float(1, cln::float_format(
Digits));
248 res = res + num / den;
249 }
while (res != resbuf);
255 cln::cl_N Li2_do_sum_Xn(
const cln::cl_N&
x)
257 std::vector<cln::cl_N>::const_iterator it = Xn[0].
begin();
258 std::vector<cln::cl_N>::const_iterator xend = Xn[0].end();
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;
268 res = res + (*it) *
factor;
272 it = Xn[0].begin() + (i-1);
275 }
while (res != resbuf);
281 cln::cl_N Lin_do_sum(
int n,
const cln::cl_N&
x)
283 cln::cl_N
factor =
x * cln::cl_float(1, cln::float_format(
Digits));
290 res = res +
factor / cln::expt(cln::cl_I(i),
n);
292 }
while (res != resbuf);
298 cln::cl_N Lin_do_sum_Xn(
int n,
const cln::cl_N&
x)
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();
303 cln::cl_N
factor = u * cln::cl_float(1, cln::float_format(
Digits));
310 res = res + (*it) *
factor;
314 it = Xn[
n-2].begin() + (i-2);
315 xend = Xn[
n-2].end();
317 }
while (res != resbuf);
323 const cln::cl_N S_num(
int n,
int p,
const cln::cl_N&
x);
327 cln::cl_N Li_projection(
int n,
const cln::cl_N&
x,
const cln::float_format_t& prec)
336 if (cln::realpart(
x) < 0.5) {
342 return Li2_do_sum(
x);
345 return Li2_do_sum_Xn(
x);
362 for (
int i=xnsize; i<
n-1; i++) {
367 if (cln::realpart(
x) < 0.5) {
371 return Lin_do_sum(
n,
x);
374 return Lin_do_sum_Xn(
n,
x);
377 cln::cl_N result = 0;
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))
389 const cln::cl_N Lin_numeric(
const int n,
const cln::cl_N&
x)
404 return -(1-cln::expt(cln::cl_I(2),1-
n)) *
cln::zeta(
n);
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))
417 cln::float_format_t prec = cln::default_float_format;
418 const cln::cl_N
value =
x;
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)));
429 if (cln::zerop(cln::imagpart(
value))) {
431 result = result +
conjugate(Li_projection(
n, cln::recip(
value), prec));
434 result = result -
conjugate(Li_projection(
n, cln::recip(
value), prec));
439 result = result + Li_projection(
n, cln::recip(
value), prec);
442 result = result - Li_projection(
n, cln::recip(
value), prec);
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))
450 result = result - add;
454 return Li_projection(
n,
value, prec);
476 cln::cl_N multipleLi_do_sum(
const std::vector<int>& s,
const std::vector<cln::cl_N>&
x)
479 for (
const auto & it :
x) {
480 if (it == 0)
return cln::cl_float(0, cln::float_format(
Digits));
483 const int j = s.size();
484 bool flag_accidental_zero =
false;
486 std::vector<cln::cl_N> t(j);
487 cln::cl_F
one = cln::cl_float(1, cln::float_format(
Digits));
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]);
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]);
504 }
while ( (t[0] != t0buf) || cln::zerop(t[0]) || flag_accidental_zero );
511 lst convert_parameter_Li_to_H(
const lst&
m,
const lst&
x, ex& pf);
515 typedef std::vector<int> Gparameter;
519 ex G_eval1(
int a,
int scale,
const exvector& gsyms)
522 const ex& scs = gsyms[
std::abs(scale)];
525 return -
log(1 - scs/as);
536 ex G_eval(
const Gparameter& a,
int scale,
const exvector& gsyms)
541 bool all_zero =
true;
542 bool all_ones =
true;
544 for (
const auto & it : a) {
562 if (newa.nops() > 1 && newa.op(0) == sc && !all_ones && a.front()!=0) {
564 Gparameter short_a(a.begin()+1, a.end());
565 ex result = G_eval1(a.front(), scale, gsyms) * G_eval(short_a, scale, gsyms);
567 auto it = short_a.begin();
568 advance(it, count_ones-1);
569 for (; it != short_a.end(); ++it) {
571 Gparameter newa(short_a.begin(), it);
573 newa.push_back(a[0]);
574 newa.insert(newa.end(), it+1, short_a.end());
575 result -= G_eval(newa, scale, gsyms);
577 return result / count_ones;
581 if (all_ones && a.size() > 1) {
582 return pow(G_eval1(a.front(),scale, gsyms), count_ones) /
factorial(count_ones);
595 for (
const auto & it : a) {
597 const ex& sym = gsyms[
std::abs(it)];
598 x.append(argbuf / sym);
610 ex G_eval_to_G(
const Gparameter& a,
int scale,
const exvector& gsyms)
614 for (
const auto & it : a) {
632 Gparameter convert_pending_integrals_G(
const Gparameter& pending_integrals)
636 if (pending_integrals.size() > 0) {
638 Gparameter new_a(pending_integrals.begin()+1, pending_integrals.end());
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)
660 auto lastnonzero = a.end();
661 for (
auto it = a.begin(); it != a.end(); ++it) {
676 if (lastnonzero == a.end())
678 return ++lastnonzero;
683 Gparameter prepare_pending_integrals(
const Gparameter& pending_integrals,
int scale)
687 if (pending_integrals.size() > 0) {
688 return pending_integrals;
690 Gparameter new_pending_integrals;
691 new_pending_integrals.push_back(scale);
692 return new_pending_integrals;
698 ex trailing_zeros_G(
const Gparameter& a,
int scale,
const exvector& gsyms)
701 int depth, trailing_zeros;
702 Gparameter::const_iterator
last, dummyit;
703 last = check_parameter_G(a, scale, convergent, depth, trailing_zeros, dummyit);
707 if ((trailing_zeros > 0) && (depth > 0)) {
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);
714 new_a.insert(new_a.end(), it, a.end()-1);
715 result -= trailing_zeros_G(new_a, scale, gsyms);
718 return result / trailing_zeros;
720 return G_eval(a, scale, gsyms);
726 ex depth_one_trafo_G(
const Gparameter& pending_integrals,
const Gparameter& a,
int scale,
const exvector& gsyms)
739 Gparameter new_pending_integrals = prepare_pending_integrals(pending_integrals,
std::abs(a.back()));
740 const int psize = pending_integrals.size();
748 result +=
log(gsyms[ex_to<numeric>(scale).
to_int()]);
750 new_pending_integrals.push_back(-scale);
753 new_pending_integrals.push_back(scale);
757 result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals),
758 pending_integrals.front(),
763 result += trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals),
764 new_pending_integrals.front(),
768 new_pending_integrals.back() = 0;
769 result -= trailing_zeros_G(convert_pending_integrals_G(new_pending_integrals),
770 new_pending_integrals.front(),
781 result -=
zeta(a.size());
783 result *= trailing_zeros_G(convert_pending_integrals_G(pending_integrals),
784 pending_integrals.front(),
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);
796 Gparameter new_pending_integrals_2;
797 new_pending_integrals_2.push_back(scale);
798 new_pending_integrals_2.push_back(0);
800 result += trailing_zeros_G(convert_pending_integrals_G(pending_integrals),
801 pending_integrals.front(),
803 * depth_one_trafo_G(new_pending_integrals_2, new_a, scale, gsyms);
805 result += depth_one_trafo_G(new_pending_integrals_2, new_a, scale, gsyms);
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);
819 ex G_transform(
const Gparameter& pendint,
const Gparameter& a,
int scale,
820 const exvector& gsyms,
bool flag_trailing_zeros_only)
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);
845 result = G_eval(a, scale, gsyms);
847 if (pendint.size() > 0) {
848 result *= trailing_zeros_G(convert_pending_integrals_G(pendint),
856 if (trailing_zeros > 0) {
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);
863 new_a.insert(new_a.end(), it, a.end()-1);
864 result -= G_transform(pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
866 return result / trailing_zeros;
870 if (flag_trailing_zeros_only)
871 return G_eval_to_G(a, scale, gsyms);
875 if (pendint.size() > 0) {
876 return G_eval(convert_pending_integrals_G(pendint),
877 pendint.front(), gsyms) *
878 G_eval(a, scale, gsyms);
880 return G_eval(a, scale, gsyms);
886 return depth_one_trafo_G(pendint, a, scale, gsyms);
895 if (min_it + 1 == a.end()) {
896 do { --min_it; }
while (*min_it == 0);
898 Gparameter a1(a.begin(),min_it+1);
899 Gparameter a2(min_it+1,a.end());
901 ex result = G_transform(pendint, a2, scale, gsyms, flag_trailing_zeros_only)*
902 G_transform(empty, a1, scale, gsyms, flag_trailing_zeros_only);
904 result -= shuffle_G(empty, a1, a2, pendint, a, scale, gsyms, flag_trailing_zeros_only);
909 Gparameter::iterator changeit;
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);
922 changeit = new_a.begin() + min_it_pos;
923 changeit = new_a.erase(changeit);
924 if (changeit != new_a.begin()) {
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;
932 result += G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
934 new_pendint.pop_back();
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);
941 result -= G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
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);
953 result += G_transform(new_pendint, new_a, scale, gsyms, flag_trailing_zeros_only);
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)
965 if (a1.size()==0 && a2.size()==0) {
967 if ( a0 == a_old )
return 0;
969 return G_transform(pendint, a0, scale, gsyms, flag_trailing_zeros_only);
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);
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);
986 Gparameter a1_removed(a1.begin()+1,a1.end());
987 Gparameter a2_removed(a2.begin()+1,a2.end());
992 a01.push_back( a1[0] );
993 a02.push_back( a2[0] );
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);
1002 G_numeric(
const std::vector<cln::cl_N>&
x,
const std::vector<int>& s,
1003 const cln::cl_N& y);
1008 G_do_hoelder(std::vector<cln::cl_N>
x,
1009 const std::vector<int>& s,
const cln::cl_N& y)
1012 const std::size_t size =
x.size();
1013 for (std::size_t i = 0; i < size; ++i)
1016 for (std::size_t
r = 0;
r <= size; ++
r) {
1017 cln::cl_N buffer(1 &
r ? -1 : 1);
1022 for (std::size_t i = 0; i < size; ++i) {
1023 if (
x[i] == cln::cl_RA(1)/p) {
1024 p = p/2 + cln::cl_RA(3)/2;
1030 cln::cl_RA q = p/(p-1);
1031 std::vector<cln::cl_N> qlstx;
1032 std::vector<int> qlsts;
1033 for (std::size_t j =
r; j >= 1; --j) {
1034 qlstx.push_back(cln::cl_N(1) -
x[j-1]);
1035 if (imagpart(
x[j-1])==0 && realpart(
x[j-1]) >= 1) {
1038 qlsts.push_back(-s[j-1]);
1041 if (qlstx.size() > 0) {
1042 buffer = buffer*G_numeric(qlstx, qlsts, 1/q);
1044 std::vector<cln::cl_N> plstx;
1045 std::vector<int> plsts;
1046 for (std::size_t j =
r+1; j <= size; ++j) {
1047 plstx.push_back(
x[j-1]);
1048 plsts.push_back(s[j-1]);
1050 if (plstx.size() > 0) {
1051 buffer = buffer*G_numeric(plstx, plsts, 1/p);
1053 result = result + buffer;
1058 class less_object_for_cl_N
1061 bool operator() (
const cln::cl_N & a,
const cln::cl_N & b)
const
1065 return (
abs(a) <
abs(b)) ? true :
false;
1068 if (phase(a) != phase(b))
1069 return (phase(a) < phase(b)) ? true :
false;
1080 G_do_trafo(
const std::vector<cln::cl_N>&
x,
const std::vector<int>& s,
1081 const cln::cl_N& y,
bool flag_trailing_zeros_only)
1084 typedef std::multimap<cln::cl_N, std::size_t, less_object_for_cl_N> sortmap_t;
1086 std::size_t size = 0;
1087 for (std::size_t i = 0; i <
x.size(); ++i) {
1089 sortmap.insert(std::make_pair(
x[i], i));
1094 sortmap.insert(std::make_pair(y,
x.size()));
1100 gsyms.push_back(symbol(
"GSYMS_ERROR"));
1101 cln::cl_N lastentry(0);
1102 for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
1103 if (it != sortmap.begin()) {
1104 if (it->second <
x.size()) {
1105 if (
x[it->second] == lastentry) {
1106 gsyms.push_back(gsyms.back());
1110 if (y == lastentry) {
1111 gsyms.push_back(gsyms.back());
1116 std::ostringstream os;
1118 gsyms.push_back(symbol(os.str()));
1120 if (it->second <
x.size()) {
1121 lastentry =
x[it->second];
1128 Gparameter a(
x.size());
1130 std::size_t pos = 1;
1132 for (sortmap_t::const_iterator it = sortmap.begin(); it != sortmap.end(); ++it) {
1133 if (it->second <
x.size()) {
1134 if (s[it->second] > 0) {
1135 a[it->second] = pos;
1137 a[it->second] = -int(pos);
1139 subslst[gsyms[pos]] = numeric(
x[it->second]);
1142 subslst[gsyms[pos]] = numeric(y);
1149 ex result = G_transform(pendint, a, scale, gsyms, flag_trailing_zeros_only);
1151 result = result.expand();
1152 result = result.subs(subslst).evalf();
1153 if (!is_a<numeric>(result))
1154 throw std::logic_error(
"G_do_trafo: G_transform returned non-numeric result");
1156 cln::cl_N ret = ex_to<numeric>(result).to_cl_N();
1163 G_numeric(
const std::vector<cln::cl_N>&
x,
const std::vector<int>& s,
1167 bool need_trafo =
false;
1168 bool need_hoelder =
false;
1169 bool have_trailing_zero =
false;
1170 std::size_t depth = 0;
1171 for (
auto & xi :
x) {
1174 const cln::cl_N x_y =
abs(xi) - y;
1175 if (instanceof(x_y, cln::cl_R_ring) &&
1176 realpart(x_y) < cln::least_negative_float(cln::float_format(
Digits - 2)))
1179 if (
abs(
abs(xi/y) - 1) < 0.01)
1180 need_hoelder =
true;
1183 if (zerop(
x.back())) {
1184 have_trailing_zero =
true;
1188 if (depth == 1 &&
x.size() == 2 && !need_trafo)
1189 return - Li_projection(2, y/
x[1], cln::float_format(
Digits));
1192 if (need_hoelder && !have_trailing_zero)
1193 return G_do_hoelder(
x, s, y);
1197 return G_do_trafo(
x, s, y, have_trailing_zero);
1200 std::vector<cln::cl_N> newx;
1201 newx.reserve(
x.size());
1203 m.reserve(
x.size());
1207 for (
auto & xi :
x) {
1211 newx.push_back(
factor/xi);
1213 m.push_back(mcount);
1219 return sign*multipleLi_do_sum(
m, newx);
1223 ex mLi_numeric(
const lst&
m,
const lst&
x)
1226 std::vector<cln::cl_N> newx;
1227 newx.reserve(
x.
nops());
1229 s.reserve(
x.
nops());
1231 for (
auto itm =
m.begin(), itx =
x.
begin(); itm !=
m.end(); ++itm, ++itx) {
1232 for (
int i = 1; i < *itm; ++i) {
1233 newx.push_back(cln::cl_N(0));
1236 const cln::cl_N xi = ex_to<numeric>(*itx).to_cl_N();
1239 if ( !instanceof(
factor, cln::cl_R_ring) && imagpart(
factor) < 0 ) {
1246 return numeric(cln::cl_N(1 &
m.nops() ? - 1 : 1)*G_numeric(newx, s, cln::cl_N(1)));
1265 return G(x_, y).
hold();
1267 lst x = is_a<lst>(x_) ? ex_to<lst>(x_) :
lst{x_};
1268 if (
x.
nops() == 0) {
1272 return G(x_, y).
hold();
1275 s.reserve(
x.
nops());
1276 bool all_zero =
true;
1277 for (
const auto & it :
x) {
1279 return G(x_, y).
hold();
1284 if ( !ex_to<numeric>(it).
is_real() && ex_to<numeric>(it).
imag() < 0 ) {
1294 std::vector<cln::cl_N> xv;
1295 xv.reserve(
x.
nops());
1296 for (
const auto & it :
x)
1297 xv.push_back(ex_to<numeric>(it).to_cl_N());
1298 cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
1308 return G(x_, y).
hold();
1310 lst x = is_a<lst>(x_) ? ex_to<lst>(x_) :
lst{x_};
1311 if (
x.
nops() == 0) {
1315 return G(x_, y).
hold();
1318 s.reserve(
x.
nops());
1319 bool all_zero =
true;
1320 bool crational =
true;
1321 for (
const auto & it :
x) {
1323 return G(x_, y).
hold();
1331 if ( !ex_to<numeric>(it).
is_real() && ex_to<numeric>(it).
imag() < 0 ) {
1345 return G(x_, y).
hold();
1347 std::vector<cln::cl_N> xv;
1348 xv.reserve(
x.
nops());
1349 for (
const auto & it :
x)
1350 xv.push_back(ex_to<numeric>(it).to_cl_N());
1351 cln::cl_N result = G_numeric(xv, s, ex_to<numeric>(y).to_cl_N());
1369 return G(x_, s_, y).
hold();
1371 lst x = is_a<lst>(x_) ? ex_to<lst>(x_) :
lst{x_};
1372 lst s = is_a<lst>(s_) ? ex_to<lst>(s_) :
lst{s_};
1374 return G(x_, s_, y).
hold();
1376 if (
x.
nops() == 0) {
1380 return G(x_, s_, y).
hold();
1382 std::vector<int> sn;
1383 sn.reserve(s.
nops());
1384 bool all_zero =
true;
1385 for (
auto itx =
x.
begin(), its = s.
begin(); itx !=
x.
end(); ++itx, ++its) {
1387 return G(x_, y).
hold();
1390 return G(x_, y).
hold();
1395 if ( ex_to<numeric>(*itx).is_real() ) {
1396 if ( ex_to<numeric>(*itx).is_positive() ) {
1408 if ( ex_to<numeric>(*itx).imag() > 0 ) {
1419 std::vector<cln::cl_N> xn;
1420 xn.reserve(
x.
nops());
1421 for (
const auto & it :
x)
1422 xn.push_back(ex_to<numeric>(it).to_cl_N());
1423 cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
1433 return G(x_, s_, y).
hold();
1435 lst x = is_a<lst>(x_) ? ex_to<lst>(x_) :
lst{x_};
1436 lst s = is_a<lst>(s_) ? ex_to<lst>(s_) :
lst{s_};
1438 return G(x_, s_, y).
hold();
1440 if (
x.
nops() == 0) {
1444 return G(x_, s_, y).
hold();
1446 std::vector<int> sn;
1447 sn.reserve(s.
nops());
1448 bool all_zero =
true;
1449 bool crational =
true;
1450 for (
auto itx =
x.
begin(), its = s.
begin(); itx !=
x.
end(); ++itx, ++its) {
1452 return G(x_, s_, y).
hold();
1455 return G(x_, s_, y).
hold();
1463 if ( ex_to<numeric>(*itx).is_real() ) {
1464 if ( ex_to<numeric>(*itx).is_positive() ) {
1476 if ( ex_to<numeric>(*itx).imag() > 0 ) {
1491 return G(x_, s_, y).
hold();
1493 std::vector<cln::cl_N> xn;
1494 xn.reserve(
x.
nops());
1495 for (
const auto & it :
x)
1496 xn.push_back(ex_to<numeric>(it).to_cl_N());
1497 cln::cl_N result = G_numeric(xn, sn, ex_to<numeric>(y).to_cl_N());
1528 int m__ = ex_to<numeric>(m_).to_int();
1529 const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
1530 const cln::cl_N result = Lin_numeric(m__, x__);
1536 int m__ = ex_to<numeric>(m_).to_int();
1537 const cln::cl_N x__ = ex_to<numeric>(x_val).to_cl_N();
1538 const cln::cl_N result = Lin_numeric(m__, x__);
1544 if (is_a<lst>(m_) && is_a<lst>(x_)) {
1546 const lst&
m = ex_to<lst>(m_);
1547 const lst&
x = ex_to<lst>(x_);
1548 if (
m.nops() !=
x.
nops()) {
1549 return Li(m_,x_).hold();
1551 if (
x.
nops() == 0) {
1555 return Li(m_,x_).hold();
1558 for (
auto itm =
m.begin(), itx =
x.
begin(); itm !=
m.end(); ++itm, ++itx) {
1560 return Li(m_, x_).hold();
1563 return Li(m_, x_).hold();
1570 return mLi_numeric(
m,
x);
1573 return Li(m_,x_).hold();
1579 if (is_a<lst>(m_)) {
1580 if (is_a<lst>(x_)) {
1582 const lst&
m = ex_to<lst>(m_);
1583 const lst&
x = ex_to<lst>(x_);
1584 if (
m.nops() !=
x.
nops()) {
1585 return Li(m_,x_).hold();
1587 if (
x.
nops() == 0) {
1591 bool is_zeta =
true;
1592 bool do_evalf =
true;
1593 bool crational =
true;
1594 for (
auto itm =
m.begin(), itx =
x.
begin(); itm !=
m.end(); ++itm, ++itx) {
1596 return Li(m_,x_).hold();
1598 if ((*itx !=
_ex1) && (*itx !=
_ex_1)) {
1616 for (
const auto & itx :
x) {
1625 return zeta(m_, newx);
1629 lst newm = convert_parameter_Li_to_H(
m,
x, prefactor);
1630 return prefactor * H(newm,
x[0]);
1632 if (do_evalf && !crational) {
1633 return mLi_numeric(
m,
x);
1636 return Li(m_, x_).hold();
1637 }
else if (is_a<lst>(x_)) {
1638 return Li(m_, x_).hold();
1649 return (
pow(2,1-m_)-1) *
zeta(m_);
1663 int m__ = ex_to<numeric>(m_).
to_int();
1664 const cln::cl_N x__ = ex_to<numeric>(x_).to_cl_N();
1665 const cln::cl_N result = Lin_numeric(m__, x__);
1669 return Li(m_, x_).hold();
1675 if (is_a<lst>(
m) || is_a<lst>(
x)) {
1678 return pseries(rel, std::move(seq));
1689 for (
int i=1; i<
order; ++i)
1695 ser +=
pseries(rel, std::move(nseq));
1700 throw std::runtime_error(
"Li_series: don't know how to do the series expansion at this point!");
1710 if (deriv_param == 0) {
1713 if (m_.
nops() > 1) {
1714 throw std::runtime_error(
"don't know how to derivate multiple polylogarithm!");
1717 if (is_a<lst>(m_)) {
1723 if (is_a<lst>(x_)) {
1729 return Li(
m-1,
x) /
x;
1739 if (is_a<lst>(m_)) {
1745 if (is_a<lst>(x_)) {
1750 c.s <<
"\\mathrm{Li}_{";
1751 auto itm =
m.begin();
1754 for (; itm !=
m.end(); itm++) {
1762 for (; itx !=
x.
end(); itx++) {
1776 do_not_evalf_params());
1794 std::vector<std::vector<cln::cl_N>> Yn;
1807 void fill_Yn(
int n,
const cln::float_format_t& prec)
1809 const int initsize = ynlength;
1811 cln::cl_N
one = cln::cl_float(1, prec);
1814 std::vector<cln::cl_N> buf(initsize);
1815 auto it = buf.begin();
1816 auto itprev = Yn[
n-1].begin();
1817 *it = (*itprev) / cln::cl_N(
n+1) *
one;
1822 for (
int i=
n+2; i<=initsize+
n; i++) {
1823 *it = *(it-1) + (*itprev) / cln::cl_N(i) *
one;
1829 std::vector<cln::cl_N> buf(initsize);
1830 auto it = buf.begin();
1833 for (
int i=2; i<=initsize; i++) {
1834 *it = *(it-1) + 1 / cln::cl_N(i) *
one;
1844 void make_Yn_longer(
int newsize,
const cln::float_format_t& prec)
1847 cln::cl_N
one = cln::cl_float(1, prec);
1849 Yn[0].resize(newsize);
1850 auto it = Yn[0].begin();
1852 for (
int i=ynlength+1; i<=newsize; i++) {
1853 *it = *(it-1) + 1 / cln::cl_N(i) *
one;
1857 for (
int n=1;
n<ynsize;
n++) {
1858 Yn[
n].resize(newsize);
1859 auto it = Yn[
n].begin();
1860 auto itprev = Yn[
n-1].begin();
1863 for (
int i=ynlength+
n+1; i<=newsize+
n; i++) {
1864 *it = *(it-1) + (*itprev) / cln::cl_N(i) *
one;
1876 cln::cl_N C(
int n,
int p)
1880 for (
int k=0;
k<p;
k++) {
1881 for (
int j=0; j<=(
n+
k-1)/2; j++) {
1885 result = result - 2 * cln::expt(cln::pi(),2*j) * S_num(
n-2*j,p,1) /
cln::factorial(2*j);
1888 result = result + 2 * cln::expt(cln::pi(),2*j) * S_num(
n-2*j,p,1) /
cln::factorial(2*j);
1896 * cln::expt(cln::pi(),2*j) * S_num(
n+
k-2*j,p-
k,1)
1901 * cln::expt(cln::pi(),2*j) * S_num(
n+
k-2*j,p-
k,1)
1912 * cln::expt(cln::pi(),2*j) * S_num(
n+
k-2*j,p-
k,1)
1921 if (((np)/2+
n) & 1) {
1935 cln::cl_N a_k(
int k)
1944 for (
int m=2;
m<=
k;
m++) {
1945 result = result + cln::expt(cln::cl_N(-1),
m) *
cln::zeta(
m) * a_k(
k-
m);
1954 cln::cl_N b_k(
int k)
1963 for (
int m=2;
m<=
k;
m++) {
1964 result = result + cln::expt(cln::cl_N(-1),
m) *
cln::zeta(
m) * b_k(
k-
m);
1972 cln::cl_N S_do_sum(
int n,
int p,
const cln::cl_N&
x,
const cln::float_format_t& prec)
1974 static cln::float_format_t oldprec = cln::default_float_format;
1977 return Li_projection(
n+1,
x, prec);
1981 if ( oldprec != prec ) {
1990 for (
int i=ynsize; i<p-1; i++) {
1996 cln::cl_F
one = cln::cl_float(1, cln::float_format(
Digits));
1997 cln::cl_N xf =
x *
one;
2002 cln::cl_N
factor = cln::expt(xf, p);
2006 if (i-p >= ynlength) {
2008 make_Yn_longer(ynlength*2, prec);
2010 res = res +
factor / cln::expt(cln::cl_I(i),
n+1) * Yn[p-2][i-p];
2014 }
while (res != resbuf);
2021 cln::cl_N S_projection(
int n,
int p,
const cln::cl_N&
x,
const cln::float_format_t& prec)
2024 if (
cln::abs(cln::realpart(
x)) > cln::cl_F(
"0.5")) {
2026 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(
cln::log(
x),
n)
2029 for (
int s=0; s<
n; s++) {
2031 for (
int r=0;
r<p;
r++) {
2032 res2 = res2 + cln::expt(cln::cl_I(-1),
r) * cln::expt(
cln::log(1-
x),
r)
2041 return S_do_sum(
n, p,
x, prec);
2046 const cln::cl_N S_num(
int n,
int p,
const cln::cl_N&
x)
2061 for (
int nu=0; nu<
n; nu++) {
2062 for (
int rho=0; rho<=p; rho++) {
2063 result = result + b_k(
n-nu-1) * b_k(p-rho) * a_k(nu+rho+1)
2067 result = result * cln::expt(cln::cl_I(-1),
n+p-1);
2074 return -(1-cln::expt(cln::cl_I(2),-
n)) *
cln::zeta(
n+1);
2081 cln::float_format_t prec = cln::default_float_format;
2082 const cln::cl_N
value =
x;
2084 if (!instanceof(realpart(
value), cln::cl_RA_ring))
2085 prec = cln::float_format(cln::the<cln::cl_F>(cln::realpart(
value)));
2086 else if (!instanceof(imagpart(
value), cln::cl_RA_ring))
2087 prec = cln::float_format(cln::the<cln::cl_F>(cln::imagpart(
value)));
2094 cln::cl_N result = cln::expt(cln::cl_I(-1),p) * cln::expt(
cln::log(
value),
n)
2097 for (
int s=0; s<
n; s++) {
2099 for (
int r=0;
r<p;
r++) {
2100 res2 = res2 + cln::expt(cln::cl_I(-1),
r) * cln::expt(
cln::log(1-
value),
r)
2114 for (
int s=0; s<p; s++) {
2115 for (
int r=0;
r<=s;
r++) {
2118 * S_num(
n+s-
r,p-s,cln::recip(
value));
2121 result = result * cln::expt(cln::cl_I(-1),
n);
2124 for (
int r=0;
r<
n;
r++) {
2129 result = result + cln::expt(cln::cl_I(-1),p) * res2;
2137 for (
int s=0; s<p-1; s++)
2140 ex res = H(
m,numeric(
value)).evalf();
2141 return ex_to<numeric>(res).to_cl_N();
2144 return S_projection(
n, p,
value, prec);
2164 const int n_ = ex_to<numeric>(
n).to_int();
2165 const int p_ = ex_to<numeric>(p).to_int();
2166 if (is_a<numeric>(
x)) {
2167 const cln::cl_N x_ = ex_to<numeric>(
x).to_cl_N();
2168 const cln::cl_N result = S_num(n_, p_, x_);
2172 if (is_a<numeric>(x_val)) {
2173 const cln::cl_N x_val_ = ex_to<numeric>(x_val).to_cl_N();
2174 const cln::cl_N result = S_num(n_, p_, x_val_);
2179 return S(
n, p,
x).hold();
2191 for (
int i=ex_to<numeric>(p).
to_int()-1; i>0; i--) {
2200 int n_ = ex_to<numeric>(
n).to_int();
2201 int p_ = ex_to<numeric>(p).to_int();
2202 const cln::cl_N x_ = ex_to<numeric>(
x).to_cl_N();
2203 const cln::cl_N result = S_num(n_, p_, x_);
2211 return S(
n, p,
x).hold();
2230 std::vector<ex> presubsum, subsum;
2231 subsum.push_back(0);
2232 for (
int i=1; i<
order-1; ++i) {
2233 subsum.push_back(subsum[i-1] +
numeric(1, i));
2235 for (
int depth=2; depth<p; ++depth) {
2237 for (
int i=1; i<
order-1; ++i) {
2238 subsum[i] = subsum[i-1] +
numeric(1, i) * presubsum[i-1];
2242 for (
int i=1; i<
order; ++i) {
2249 ser +=
pseries(rel, std::move(nseq));
2254 throw std::runtime_error(
"S_series: don't know how to do the series expansion at this point!");
2264 if (deriv_param < 2) {
2268 return S(
n-1, p,
x) /
x;
2270 return S(
n, p-1,
x) / (1-
x);
2277 c.s <<
"\\mathrm{S}_{";
2293 do_not_evalf_params());
2310 symbol H_polesign(
"IMSIGN");
2316 bool convert_parameter_H_to_Li(
const lst& l,
lst&
m,
lst& s,
ex& pf)
2320 for (
const auto & it : l) {
2322 for (
ex count=it-1; count > 0; count--) {
2326 }
else if (it < -1) {
2327 for (ex count=it+1; count < 0; count++) {
2338 bool has_negative_parameters =
false;
2340 for (
const auto & it : mexp) {
2346 m.append((it+acc-1) * signum);
2348 m.append((it-acc+1) * signum);
2354 has_negative_parameters =
true;
2357 if (has_negative_parameters) {
2358 for (std::size_t i=0; i<
m.nops(); i++) {
2360 m.let_op(i) = -
m.op(i);
2368 return has_negative_parameters;
2373 struct map_trafo_H_convert_to_Li :
public map_function
2375 ex operator()(
const ex& e)
override
2377 if (is_a<add>(e) || is_a<mul>(e)) {
2378 return e.
map(*
this);
2380 if (is_a<function>(e)) {
2381 std::string name = ex_to<function>(e).get_name();
2384 if (is_a<lst>(e.op(0))) {
2385 parameter = ex_to<lst>(e.op(0));
2387 parameter =
lst{e.op(0)};
2394 if (convert_parameter_H_to_Li(parameter,
m, s, pf)) {
2396 return pf * Li(
m, s).hold();
2398 for (std::size_t i=0; i<
m.nops(); i++) {
2402 return Li(
m, s).hold();
2412 struct map_trafo_H_convert_to_zeta :
public map_function
2414 ex operator()(
const ex& e)
override
2416 if (is_a<add>(e) || is_a<mul>(e)) {
2417 return e.
map(*
this);
2419 if (is_a<function>(e)) {
2420 std::string name = ex_to<function>(e).get_name();
2423 if (is_a<lst>(e.op(0))) {
2424 parameter = ex_to<lst>(e.op(0));
2426 parameter =
lst{e.op(0)};
2432 if (convert_parameter_H_to_Li(parameter,
m, s, pf)) {
2433 return pf *
zeta(
m, s);
2445 struct map_trafo_H_reduce_trailing_zeros :
public map_function
2447 ex operator()(
const ex& e)
override
2449 if (is_a<add>(e) || is_a<mul>(e)) {
2450 return e.map(*
this);
2452 if (is_a<function>(e)) {
2453 std::string name = ex_to<function>(e).get_name();
2456 if (is_a<lst>(e.op(0))) {
2457 parameter = ex_to<lst>(e.op(0));
2459 parameter =
lst{e.op(0)};
2462 if (parameter.op(parameter.nops()-1) == 0) {
2465 if (parameter.nops() == 1) {
2470 auto it = parameter.begin();
2471 while ((it != parameter.end()) && (*it == 0)) {
2474 if (it == parameter.end()) {
2479 parameter.remove_last();
2480 std::size_t lastentry = parameter.nops();
2481 while ((lastentry > 0) && (parameter[lastentry-1] == 0)) {
2486 ex result =
log(arg) * H(parameter,arg).
hold();
2488 for (ex i=0; i<lastentry; i++) {
2489 if (parameter[i] > 0) {
2491 result -= (acc + parameter[i]-1) * H(parameter, arg).
hold();
2494 }
else if (parameter[i] < 0) {
2496 result -= (acc +
abs(parameter[i]+1)) * H(parameter, arg).
hold();
2504 if (lastentry < parameter.nops()) {
2505 result = result / (parameter.nops()-lastentry+1);
2506 return result.
map(*
this);
2519 ex convert_H_to_zeta(
const lst&
m)
2521 symbol xtemp(
"xtemp");
2522 map_trafo_H_reduce_trailing_zeros filter;
2523 map_trafo_H_convert_to_zeta filter2;
2524 return filter2(filter(H(
m, xtemp).hold())).subs(xtemp == 1);
2529 lst convert_parameter_Li_to_H(
const lst&
m,
const lst&
x, ex& pf)
2532 auto itm =
m.begin();
2538 while (itx !=
x.
end()) {
2544 signum *= (*itx !=
_ex_1) ? 1 : -1;
2546 res.append((*itm) * signum);
2556 ex trafo_H_mult(
const ex& h1,
const ex& h2)
2561 ex h1nops = h1.op(0).nops();
2562 ex h2nops = h2.op(0).nops();
2564 hshort = h2.op(0).op(0);
2565 hlong = ex_to<lst>(h1.op(0));
2567 hshort = h1.op(0).op(0);
2569 hlong = ex_to<lst>(h2.op(0));
2571 hlong =
lst{h2.op(0).op(0)};
2574 for (std::size_t i=0; i<=hlong.nops(); i++) {
2578 newparameter.append(hlong[j]);
2580 newparameter.append(hshort);
2581 for (; j<hlong.nops(); j++) {
2582 newparameter.append(hlong[j]);
2584 res += H(newparameter, h1.op(1)).hold();
2591 struct map_trafo_H_mult :
public map_function
2593 ex operator()(
const ex& e)
override
2596 return e.map(*
this);
2604 for (std::size_t pos=0; pos<e.nops(); pos++) {
2605 if (is_a<power>(e.op(pos)) && is_a<function>(e.op(pos).op(0))) {
2606 std::string name = ex_to<function>(e.op(pos).op(0)).get_name();
2608 for (ex i=0; i<e.op(pos).
op(1); i++) {
2609 Hlst.append(e.op(pos).op(0));
2613 }
else if (is_a<function>(e.op(pos))) {
2614 std::string name = ex_to<function>(e.op(pos)).get_name();
2616 if (e.op(pos).op(0).nops() > 1) {
2619 Hlst.append(e.op(pos));
2624 result *= e.op(pos);
2627 if (Hlst.nops() > 0) {
2628 firstH = Hlst[Hlst.nops()-1];
2635 if (Hlst.nops() > 0) {
2636 ex buffer = trafo_H_mult(firstH, Hlst.op(0));
2638 for (std::size_t i=1; i<Hlst.nops(); i++) {
2639 result *= Hlst.op(i);
2641 result = result.expand();
2642 map_trafo_H_mult recursion;
2643 return recursion(result);
2656 ex trafo_H_1tx_prepend_zero(
const ex& e,
const ex& arg)
2660 if (is_a<function>(e)) {
2661 name = ex_to<function>(e).get_name();
2666 for (std::size_t i=0; i<e.nops(); i++) {
2667 if (is_a<function>(e.op(i))) {
2668 std::string name = ex_to<function>(e.op(i)).get_name();
2676 lst newparameter = ex_to<lst>(h.op(0));
2677 newparameter.prepend(0);
2678 ex addzeta = convert_H_to_zeta(newparameter);
2679 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
2681 return e * (-H(
lst{ex(0)},1/arg).hold());
2688 ex trafo_H_prepend_one(
const ex& e,
const ex& arg)
2692 if (is_a<function>(e)) {
2693 name = ex_to<function>(e).get_name();
2698 for (std::size_t i=0; i<e.nops(); i++) {
2699 if (is_a<function>(e.op(i))) {
2700 std::string name = ex_to<function>(e.op(i)).get_name();
2708 lst newparameter = ex_to<lst>(h.op(0));
2709 newparameter.prepend(1);
2710 return e.subs(h == H(newparameter, h.op(1)).hold());
2712 return e * H(
lst{ex(1)},1-arg).hold();
2719 ex trafo_H_1tx_prepend_minusone(
const ex& e,
const ex& arg)
2723 if (is_a<function>(e)) {
2724 name = ex_to<function>(e).get_name();
2729 for (std::size_t i=0; i<e.nops(); i++) {
2730 if (is_a<function>(e.op(i))) {
2731 std::string name = ex_to<function>(e.op(i)).get_name();
2739 lst newparameter = ex_to<lst>(h.op(0));
2740 newparameter.prepend(-1);
2741 ex addzeta = convert_H_to_zeta(newparameter);
2742 return e.subs(h == (addzeta-H(newparameter, h.op(1)).hold())).expand();
2744 ex addzeta = convert_H_to_zeta(
lst{ex(-1)});
2745 return (e * (addzeta - H(
lst{ex(-1)},1/arg).hold())).expand();
2752 ex trafo_H_1mxt1px_prepend_minusone(
const ex& e,
const ex& arg)
2756 if (is_a<function>(e)) {
2757 name = ex_to<function>(e).get_name();
2762 for (std::size_t i = 0; i < e.nops(); i++) {
2763 if (is_a<function>(e.op(i))) {
2764 std::string name = ex_to<function>(e.op(i)).get_name();
2772 lst newparameter = ex_to<lst>(h.op(0));
2773 newparameter.prepend(-1);
2774 return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
2776 return (e * H(
lst{ex(-1)},(1-arg)/(1+arg)).hold()).
expand();
2783 ex trafo_H_1mxt1px_prepend_one(
const ex& e,
const ex& arg)
2787 if (is_a<function>(e)) {
2788 name = ex_to<function>(e).get_name();
2793 for (std::size_t i = 0; i < e.nops(); i++) {
2794 if (is_a<function>(e.op(i))) {
2795 std::string name = ex_to<function>(e.op(i)).get_name();
2803 lst newparameter = ex_to<lst>(h.op(0));
2804 newparameter.prepend(1);
2805 return e.subs(h == H(newparameter, h.op(1)).hold()).expand();
2807 return (e * H(
lst{ex(1)},(1-arg)/(1+arg)).hold()).
expand();
2813 struct map_trafo_H_1mx :
public map_function
2815 ex operator()(
const ex& e)
override
2817 if (is_a<add>(e) || is_a<mul>(e)) {
2818 return e.
map(*
this);
2821 if (is_a<function>(e)) {
2822 std::string name = ex_to<function>(e).get_name();
2825 lst parameter = ex_to<lst>(e.op(0));
2829 bool allthesame =
true;
2830 if (parameter.op(0) == 0) {
2831 for (std::size_t i = 1; i < parameter.nops(); i++) {
2832 if (parameter.op(i) != 0) {
2839 for (
int i=parameter.nops(); i>0; i--) {
2840 newparameter.append(1);
2842 return pow(-1, parameter.nops()) * H(newparameter, 1-arg).
hold();
2844 }
else if (parameter.op(0) == -1) {
2845 throw std::runtime_error(
"map_trafo_H_1mx: cannot handle weights equal -1!");
2847 for (std::size_t i = 1; i < parameter.nops(); i++) {
2848 if (parameter.op(i) != 1) {
2855 for (
int i=parameter.nops(); i>0; i--) {
2856 newparameter.append(0);
2858 return pow(-1, parameter.nops()) * H(newparameter, 1-arg).
hold();
2862 lst newparameter = parameter;
2863 newparameter.remove_first();
2865 if (parameter.op(0) == 0) {
2868 ex res = convert_H_to_zeta(parameter);
2870 map_trafo_H_1mx recursion;
2871 ex buffer = recursion(H(newparameter, arg).hold());
2872 if (is_a<add>(buffer)) {
2873 for (std::size_t i = 0; i < buffer.nops(); i++) {
2874 res -= trafo_H_prepend_one(buffer.op(i), arg);
2877 res -= trafo_H_prepend_one(buffer, arg);
2884 map_trafo_H_1mx recursion;
2885 map_trafo_H_mult unify;
2886 ex res = H(
lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
2887 std::size_t firstzero = 0;
2888 while (parameter.op(firstzero) == 1) {
2891 for (std::size_t i = firstzero-1; i < parameter.nops()-1; i++) {
2895 newparameter.append(parameter[j+1]);
2897 newparameter.append(1);
2898 for (; j<parameter.nops()-1; j++) {
2899 newparameter.append(parameter[j+1]);
2901 res -= H(newparameter, arg).hold();
2903 res = recursion(res).expand() / firstzero;
2914 struct map_trafo_H_1overx :
public map_function
2916 ex operator()(
const ex& e)
override
2918 if (is_a<add>(e) || is_a<mul>(e)) {
2919 return e.map(*
this);
2922 if (is_a<function>(e)) {
2923 std::string name = ex_to<function>(e).get_name();
2926 lst parameter = ex_to<lst>(e.op(0));
2930 bool allthesame =
true;
2931 if (parameter.op(0) == 0) {
2932 for (std::size_t i = 1; i < parameter.nops(); i++) {
2933 if (parameter.op(i) != 0) {
2939 return pow(-1, parameter.nops()) * H(parameter, 1/arg).
hold();
2941 }
else if (parameter.op(0) == -1) {
2942 for (std::size_t i = 1; i < parameter.nops(); i++) {
2943 if (parameter.op(i) != -1) {
2949 map_trafo_H_mult unify;
2950 return unify((
pow(H(
lst{ex(-1)},1/arg).hold() - H(
lst{ex(0)},1/arg).hold(), parameter.nops())
2951 /
factorial(parameter.nops())).expand());
2954 for (std::size_t i = 1; i < parameter.nops(); i++) {
2955 if (parameter.op(i) != 1) {
2961 map_trafo_H_mult unify;
2962 return unify((
pow(H(
lst{ex(1)},1/arg).hold() + H(
lst{ex(0)},1/arg).hold() + H_polesign, parameter.nops())
2963 /
factorial(parameter.nops())).expand());
2967 lst newparameter = parameter;
2968 newparameter.remove_first();
2970 if (parameter.op(0) == 0) {
2973 ex res = convert_H_to_zeta(parameter);
2974 map_trafo_H_1overx recursion;
2975 ex buffer = recursion(H(newparameter, arg).hold());
2976 if (is_a<add>(buffer)) {
2977 for (std::size_t i = 0; i < buffer.nops(); i++) {
2978 res += trafo_H_1tx_prepend_zero(buffer.op(i), arg);
2981 res += trafo_H_1tx_prepend_zero(buffer, arg);
2985 }
else if (parameter.op(0) == -1) {
2988 ex res = convert_H_to_zeta(parameter);
2989 map_trafo_H_1overx recursion;
2990 ex buffer = recursion(H(newparameter, arg).hold());
2991 if (is_a<add>(buffer)) {
2992 for (std::size_t i = 0; i < buffer.nops(); i++) {
2993 res += trafo_H_1tx_prepend_zero(buffer.op(i), arg) - trafo_H_1tx_prepend_minusone(buffer.op(i), arg);
2996 res += trafo_H_1tx_prepend_zero(buffer, arg) - trafo_H_1tx_prepend_minusone(buffer, arg);
3003 map_trafo_H_1overx recursion;
3004 map_trafo_H_mult unify;
3005 ex res = H(
lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
3006 std::size_t firstzero = 0;
3007 while (parameter.op(firstzero) == 1) {
3010 for (std::size_t i = firstzero-1; i < parameter.nops() - 1; i++) {
3014 newparameter.append(parameter[j+1]);
3016 newparameter.append(1);
3017 for (; j<parameter.nops()-1; j++) {
3018 newparameter.append(parameter[j+1]);
3020 res -= H(newparameter, arg).hold();
3022 res = recursion(res).expand() / firstzero;
3035 struct map_trafo_H_1mxt1px :
public map_function
3037 ex operator()(
const ex& e)
override
3039 if (is_a<add>(e) || is_a<mul>(e)) {
3040 return e.map(*
this);
3043 if (is_a<function>(e)) {
3044 std::string name = ex_to<function>(e).get_name();
3047 lst parameter = ex_to<lst>(e.op(0));
3051 bool allthesame =
true;
3052 if (parameter.op(0) == 0) {
3053 for (std::size_t i = 1; i < parameter.nops(); i++) {
3054 if (parameter.op(i) != 0) {
3060 map_trafo_H_mult unify;
3061 return unify((
pow(-H(
lst{ex(1)},(1-arg)/(1+arg)).hold() - H(
lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.nops())
3062 /
factorial(parameter.nops())).expand());
3064 }
else if (parameter.op(0) == -1) {
3065 for (std::size_t i = 1; i < parameter.nops(); i++) {
3066 if (parameter.op(i) != -1) {
3072 map_trafo_H_mult unify;
3073 return unify((
pow(
log(2) - H(
lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.
nops())
3074 /
factorial(parameter.nops())).expand());
3077 for (std::size_t i = 1; i < parameter.nops(); i++) {
3078 if (parameter.op(i) != 1) {
3084 map_trafo_H_mult unify;
3085 return unify((
pow(-
log(2) - H(
lst{ex(0)},(1-arg)/(1+arg)).hold() + H(
lst{ex(-1)},(1-arg)/(1+arg)).hold(), parameter.
nops())
3086 /
factorial(parameter.nops())).expand());
3090 lst newparameter = parameter;
3091 newparameter.remove_first();
3093 if (parameter.op(0) == 0) {
3096 ex res = convert_H_to_zeta(parameter);
3097 map_trafo_H_1mxt1px recursion;
3098 ex buffer = recursion(H(newparameter, arg).hold());
3099 if (is_a<add>(buffer)) {
3100 for (std::size_t i = 0; i < buffer.nops(); i++) {
3101 res -= trafo_H_1mxt1px_prepend_one(buffer.op(i), arg) + trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
3104 res -= trafo_H_1mxt1px_prepend_one(buffer, arg) + trafo_H_1mxt1px_prepend_minusone(buffer, arg);
3108 }
else if (parameter.op(0) == -1) {
3111 ex res = convert_H_to_zeta(parameter);
3112 map_trafo_H_1mxt1px recursion;
3113 ex buffer = recursion(H(newparameter, arg).hold());
3114 if (is_a<add>(buffer)) {
3115 for (std::size_t i = 0; i < buffer.nops(); i++) {
3116 res -= trafo_H_1mxt1px_prepend_minusone(buffer.op(i), arg);
3119 res -= trafo_H_1mxt1px_prepend_minusone(buffer, arg);
3126 map_trafo_H_1mxt1px recursion;
3127 map_trafo_H_mult unify;
3128 ex res = H(
lst{ex(1)}, arg).hold() * H(newparameter, arg).hold();
3129 std::size_t firstzero = 0;
3130 while (parameter.op(firstzero) == 1) {
3133 for (std::size_t i = firstzero - 1; i < parameter.nops() - 1; i++) {
3137 newparameter.append(parameter[j+1]);
3139 newparameter.append(1);
3140 for (; j<parameter.nops()-1; j++) {
3141 newparameter.append(parameter[j+1]);
3143 res -= H(newparameter, arg).hold();
3145 res = recursion(res).expand() / firstzero;
3158 cln::cl_N H_do_sum(
const std::vector<int>&
m,
const cln::cl_N&
x)
3160 const int j =
m.size();
3162 std::vector<cln::cl_N> t(j);
3164 cln::cl_F
one = cln::cl_float(1, cln::float_format(
Digits));
3171 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),
m[j-1]);
3172 for (
int k=j-2;
k>=1;
k--) {
3173 t[
k] = t[
k] + t[
k+1] / cln::expt(cln::cl_I(q+j-1-
k),
m[
k]);
3175 t[0] = t[0] + t[1] *
factor / cln::expt(cln::cl_I(q+j-1),
m[0]);
3177 }
while (t[0] != t0buf);
3197 if (is_a<lst>(x1)) {
3200 if (is_a<numeric>(x2)) {
3201 x = ex_to<numeric>(x2).to_cl_N();
3204 if (is_a<numeric>(x2_val)) {
3205 x = ex_to<numeric>(x2_val).to_cl_N();
3209 for (std::size_t i = 0; i < x1.
nops(); i++) {
3211 return H(x1, x2).hold();
3214 if (x1.
nops() < 1) {
3215 return H(x1, x2).hold();
3218 const lst& morg = ex_to<lst>(x1);
3220 if (*(--morg.
end()) == 0) {
3222 map_trafo_H_reduce_trailing_zeros filter;
3223 return filter(H(x1, xtemp).hold()).subs(xtemp==x2).evalf();
3227 for (
const auto & it : morg) {
3229 for (
ex count=it-1; count > 0; count--) {
3233 }
else if (it <= -1) {
3234 for (
ex count=it+1; count < 0; count++) {
3248 if (convert_parameter_H_to_Li(
m, m_lst, s_lst, pf)) {
3250 std::vector<int> m_int;
3251 std::vector<cln::cl_N> x_cln;
3252 for (
auto it_int = m_lst.
begin(), it_cln = s_lst.
begin();
3253 it_int != m_lst.
end(); it_int++, it_cln++) {
3254 m_int.push_back(ex_to<numeric>(*it_int).to_int());
3255 x_cln.push_back(ex_to<numeric>(*it_cln).to_cl_N());
3257 x_cln.front() = x_cln.front() *
x;
3258 return pf *
numeric(multipleLi_do_sum(m_int, x_cln));
3262 if (m_lst.
nops() == 1) {
3263 return Li(m_lst.
op(0), x2).evalf();
3265 std::vector<int> m_int;
3266 for (
const auto & it : m_lst) {
3267 m_int.push_back(ex_to<numeric>(it).
to_int());
3269 return numeric(H_do_sum(m_int,
x));
3277 if (cln::realpart(
x) < 0) {
3279 for (std::size_t i = 0; i <
m.nops(); i++) {
3281 m.let_op(i) = -
m.op(i);
3289 map_trafo_H_1overx trafo;
3290 res *= trafo(H(
m, xtemp).hold());
3291 if (cln::imagpart(
x) <= 0) {
3292 res = res.
subs(H_polesign == -
I*
Pi);
3294 res = res.
subs(H_polesign ==
I*
Pi);
3300 bool has_minus_one =
false;
3301 for (
const auto & it :
m) {
3303 has_minus_one =
true;
3311 map_trafo_H_1mxt1px trafo;
3312 res *= trafo(H(
m, xtemp).hold());
3315 if (has_minus_one) {
3316 map_trafo_H_convert_to_Li filter;
3317 return filter(H(
m,
numeric(
x)).hold()).evalf();
3319 map_trafo_H_1mx trafo;
3320 res *= trafo(H(
m, xtemp).hold());
3326 return H(x1,x2).hold();
3333 if (is_a<lst>(m_)) {
3338 if (
m.nops() == 0) {
3346 if (*
m.begin() >
_ex1) {
3352 }
else if (*
m.begin() <
_ex_1) {
3358 }
else if (*
m.begin() ==
_ex0) {
3365 for (
auto it = ++
m.begin(); it !=
m.end(); it++) {
3377 }
else if (*it <
_ex_1) {
3397 }
else if (
step == 1) {
3410 return H(m_,
x).hold();
3414 return convert_H_to_zeta(
m);
3420 return H(m_,
x).hold();
3427 }
else if ((
step == 1) && (pos1 ==
_ex0)){
3432 return pow(-1, p) * S(
n, p, -
x);
3441 return H(m_,
x).hold();
3448 return pseries(rel, std::move(seq));
3455 if (deriv_param == 0) {
3459 if (is_a<lst>(m_)) {
3475 return 1/(1-
x) * H(
m,
x);
3476 }
else if (mb ==
_ex_1) {
3477 return 1/(1+
x) * H(
m,
x);
3487 if (is_a<lst>(m_)) {
3492 c.s <<
"\\mathrm{H}_{";
3493 auto itm =
m.begin();
3496 for (; itm !=
m.end(); itm++) {
3512 do_not_evalf_params());
3518 map_trafo_H_reduce_trailing_zeros filter;
3519 map_trafo_H_convert_to_Li filter2;
3521 return filter2(filter(H(
m,
x).hold()));
3523 return filter2(filter(H(
lst{
m},
x).hold()));
3542 const cln::cl_N lambda = cln::cl_N(
"319/320");
3544 void halfcyclic_convolute(
const std::vector<cln::cl_N>& a,
const std::vector<cln::cl_N>& b, std::vector<cln::cl_N>&
c)
3546 const int size = a.size();
3547 for (
int n=0;
n<size;
n++) {
3549 for (
int m=0;
m<=
n;
m++) {
3557 static void initcX(std::vector<cln::cl_N>& crX,
3558 const std::vector<int>& s,
3561 std::vector<cln::cl_N> crB(L2 + 1);
3562 for (
int i=0; i<=L2; i++)
3567 std::vector<std::vector<cln::cl_N>> crG(s.size() - 1, std::vector<cln::cl_N>(L2 + 1));
3568 for (
int m=0;
m < (int)s.size() - 1;
m++) {
3571 for (
int i = 0; i <= L2; i++)
3577 for (std::size_t
m = 0;
m < s.size() - 1;
m++) {
3578 std::vector<cln::cl_N> Xbuf(L2 + 1);
3579 for (
int i = 0; i <= L2; i++)
3580 Xbuf[i] = crX[i] * crG[
m][i];
3582 halfcyclic_convolute(Xbuf, crB, crX);
3588 static cln::cl_N crandall_Y_loop(
const cln::cl_N& Sqk,
3589 const std::vector<cln::cl_N>& crX)
3591 cln::cl_F
one = cln::cl_float(1, cln::float_format(
Digits));
3592 cln::cl_N
factor = cln::expt(lambda, Sqk);
3593 cln::cl_N res =
factor / Sqk * crX[0] *
one;
3600 res = res + crX[N] *
factor / (N+Sqk);
3601 }
while (((res != resbuf) || cln::zerop(crX[N])) && (N+1 < crX.size()));
3607 static void calc_f(std::vector<std::vector<cln::cl_N>>& f_kj,
3608 const int maxr,
const int L1)
3610 cln::cl_N t0, t1, t2, t3, t4;
3612 auto it = f_kj.begin();
3613 cln::cl_F
one = cln::cl_float(1, cln::float_format(
Digits));
3617 for (
k=1;
k<=L1;
k++) {
3620 for (j=1; j<=maxr; j++) {
3623 for (i=2; i<=j; i++) {
3627 (*it).push_back(t2 * t3 * cln::expt(cln::cl_I(
k),-j) *
one);
3635 static cln::cl_N crandall_Z(
const std::vector<int>& s,
3636 const std::vector<std::vector<cln::cl_N>>& f_kj)
3638 const int j = s.size();
3647 t0 = t0 + f_kj[q+j-2][s[0]-1];
3648 }
while ((t0 != t0buf) && (q+j-1 < f_kj.size()));
3653 std::vector<cln::cl_N> t(j);
3660 t[j-1] = t[j-1] + 1 / cln::expt(cln::cl_I(q),s[j-1]);
3661 for (
int k=j-2;
k>=1;
k--) {
3662 t[
k] = t[
k] + t[
k+1] / cln::expt(cln::cl_I(q+j-1-
k), s[
k]);
3664 t[0] = t[0] + t[1] * f_kj[q+j-2][s[0]-1];
3665 }
while ((t[0] != t0buf) && (q+j-1 < f_kj.size()));
3672 cln::cl_N zeta_do_sum_Crandall(
const std::vector<int>& s)
3674 std::vector<int>
r = s;
3675 const int j =
r.size();
3690 }
else if (
Digits < 86) {
3692 }
else if (
Digits < 192) {
3694 }
else if (
Digits < 394) {
3696 }
else if (
Digits < 808) {
3698 }
else if (
Digits < 1636) {
3709 for (
int i=0; i<j; i++) {
3716 std::vector<std::vector<cln::cl_N>> f_kj(L1);
3717 calc_f(f_kj, maxr, L1);
3721 std::vector<int> rz;
3724 for (
int k=
r.size()-1;
k>0;
k--) {
3726 rz.insert(rz.begin(),
r.back());
3727 skp1buf = rz.front();
3731 std::vector<cln::cl_N> crX;
3734 for (
int q=0; q<skp1buf; q++) {
3736 cln::cl_N pp1 = crandall_Y_loop(Srun+q-
k, crX);
3737 cln::cl_N pp2 = crandall_Z(rz, f_kj);
3747 rz.front() = skp1buf;
3749 rz.insert(rz.begin(),
r.back());
3751 std::vector<cln::cl_N> crX;
3752 initcX(crX, rz, L2);
3754 res = (res + crandall_Y_loop(S-j, crX)) / r0factorial
3755 + crandall_Z(rz, f_kj);
3761 cln::cl_N zeta_do_sum_simple(
const std::vector<int>&
r)
3763 const int j =
r.size();
3766 std::vector<cln::cl_N> t(j);
3767 cln::cl_F
one = cln::cl_float(1, cln::float_format(
Digits));
3774 t[j-1] = t[j-1] +
one / cln::expt(cln::cl_I(q),
r[j-1]);
3775 for (
int k=j-2;
k>=0;
k--) {
3776 t[
k] = t[
k] +
one * t[
k+1] / cln::expt(cln::cl_I(q+j-1-
k),
r[
k]);
3778 }
while (t[0] != t0buf);
3785 cln::cl_N zeta_do_Hoelder_convolution(
const std::vector<int>& m_,
const std::vector<int>& s_)
3789 std::vector<int> s = s_;
3790 std::vector<int> m_p = m_;
3791 std::vector<int> m_q;
3793 std::vector<cln::cl_N> s_p(s.size(), cln::cl_N(1));
3794 s_p[0] = s_p[0] * cln::cl_N(
"1/2");
3797 for (std::size_t i = 0; i < s_.size(); i++) {
3804 std::vector<cln::cl_N> s_q;
3805 cln::cl_N signum = 1;
3808 cln::cl_N res = multipleLi_do_sum(m_p, s_p);
3814 if (s.front() > 0) {
3815 if (m_p.front() == 1) {
3816 m_p.erase(m_p.begin());
3817 s_p.erase(s_p.begin());
3818 if (s_p.size() > 0) {
3819 s_p.front() = s_p.front() * cln::cl_N(
"1/2");
3825 m_q.insert(m_q.begin(), 1);
3826 if (s_q.size() > 0) {
3827 s_q.front() = s_q.front() * 2;
3829 s_q.insert(s_q.begin(), cln::cl_N(
"1/2"));
3832 if (m_p.front() == 1) {
3833 m_p.erase(m_p.begin());
3834 cln::cl_N spbuf = s_p.front();
3835 s_p.erase(s_p.begin());
3836 if (s_p.size() > 0) {
3837 s_p.front() = s_p.front() * spbuf;
3840 m_q.insert(m_q.begin(), 1);
3841 if (s_q.size() > 0) {
3842 s_q.front() = s_q.front() * 4;
3844 s_q.insert(s_q.begin(), cln::cl_N(
"1/4"));
3848 m_q.insert(m_q.begin(), 1);
3849 if (s_q.size() > 0) {
3850 s_q.front() = s_q.front() * 2;
3852 s_q.insert(s_q.begin(), cln::cl_N(
"1/2"));
3857 if (m_p.size() == 0)
break;
3859 res = res + signum * multipleLi_do_sum(m_p, s_p) * multipleLi_do_sum(m_q, s_q);
3864 res = res + signum * multipleLi_do_sum(m_q, s_q);
3884 if (is_exactly_a<lst>(
x) && (
x.
nops()>1)) {
3887 const int count =
x.
nops();
3888 const lst& xlst = ex_to<lst>(
x);
3889 std::vector<int>
r(count);
3890 std::vector<int> si(count);
3893 auto it1 = xlst.
begin();
3894 auto it2 =
r.begin();
3895 auto it_swrite = si.begin();
3900 *it2 = ex_to<numeric>(*it1).to_int();
3905 }
while (it2 !=
r.end());
3914 return numeric(zeta_do_Hoelder_convolution(
r, si));
3918 int limit = (
Digits>17) ? 10 : 6;
3919 if ((
r[0] < limit) || ((count > 3) && (
r[1] < limit/2))) {
3920 return numeric(zeta_do_sum_Crandall(
r));
3922 return numeric(zeta_do_sum_simple(
r));
3927 if (is_exactly_a<numeric>(
x) && (
x != 1)) {
3929 return zeta(ex_to<numeric>(
x));
3930 }
catch (
const dunno &e) { }
3939 if (is_exactly_a<lst>(
m)) {
3940 if (
m.nops() == 1) {
3941 return zeta(
m.op(0));
3947 const numeric& y = ex_to<numeric>(
m);
3983 if (is_exactly_a<lst>(
m)) {
3986 return zetaderiv(
_ex1,
m);
3994 if (is_a<lst>(m_)) {
3995 const lst&
m = ex_to<lst>(m_);
3996 auto it =
m.begin();
3999 for (; it !=
m.end(); it++) {
4015 do_not_evalf_params().
4030 if (is_exactly_a<lst>(
x)) {
4033 const int count =
x.
nops();
4034 const lst& xlst = ex_to<lst>(
x);
4035 const lst& slst = ex_to<lst>(s);
4036 std::vector<int> xi(count);
4037 std::vector<int> si(count);
4040 auto it_xread = xlst.
begin();
4041 auto it_sread = slst.
begin();
4042 auto it_xwrite = xi.begin();
4043 auto it_swrite = si.begin();
4048 *it_xwrite = ex_to<numeric>(*it_xread).to_int();
4049 if (*it_sread > 0) {
4058 }
while (it_xwrite != xi.end());
4061 if ((xi[0] == 1) && (si[0] == 1)) {
4066 return numeric(zeta_do_Hoelder_convolution(xi, si));
4076 if (is_exactly_a<lst>(s_)) {
4077 const lst& s = ex_to<lst>(s_);
4078 for (
const auto & it : s) {
4097 if (is_exactly_a<lst>(
m)) {
4101 return zetaderiv(
_ex1,
m);
4111 if (is_a<lst>(m_)) {
4117 if (is_a<lst>(s_)) {
4123 auto itm =
m.begin();
4124 auto its = s.
begin();
4126 c.s <<
"\\overline{";
4134 for (; itm !=
m.end(); itm++, its++) {
4137 c.s <<
"\\overline{";
4153 do_not_evalf_params().