elliptic multiple polylogarithms or iterated integrals of modular forms.
Changes to be committed:
modified: check/CMakeLists.txt
modified: check/Makefile.am
new file: check/exam_inifcns_elliptic.cpp
modified: doc/tutorial/ginac.texi
modified: ginac/CMakeLists.txt
modified: ginac/Makefile.am
modified: ginac/ginac.h
modified: ginac/inifcns.h
new file: ginac/inifcns_elliptic.cpp
new file: ginac/integration_kernel.cpp
new file: ginac/integration_kernel.h
new file: ginac/utils_multi_iterator.h
modified: ginsh/ginsh_parser.ypp
exam_powerlaws
exam_inifcns
exam_inifcns_nstdsums
+ exam_inifcns_elliptic
exam_differentiation
exam_polygcd
exam_collect_common_factors
exam_powerlaws \
exam_inifcns \
exam_inifcns_nstdsums \
+ exam_inifcns_elliptic \
exam_differentiation \
exam_polygcd \
exam_collect_common_factors \
exam_inifcns_nstdsums.h
exam_inifcns_nstdsums_LDADD = ../ginac/libginac.la
+exam_inifcns_elliptic_SOURCES = exam_inifcns_elliptic.cpp
+exam_inifcns_elliptic_LDADD = ../ginac/libginac.la
+
exam_differentiation_SOURCES = exam_differentiation.cpp
exam_differentiation_LDADD = ../ginac/libginac.la
--- /dev/null
+/** @file exam_inifcns_nstdsums.cpp
+ *
+ * This test routine applies assorted tests on initially known higher level
+ * functions. */
+
+/*
+ * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "ginac.h"
+using namespace GiNaC;
+
+#include <iostream>
+#include <fstream>
+using namespace std;
+
+
+static unsigned check_q_expansion()
+{
+ unsigned err = 0;
+
+ symbol q("q");
+ int order = 200;
+
+ ex res;
+
+ // q-expansions from Sage up to order 200
+ // notation Ek_N_|a|_|b|_K
+ ex E1_12_1_4_1 = pow(q,2)+pow(q,4)+2*pow(q,5)+pow(q,8)+pow(q,9)+2*pow(q,10)+2*pow(q,13)+pow(q,16)+2*pow(q,17)+pow(q,18)+2*pow(q,20)+3*pow(q,25)+2*pow(q,26)
+ +2*pow(q,29)+pow(q,32)+2*pow(q,34)+pow(q,36)+2*pow(q,37)+2*pow(q,40)+2*pow(q,41)+2*pow(q,45)+pow(q,49)+3*pow(q,50)+2*pow(q,52)+2*pow(q,53)+2*pow(q,58)
+ +2*pow(q,61)+pow(q,64)+4*pow(q,65)+2*pow(q,68)+pow(q,72)+2*pow(q,73)+2*pow(q,74)+2*pow(q,80)+pow(q,81)+2*pow(q,82)+4*pow(q,85)+2*pow(q,89)+2*pow(q,90)
+ +2*pow(q,97)+pow(q,98)+3*pow(q,100)+2*pow(q,101)+2*pow(q,104)+2*pow(q,106)+2*pow(q,109)+2*pow(q,113)+2*pow(q,116)+2*pow(q,117)+pow(q,121)+2*pow(q,122)
+ +4*pow(q,125)+pow(q,128)+4*pow(q,130)+2*pow(q,136)+2*pow(q,137)+pow(q,144)+4*pow(q,145)+2*pow(q,146)+2*pow(q,148)+2*pow(q,149)+2*pow(q,153)+2*pow(q,157)
+ +2*pow(q,160)+pow(q,162)+2*pow(q,164)+3*pow(q,169)+4*pow(q,170)+2*pow(q,173)+2*pow(q,178)+2*pow(q,180)+2*pow(q,181)+4*pow(q,185)+2*pow(q,193)
+ +2*pow(q,194)+pow(q,196)+2*pow(q,197)+numeric(1,4)+q;
+
+ ex E1_12_1_4_3 = numeric(1,4)+pow(q,3)+pow(q,6)+pow(q,12)+2*pow(q,15)+pow(q,24)+pow(q,27)+2*pow(q,30)+2*pow(q,39)+pow(q,48)+2*pow(q,51)+pow(q,54)+2*pow(q,60)
+ +3*pow(q,75)+2*pow(q,78)+2*pow(q,87)+pow(q,96)+2*pow(q,102)+pow(q,108)+2*pow(q,111)+2*pow(q,120)+2*pow(q,123)+2*pow(q,135)+pow(q,147)+3*pow(q,150)
+ +2*pow(q,156)+2*pow(q,159)+2*pow(q,174)+2*pow(q,183)+pow(q,192)+4*pow(q,195);
+
+ ex E1_12_1_3_1 = 2*pow(q,175)+2*pow(q,189)+2*pow(q,199)+numeric(1,6)+2*pow(q,129)+4*pow(q,133)+2*pow(q,139)+2*pow(q,151)+2*pow(q,163)+2*pow(q,171)+2*pow(q,172)
+ +2*pow(q,31)+2*pow(q,43)+2*pow(q,57)+2*pow(q,63)+2*pow(q,67)+2*pow(q,76)+2*pow(q,79)+2*pow(q,84)+4*pow(q,91)+2*pow(q,93)+2*pow(q,103)+2*pow(q,112)
+ +2*pow(q,124)+2*pow(q,127)+pow(q,4)+pow(q,9)+2*pow(q,13)+pow(q,16)+pow(q,25)+pow(q,36)+2*pow(q,37)+3*pow(q,49)+2*pow(q,52)+2*pow(q,61)+pow(q,64)
+ +2*pow(q,73)+pow(q,81)+2*pow(q,97)+pow(q,100)+2*pow(q,109)+2*pow(q,117)+pow(q,121)+pow(q,144)+2*pow(q,148)+2*pow(q,157)+3*pow(q,169)+2*pow(q,181)
+ +2*pow(q,193)+3*pow(q,196)+pow(q,3)+pow(q,12)+pow(q,27)+2*pow(q,39)+pow(q,48)+pow(q,75)+pow(q,108)+2*pow(q,111)+3*pow(q,147)+2*pow(q,156)+2*pow(q,183)
+ +pow(q,192)+2*pow(q,7)+2*pow(q,19)+2*pow(q,21)+2*pow(q,28)+q;
+
+ ex E1_12_1_3_2 = numeric(1,6)+pow(q,2)+pow(q,8)+pow(q,18)+2*pow(q,26)+pow(q,32)+pow(q,50)+pow(q,72)+2*pow(q,74)+3*pow(q,98)+2*pow(q,104)+2*pow(q,122)+pow(q,128)
+ +2*pow(q,146)+pow(q,162)+2*pow(q,194)+pow(q,6)+pow(q,24)+pow(q,54)+2*pow(q,78)+pow(q,96)+pow(q,150)+2*pow(q,14)+2*pow(q,38)+2*pow(q,42)+2*pow(q,56)
+ +2*pow(q,62)+2*pow(q,86)+2*pow(q,114)+2*pow(q,126)+2*pow(q,134)+2*pow(q,152)+2*pow(q,158)+2*pow(q,168)+4*pow(q,182)+2*pow(q,186);
+
+ ex E1_12_1_3_4 = numeric(1,6)+pow(q,4)+pow(q,12)+pow(q,16)+2*pow(q,28)+pow(q,36)+pow(q,48)+2*pow(q,52)+pow(q,64)+2*pow(q,76)+2*pow(q,84)+pow(q,100)+pow(q,108)
+ +2*pow(q,112)+2*pow(q,124)+pow(q,144)+2*pow(q,148)+2*pow(q,156)+2*pow(q,172)+pow(q,192)+3*pow(q,196);
+
+
+ ex E2_12_1_1_2 = 248*pow(q,175)+320*pow(q,189)+200*pow(q,199)+176*pow(q,129)+160*pow(q,133)+140*pow(q,139)+152*pow(q,151)+164*pow(q,163)+260*pow(q,171)
+ +44*pow(q,172)+32*pow(q,31)+44*pow(q,43)+80*pow(q,57)+104*pow(q,63)+68*pow(q,67)+20*pow(q,76)+80*pow(q,79)+32*pow(q,84)+112*pow(q,91)+128*pow(q,93)
+ +104*pow(q,103)+8*pow(q,112)+32*pow(q,124)+128*pow(q,127)+pow(q,2)+pow(q,4)+6*pow(q,5)+pow(q,8)+13*pow(q,9)+6*pow(q,10)+14*pow(q,13)+pow(q,16)
+ +18*pow(q,17)+13*pow(q,18)+6*pow(q,20)+31*pow(q,25)+14*pow(q,26)+30*pow(q,29)+pow(q,32)+18*pow(q,34)+13*pow(q,36)+38*pow(q,37)+6*pow(q,40)+42*pow(q,41)
+ +78*pow(q,45)+57*pow(q,49)+31*pow(q,50)+14*pow(q,52)+54*pow(q,53)+30*pow(q,58)+62*pow(q,61)+pow(q,64)+84*pow(q,65)+18*pow(q,68)+13*pow(q,72)+74*pow(q,73)
+ +38*pow(q,74)+6*pow(q,80)+121*pow(q,81)+42*pow(q,82)+108*pow(q,85)+90*pow(q,89)+78*pow(q,90)+98*pow(q,97)+57*pow(q,98)+31*pow(q,100)+102*pow(q,101)
+ +14*pow(q,104)+54*pow(q,106)+110*pow(q,109)+114*pow(q,113)+30*pow(q,116)+182*pow(q,117)+133*pow(q,121)+62*pow(q,122)+156*pow(q,125)+pow(q,128)+84*pow(q,130)
+ +18*pow(q,136)+138*pow(q,137)+13*pow(q,144)+180*pow(q,145)+74*pow(q,146)+38*pow(q,148)+150*pow(q,149)+234*pow(q,153)+158*pow(q,157)+6*pow(q,160)+121*pow(q,162)
+ +42*pow(q,164)+183*pow(q,169)+108*pow(q,170)+174*pow(q,173)+90*pow(q,178)+78*pow(q,180)+182*pow(q,181)+228*pow(q,185)+194*pow(q,193)+98*pow(q,194)+57*pow(q,196)
+ +198*pow(q,197)+4*pow(q,3)+4*pow(q,6)+4*pow(q,12)+24*pow(q,15)+4*pow(q,24)+40*pow(q,27)+24*pow(q,30)+56*pow(q,39)+4*pow(q,48)+72*pow(q,51)+40*pow(q,54)
+ +24*pow(q,60)+124*pow(q,75)+56*pow(q,78)+120*pow(q,87)+4*pow(q,96)+72*pow(q,102)+40*pow(q,108)+152*pow(q,111)+24*pow(q,120)+168*pow(q,123)+240*pow(q,135)
+ +228*pow(q,147)+124*pow(q,150)+56*pow(q,156)+216*pow(q,159)+120*pow(q,174)+248*pow(q,183)+4*pow(q,192)+336*pow(q,195)+8*pow(q,7)+20*pow(q,19)+32*pow(q,21)
+ +8*pow(q,28)+8*pow(q,14)+20*pow(q,38)+32*pow(q,42)+8*pow(q,56)+32*pow(q,62)+44*pow(q,86)+80*pow(q,114)+104*pow(q,126)+68*pow(q,134)+20*pow(q,152)+80*pow(q,158)
+ +32*pow(q,168)+112*pow(q,182)+128*pow(q,186)+12*pow(q,11)+12*pow(q,22)+24*pow(q,23)+48*pow(q,33)+48*pow(q,35)+12*pow(q,44)+24*pow(q,46)+48*pow(q,47)+72*pow(q,55)
+ +60*pow(q,59)+48*pow(q,66)+96*pow(q,69)+48*pow(q,70)+72*pow(q,71)+96*pow(q,77)+84*pow(q,83)+12*pow(q,88)+24*pow(q,92)+48*pow(q,94)+120*pow(q,95)+156*pow(q,99)
+ +192*pow(q,105)+108*pow(q,107)+72*pow(q,110)+144*pow(q,115)+60*pow(q,118)+144*pow(q,119)+132*pow(q,131)+48*pow(q,132)+96*pow(q,138)+48*pow(q,140)+192*pow(q,141)
+ +72*pow(q,142)+168*pow(q,143)+96*pow(q,154)+192*pow(q,155)+192*pow(q,161)+288*pow(q,165)+84*pow(q,166)+168*pow(q,167)+12*pow(q,176)+240*pow(q,177)+180*pow(q,179)
+ +24*pow(q,184)+216*pow(q,187)+48*pow(q,188)+120*pow(q,190)+192*pow(q,191)+156*pow(q,198)+numeric(1,24)+q;
+
+ ex E3_12_3_1_4 = 1850*pow(q,172)+362*pow(q,76)+450*pow(q,84)+650*pow(q,112)+962*pow(q,124)+pow(q,4)+3*pow(q,8)+13*pow(q,16)+24*pow(q,20)+51*pow(q,32)+81*pow(q,36)
+ +72*pow(q,40)+170*pow(q,52)+205*pow(q,64)+288*pow(q,68)+243*pow(q,72)+312*pow(q,80)+601*pow(q,100)+510*pow(q,104)+840*pow(q,116)+819*pow(q,128)+864*pow(q,136)
+ +1053*pow(q,144)+1370*pow(q,148)+1224*pow(q,160)+1680*pow(q,164)+1944*pow(q,180)+2451*pow(q,196)+9*pow(q,12)+27*pow(q,24)+117*pow(q,48)+216*pow(q,60)
+ +459*pow(q,96)+729*pow(q,108)+648*pow(q,120)+1530*pow(q,156)+1845*pow(q,192)+50*pow(q,28)+150*pow(q,56)+1086*pow(q,152)+1350*pow(q,168)+120*pow(q,44)
+ +360*pow(q,88)+528*pow(q,92)+1080*pow(q,132)+1200*pow(q,140)+1560*pow(q,176)+1584*pow(q,184)+2208*pow(q,188);
+
+ // basis of Eisenstein space for Gamma_1(12) of weight 1
+ res = series_to_poly(Eisenstein_kernel(1, 12, 1, -3, 1).q_expansion_modular_form(q, order)) - E1_12_1_3_1;
+ if ( res != 0 ) err++;
+
+ res = series_to_poly(Eisenstein_kernel(1, 12, 1, -3, 2).q_expansion_modular_form(q, order)) - E1_12_1_3_2;
+ if ( res != 0 ) err++;
+
+ res = series_to_poly(Eisenstein_kernel(1, 12, 1, -3, 4).q_expansion_modular_form(q, order)) - E1_12_1_3_4;
+ if ( res != 0 ) err++;
+
+ res = series_to_poly(Eisenstein_kernel(1, 12, 1, -4, 1).q_expansion_modular_form(q, order)) - E1_12_1_4_1;
+ if ( res != 0 ) err++;
+
+ res = series_to_poly(Eisenstein_kernel(1, 12, 1, -4, 3).q_expansion_modular_form(q, order)) - E1_12_1_4_3;
+ if ( res != 0 ) err++;
+
+ // test one series of weight 2
+ res = series_to_poly(Eisenstein_kernel(2, 12, 1, 1, 2).q_expansion_modular_form(q, order)) - E2_12_1_1_2;
+ if ( res != 0 ) err++;
+
+ // and one of weight 3
+ res = series_to_poly(Eisenstein_kernel(3, 12, -3, 1, 4).q_expansion_modular_form(q, order)) - E3_12_3_1_4;
+ if ( res != 0 ) err++;
+
+ return err;
+}
+
+static unsigned check_polylogs()
+{
+ unsigned err = 0;
+
+ int digitsbuf = Digits;
+ Digits = 100;
+ ex prec = 5 * pow(10, -(ex)Digits);
+
+ ex y = numeric(9,10);
+
+ ex z2 = numeric(2);
+ ex z3 = numeric(3);
+
+ ex L0 = basic_log_kernel();
+ ex omega_2 = multiple_polylog_kernel(z2);
+ ex omega_3 = multiple_polylog_kernel(z3);
+
+ ex res1,res2;
+
+ res1 = G(lst{z2},y).evalf();
+ res2 = iterated_integral(lst{omega_2},y).evalf();
+ if ( abs(res1-res2) > prec ) err++;
+
+ res1 = G(lst{0},y).evalf();
+ res2 = iterated_integral(lst{L0},y).evalf();
+ if ( abs(res1-res2) > prec ) err++;
+
+ res1 = G(lst{z2,0},y).evalf();
+ res2 = iterated_integral(lst{omega_2,L0},y).evalf();
+ if ( abs(res1-res2) > prec ) err++;
+
+ res1 = G(lst{0,0},y).evalf();
+ res2 = iterated_integral(lst{L0,L0},y).evalf();
+ if ( abs(res1-res2) > prec ) err++;
+
+ res1 = G(lst{z2,0,0},y).evalf();
+ res2 = iterated_integral(lst{omega_2,L0,L0},y).evalf();
+ if ( abs(res1-res2) > prec ) err++;
+
+ Digits = digitsbuf;
+
+ return err;
+}
+
+static unsigned check_iterated_integral_modular_form_versus_Kronecker_dtau()
+{
+ unsigned err = 0;
+
+ int digitsbuf = Digits;
+ Digits = 30;
+ ex prec = 5 * pow(10, -(ex)Digits);
+
+ ex tau_6 = I;
+ ex qbar_6 = exp(2*Pi*I*tau_6);
+ ex omega_0 = basic_log_kernel();
+
+ ex eta_1 = Eisenstein_kernel(3, 6, -3, 1, 1);
+ ex eta_2 = Eisenstein_kernel(3, 6, -3, 1, 2);
+ ex omega_3 = modular_form_kernel(3, eta_1-8*eta_2);
+ ex res1 = iterated_integral(lst{omega_0,omega_3},qbar_6).evalf();
+
+ ex C_3 = I/sqrt(numeric(3));
+ ex g3_1 = Kronecker_dtau_kernel(3,numeric(1,3),1,C_3);
+ ex g3_2 = Kronecker_dtau_kernel(3,numeric(1,3),2,C_3);
+ ex expr = iterated_integral(lst{omega_0,g3_1},qbar_6) - 4*iterated_integral(lst{omega_0,g3_2},qbar_6);
+ ex res2 = expr.evalf();
+
+ if ( abs(res1-res2) > prec ) err++;
+
+ Digits = digitsbuf;
+
+ return err;
+}
+
+static unsigned check_modular_trafo()
+{
+ unsigned err = 0;
+
+ int digitsbuf = Digits;
+ Digits = 50;
+ ex prec = 5 * pow(10, -(ex)Digits);
+
+ int N_trunc = 100;
+
+ int a = 0;
+ int b = -1;
+ int c = 1;
+ int d = 0;
+
+ ex tau = numeric(1,10)+numeric(4,5)*I;
+ ex qbar = evalf(exp(2*Pi*I*tau));
+ ex qbar_2 = evalf(exp(2*Pi*I*tau*numeric(1,2)));
+ ex qbar_3 = evalf(exp(2*Pi*I*tau*numeric(1,3)));
+ ex qbar_4 = evalf(exp(2*Pi*I*tau*numeric(1,4)));
+ ex qbar_6 = evalf(exp(2*Pi*I*tau*numeric(1,6)));
+ ex qbar_12 = evalf(exp(2*Pi*I*tau*numeric(1,12)));
+
+ ex tau_prime = (a*tau+b)/(c*tau+d);
+ ex qbar_prime = evalf(exp(2*Pi*I*tau_prime));
+ ex qbar_prime_2 = evalf(exp(2*Pi*I*tau_prime*numeric(1,2)));
+ ex qbar_prime_3 = evalf(exp(2*Pi*I*tau_prime*numeric(1,3)));
+ ex qbar_prime_4 = evalf(exp(2*Pi*I*tau_prime*numeric(1,4)));
+ ex qbar_prime_6 = evalf(exp(2*Pi*I*tau_prime*numeric(1,6)));
+ ex qbar_prime_12 = evalf(exp(2*Pi*I*tau_prime*numeric(1,12)));
+
+ numeric k,N,r,s;
+ ex eta,eta_trafo,res1,res2;
+
+ k = 4;
+ N = 1;
+ eta = Eisenstein_kernel(k, N, 1, 1, 1);
+ res1 = ex_to<Eisenstein_kernel>(eta).get_numerical_value(qbar_prime,N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_kernel>(eta).get_numerical_value(qbar,N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 4;
+ N = 2;
+ r = 0;
+ s = 0;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_2,N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_2,N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 4;
+ N = 2;
+ r = 1;
+ s = 1;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_2,N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_2,N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 4;
+ N = 4;
+ r = 2;
+ s = 2;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_4,N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_4,N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 4;
+ N = 6;
+ r = 3;
+ s = 3;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_6,3*N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_6,3*N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 4;
+ N = 12;
+ r = 6;
+ s = 6;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_12,6*N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_12,6*N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 4;
+ N = 2;
+ r = 1;
+ s = 0;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_2,N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_2,N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 6;
+ N = 6;
+ r = 1;
+ s = 5;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_6,2*N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_6,2*N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 6;
+ N = 12;
+ r = 2;
+ s = 10;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_12,4*N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_12,4*N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 4;
+ N = 3;
+ r = 1;
+ s = 2;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_3,N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_3,N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ k = 1;
+ N = 6;
+ r = 1;
+ s = 5;
+ eta = Eisenstein_h_kernel(k, N, r, s);
+ eta_trafo = Eisenstein_h_kernel(k, N, mod(r*d+s*b,N), mod(r*c+s*a,N));
+ res1 = ex_to<Eisenstein_h_kernel>(eta).get_numerical_value(qbar_prime_6,2*N_trunc);
+ res2 = pow(c*tau+d,k)*ex_to<Eisenstein_h_kernel>(eta_trafo).get_numerical_value(qbar_6,2*N_trunc);
+ if ( abs(res1-res2) > prec ) err++;
+
+ Digits = digitsbuf;
+
+ return err;
+}
+
+static unsigned check_Kronecker_g()
+{
+ unsigned err = 0;
+
+ int digitsbuf = Digits;
+ Digits = 20;
+ ex prec = 5 * pow(10, -(ex)Digits);
+
+ ex tau = numeric(1,10)+numeric(2)*I;
+ ex qbar = evalf(exp(2*Pi*I*tau));
+
+ ex z = numeric(2,100)+numeric(1,10)*I;
+ ex wbar = evalf(exp(2*Pi*I*z));
+
+ ex z_j = numeric(-1,10)*I;
+
+ ex res1,res2,res3;
+
+ res1 = Kronecker_dtau_kernel(0,z).get_numerical_value(qbar);
+ res2 = Kronecker_dz_kernel(1,0,tau).get_numerical_value(z);
+ res3 = Kronecker_dz_kernel(1,z_j,tau).get_numerical_value(z+z_j);
+ if ( abs(res1-res2) > prec ) err++;
+ if ( abs(res1-res3) > prec ) err++;
+
+ res1 = Kronecker_dtau_kernel(1,z).get_numerical_value(qbar);
+ res2 = Kronecker_dz_kernel(2,0,tau).get_numerical_value(z);
+ res3 = Kronecker_dz_kernel(2,z_j,tau).get_numerical_value(z+z_j);
+ if ( abs(res1-res2) > prec ) err++;
+ if ( abs(res1-res3) > prec ) err++;
+
+ res1 = Kronecker_dtau_kernel(2,z).get_numerical_value(qbar);
+ res2 = Kronecker_dz_kernel(3,0,tau).get_numerical_value(z);
+ res3 = Kronecker_dz_kernel(3,z_j,tau).get_numerical_value(z+z_j);
+ if ( abs(res1-res2) > prec ) err++;
+ if ( abs(res1-res3) > prec ) err++;
+
+ res1 = Kronecker_dtau_kernel(3,z).get_numerical_value(qbar);
+ res2 = Kronecker_dz_kernel(4,0,tau).get_numerical_value(z);
+ res3 = Kronecker_dz_kernel(4,z_j,tau).get_numerical_value(z+z_j);
+ if ( abs(res1-res2) > prec ) err++;
+ if ( abs(res1-res3) > prec ) err++;
+
+ res1 = Kronecker_dtau_kernel(4,z).get_numerical_value(qbar);
+ res2 = Kronecker_dz_kernel(5,0,tau).get_numerical_value(z);
+ res3 = Kronecker_dz_kernel(5,z_j,tau).get_numerical_value(z+z_j);
+ if ( abs(res1-res2) > prec ) err++;
+ if ( abs(res1-res3) > prec ) err++;
+
+ res1 = Kronecker_dtau_kernel(5,z).get_numerical_value(qbar);
+ res2 = Kronecker_dz_kernel(6,0,tau).get_numerical_value(z);
+ res3 = Kronecker_dz_kernel(6,z_j,tau).get_numerical_value(z+z_j);
+ if ( abs(res1-res2) > prec ) err++;
+ if ( abs(res1-res3) > prec ) err++;
+
+ Digits = digitsbuf;
+
+ return err;
+}
+
+unsigned exam_inifcns_elliptic(void)
+{
+ unsigned result = 0;
+
+ cout << "examining consistency of iterated integrals" << flush;
+
+ result += check_q_expansion();
+ result += check_polylogs();
+ result += check_iterated_integral_modular_form_versus_Kronecker_dtau();
+ result += check_modular_trafo();
+ result += check_Kronecker_g();
+
+ return result;
+}
+
+int main(int argc, char** argv)
+{
+ return exam_inifcns_elliptic();
+}
* Symmetrization::
* Built-in functions:: List of predefined mathematical functions.
* Multiple polylogarithms::
+* Iterated integrals::
* Complex expressions::
* Solving linear systems of equations::
* Input/output:: Input and output of expressions.
@cindex @code{zeta()}
@item @code{zetaderiv(n, x)}
@tab derivatives of Riemann's zeta function
+@item @code{iterated_integral(a, y)}
+@tab iterated integral
+@cindex @code{iterated_integral()}
+@item @code{iterated_integral(a, y, N)}
+@tab iterated integral with explicit truncation parameter
+@cindex @code{iterated_integral()}
@item @code{tgamma(x)}
@tab gamma function
@cindex @code{tgamma()}
@cindex @code{psi()}
@item @code{psi(n, x)}
@tab derivatives of psi function (polygamma functions)
+@item @code{EllipticK(x)}
+@tab complete elliptic integral of the first kind
+@cindex @code{EllipticK()}
+@item @code{EllipticE(x)}
+@tab complete elliptic integral of the second kind
+@cindex @code{EllipticE()}
@item @code{factorial(n)}
@tab factorial function @math{n!}
@cindex @code{factorial()}
argument. Of course, a user can fine-tune this behavior by sequential
calls of several @code{expand()} methods with desired flags.
-@node Multiple polylogarithms, Complex expressions, Built-in functions, Methods and functions
+@node Multiple polylogarithms, Iterated integrals, Built-in functions, Methods and functions
@c node-name, next, previous, up
@subsection Multiple polylogarithms
@cite{Numerical Evaluation of Multiple Polylogarithms},
J.Vollinga, S.Weinzierl, hep-ph/0410259
-@node Complex expressions, Solving linear systems of equations, Multiple polylogarithms, Methods and functions
+@node Iterated integrals, Complex expressions, Multiple polylogarithms, Methods and functions
+@c node-name, next, previous, up
+@subsection Iterated integrals
+
+Multiple polylogarithms are a particular example of iterated integrals.
+An iterated integral is defined by the function @code{iterated_integral(a,y)}.
+The variable @code{y} gives the upper integration limit for the outermost integration, by convention the lower integration limit is always set to zero.
+The variable @code{a} must be a GiNaC @code{lst} containing sub-classes of @code{integration_kernel} as elements.
+The depth of the iterated integral corresponds to the number of elements of @code{a}.
+The available integrands for iterated integrals are
+(for a more detailed description the user is referred to the publications listed at the end of this section)
+@cartouche
+@multitable @columnfractions .40 .60
+@item @strong{Class} @tab @strong{Description}
+@item @code{integration_kernel()}
+@tab Base class, represents the one-form @math{dy}
+@cindex @code{integration_kernel()}
+@item @code{basic_log_kernel()}
+@tab Logarithmic one-form @math{dy/y}
+@cindex @code{basic_log_kernel()}
+@item @code{multiple_polylog_kernel(z_j)}
+@tab The one-form @math{dy/(y-z_j)}
+@cindex @code{multiple_polylog_kernel()}
+@item @code{ELi_kernel(n, m, x, y)}
+@tab The one form @math{ELi_{n;m}(x;y;q) dq/q}
+@cindex @code{ELi_kernel()}
+@item @code{Ebar_kernel(n, m, x, y)}
+@tab The one form @math{\overline{E}_{n;m}(x;y;q) dq/q}
+@cindex @code{Ebar_kernel()}
+@item @code{Kronecker_dtau_kernel(k, z_j, K, C_k)}
+@tab The one form @math{C_k K (k-1)/(2 \pi i)^k g^{(k)}(z_j,K \tau) dq/q}
+@cindex @code{Kronecker_dtau_kernel()}
+@item @code{Kronecker_dz_kernel(k, z_j, tau, K, C_k)}
+@tab The one form @math{C_k (2 \pi i)^{2-k} g^{(k-1)}(z-z_j,K \tau) dz}
+@cindex @code{Kronecker_dz_kernel()}
+@item @code{Eisenstein_kernel(k, N, a, b, K, C_k)}
+@tab The one form @math{C_k E_{k,N,a,b,K}(\tau) dq/q}
+@cindex @code{Eisenstein_kernel()}
+@item @code{Eisenstein_h_kernel(k, N, r, s, C_k)}
+@tab The one form @math{C_k h_{k,N,r,s}(\tau) dq/q}
+@cindex @code{Eisenstein_h_kernel()}
+@item @code{modular_form_kernel(k, P, C_k)}
+@tab The one form @math{C_k P dq/q}
+@cindex @code{modular_form_kernel()}
+@item @code{user_defined_kernel(f, y)}
+@tab The one form @math{f(y) dy}
+@cindex @code{user_defined_kernel()}
+@end multitable
+@end cartouche
+All parameters are assumed to be such that all integration kernels have a convergent Laurent expansion
+around zero with at most a simple pole at zero.
+The iterated integral may also be called with an optional third parameter
+@code{iterated_integral(a,y,N_trunc)}, in which case the numerical evaluation will truncate the series
+expansion at order @code{N_trunc}.
+
+The classes @code{Eisenstein_kernel()}, @code{Eisenstein_h_kernel()} and @code{modular_form_kernel()}
+provide a method @code{q_expansion_modular_form(q, order)}, which can used to obtain the q-expansion
+of @math{E_{k,N,a,b,K}(\tau)}, @math{h_{k,N,r,s}(\tau)} or @math{P} to the specified order.
+
+Useful publications:
+
+@cite{Numerical evaluation of iterated integrals related to elliptic Feynman integrals},
+M.Walden, S.Weinzierl, arXiv:2010.xxxxx
+
+@node Complex expressions, Solving linear systems of equations, Iterated integrals, Methods and functions
@c node-name, next, previous, up
@section Complex expressions
@c
inifcns.cpp
inifcns_gamma.cpp
inifcns_nstdsums.cpp
+ inifcns_elliptic.cpp
inifcns_trans.cpp
+ integration_kernel.cpp
integral.cpp
lst.cpp
matrix.cpp
indexed.h
inifcns.h
integral.h
+ integration_kernel.h
lst.h
matrix.h
mul.h
utils.h
crc32.h
hash_seed.h
+ utils_multi_iterator.h
parser/lexer.h
parser/debug.h
polynomial/gcd_euclid.h
libginac_la_SOURCES = add.cpp archive.cpp basic.cpp clifford.cpp color.cpp \
constant.cpp ex.cpp excompiler.cpp expair.cpp expairseq.cpp exprseq.cpp \
fail.cpp factor.cpp fderivative.cpp function.cpp idx.cpp indexed.cpp inifcns.cpp \
- inifcns_trans.cpp inifcns_gamma.cpp inifcns_nstdsums.cpp \
+ inifcns_trans.cpp inifcns_gamma.cpp inifcns_nstdsums.cpp inifcns_elliptic.cpp integration_kernel.cpp \
integral.cpp lst.cpp matrix.cpp mul.cpp ncmul.cpp normal.cpp numeric.cpp \
operators.cpp power.cpp registrar.cpp relational.cpp remember.cpp \
pseries.cpp print.cpp symbol.cpp symmetry.cpp tensor.cpp \
utils.cpp wildcard.cpp \
remember.h utils.h crc32.h hash_seed.h \
+ utils_multi_iterator.h \
parser/parse_binop_rhs.cpp \
parser/parser.cpp \
parser/parse_context.cpp \
ginacinclude_HEADERS = ginac.h add.h archive.h assertion.h basic.h class_info.h \
clifford.h color.h constant.h container.h ex.h excompiler.h expair.h expairseq.h \
exprseq.h fail.h factor.h fderivative.h flags.h function.h hash_map.h idx.h indexed.h \
- inifcns.h integral.h lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h \
+ inifcns.h integration_kernel.h integral.h lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h \
power.h print.h pseries.h ptr.h registrar.h relational.h structure.h \
symbol.h symmetry.h tensor.h version.h wildcard.h compiler.h \
parser/parser.h \
#include "factor.h"
+#include "integration_kernel.h"
+
#include "excompiler.h"
#ifndef IN_GINAC
return is_the_function<psi1_SERIAL>(x) || is_the_function<psi2_SERIAL>(x);
}
+/** Complete elliptic integral of the first kind. */
+DECLARE_FUNCTION_1P(EllipticK)
+
+/** Complete elliptic integral of the second kind. */
+DECLARE_FUNCTION_1P(EllipticE)
+
+// overloading at work: we cannot use the macros here
+/** Iterated integral. */
+class iterated_integral2_SERIAL { public: static unsigned serial; };
+template<typename T1, typename T2>
+inline function iterated_integral(const T1& kernel_lst, const T2& lambda) {
+ return function(iterated_integral2_SERIAL::serial, ex(kernel_lst), ex(lambda));
+}
+/** Iterated integral with explicit truncation. */
+class iterated_integral3_SERIAL { public: static unsigned serial; };
+template<typename T1, typename T2, typename T3>
+inline function iterated_integral(const T1& kernel_lst, const T2& lambda, const T3& N_trunc) {
+ return function(iterated_integral3_SERIAL::serial, ex(kernel_lst), ex(lambda), ex(N_trunc));
+}
+class iterated_integral_SERIAL;
+template<> inline bool is_the_function<iterated_integral_SERIAL>(const ex& x)
+{
+ return is_the_function<iterated_integral2_SERIAL>(x) || is_the_function<iterated_integral3_SERIAL>(x);
+}
+
+
/** Factorial function. */
DECLARE_FUNCTION_1P(factorial)
--- /dev/null
+/** @file inifcns_elliptic.cpp
+ *
+ * Implementation of some special functions related to elliptic curves
+ *
+ * The functions are:
+ * complete elliptic integral of the first kind EllipticK(k)
+ * complete elliptic integral of the second kind EllipticE(k)
+ * iterated integral iterated_integral(a,y) or iterated_integral(a,y,N_trunc)
+ *
+ * Some remarks:
+ *
+ * - All formulae used can be looked up in the following publication:
+ * [WW] Numerical evaluation of iterated integrals related to elliptic Feynman integrals, M.Walden, S.Weinzierl, arXiv:2010.xxxxx
+ *
+ * - When these routines and methods are used for scientific work that leads to publication in a scientific journal,
+ * please refer to this program as :
+ * M.Walden, S.Weinzierl, "Numerical evaluation of iterated integrals related to elliptic Feynman integrals", arXiv:2010.xxxxx
+ *
+ * - As these routines build on the core part of GiNaC, it is also polite to acknowledge
+ * C. Bauer, A. Frink, R. Kreckel, "Introduction to the GiNaC Framework for Symbolic Computation within the C++ Programming Language",
+ * J. Symbolic Computations 33, 1 (2002), cs.sc/0004015
+ *
+ */
+
+/*
+ * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "inifcns.h"
+
+#include "add.h"
+#include "constant.h"
+#include "lst.h"
+#include "mul.h"
+#include "numeric.h"
+#include "operators.h"
+#include "power.h"
+#include "pseries.h"
+#include "relational.h"
+#include "symbol.h"
+#include "utils.h"
+#include "wildcard.h"
+
+#include "integration_kernel.h"
+#include "utils_multi_iterator.h"
+
+#include <cln/cln.h>
+#include <sstream>
+#include <stdexcept>
+#include <vector>
+#include <cmath>
+
+namespace GiNaC {
+
+
+//////////////////////////////////////////////////////////////////////
+//
+// Complete elliptic integrals
+//
+// helper functions
+//
+//////////////////////////////////////////////////////////////////////
+
+
+// anonymous namespace for helper function
+namespace {
+
+// Computes the arithmetic geometric of two numbers a_0 and b_0
+cln::cl_N arithmetic_geometric_mean(const cln::cl_N & a_0, const cln::cl_N & b_0)
+{
+ cln::cl_N a_old = a_0 * cln::cl_float(1, cln::float_format(Digits));
+ cln::cl_N b_old = b_0 * cln::cl_float(1, cln::float_format(Digits));
+ cln::cl_N a_new;
+ cln::cl_N b_new;
+ cln::cl_N res = a_old;
+ cln::cl_N resbuf;
+ do {
+ resbuf = res;
+
+ a_new = (a_old+b_old)/2;
+ b_new = sqrt(a_old*b_old);
+
+ if ( ( abs(a_new-b_new) > abs(a_new+b_new) )
+ ||
+ ( (abs(a_new-b_new) == abs(a_new+b_new)) && (imagpart(b_new/a_new) <= 0) ) ) {
+ b_new *= -1;
+ }
+
+ res = a_new;
+ a_old = a_new;
+ b_old = b_new;
+
+ } while (res != resbuf);
+ return res;
+}
+
+// Computes
+// a0^2 - sum_{n=0}^infinity 2^{n-1}*c_n^2
+// with
+// c_{n+1} = c_n^2/4/a_{n+1}
+//
+// Needed for the complete elliptic integral of the second kind.
+//
+cln::cl_N agm_helper_second_kind(const cln::cl_N & a_0, const cln::cl_N & b_0, const cln::cl_N & c_0)
+{
+ cln::cl_N a_old = a_0 * cln::cl_float(1, cln::float_format(Digits));
+ cln::cl_N b_old = b_0 * cln::cl_float(1, cln::float_format(Digits));
+ cln::cl_N c_old = c_0 * cln::cl_float(1, cln::float_format(Digits));
+ cln::cl_N a_new;
+ cln::cl_N b_new;
+ cln::cl_N c_new;
+ cln::cl_N res = square(a_old)-square(c_old)/2;
+ cln::cl_N resbuf;
+ cln::cl_N pre = cln::cl_float(1, cln::float_format(Digits));
+ do {
+ resbuf = res;
+
+ a_new = (a_old+b_old)/2;
+ b_new = sqrt(a_old*b_old);
+
+ if ( ( abs(a_new-b_new) > abs(a_new+b_new) )
+ ||
+ ( (abs(a_new-b_new) == abs(a_new+b_new)) && (imagpart(b_new/a_new) <= 0) ) ) {
+ b_new *= -1;
+ }
+
+ c_new = square(c_old)/4/a_new;
+
+ res -= pre*square(c_new);
+
+ a_old = a_new;
+ b_old = b_new;
+ c_old = c_new;
+ pre *= 2;
+
+ } while (res != resbuf);
+ return res;
+}
+
+
+} // end of anonymous namespace
+
+
+//////////////////////////////////////////////////////////////////////
+//
+// Complete elliptic integrals
+//
+// GiNaC function
+//
+//////////////////////////////////////////////////////////////////////
+
+static ex EllipticK_evalf(const ex& k)
+{
+ if ( !k.info(info_flags::numeric) ) {
+ return EllipticK(k).hold();
+ }
+
+ cln::cl_N kbar = sqrt(1-square(ex_to<numeric>(k).to_cl_N()));
+
+ ex result = Pi/2/numeric(arithmetic_geometric_mean(1,kbar));
+
+ return result.evalf();
+}
+
+
+static ex EllipticK_eval(const ex& k)
+{
+ if (k == _ex0) {
+ return Pi/2;
+ }
+
+ if ( k.info(info_flags::numeric) && !k.info(info_flags::crational) ) {
+ return EllipticK(k).evalf();
+ }
+
+ return EllipticK(k).hold();
+}
+
+
+static ex EllipticK_deriv(const ex& k, unsigned deriv_param)
+{
+ return -EllipticK(k)/k + EllipticE(k)/k/(1-k*k);
+}
+
+
+static ex EllipticK_series(const ex& k, const relational& rel, int order, unsigned options)
+{
+ const ex k_pt = k.subs(rel, subs_options::no_pattern);
+
+ if (k_pt == _ex0) {
+ const symbol s;
+ ex ser;
+ // manually construct the primitive expansion
+ for (int i=0; i<(order+1)/2; ++i)
+ {
+ ser += Pi/2 * numeric(cln::square(cln::binomial(2*i,i))) * pow(s/4,2*i);
+ }
+ // substitute the argument's series expansion
+ ser = ser.subs(s==k.series(rel, order), subs_options::no_pattern);
+ // maybe that was terminating, so add a proper order term
+ epvector nseq { expair(Order(_ex1), order) };
+ ser += pseries(rel, std::move(nseq));
+ // reexpanding it will collapse the series again
+ return ser.series(rel, order);
+ }
+
+ if ( (k_pt == _ex1) || (k_pt == _ex_1) ) {
+ throw std::runtime_error("EllipticK_series: don't know how to do the series expansion at this point!");
+ }
+
+ // all other cases
+ throw do_taylor();
+}
+
+static void EllipticK_print_latex(const ex& k, const print_context& c)
+{
+ c.s << "\\mathrm{K}(";
+ k.print(c);
+ c.s << ")";
+}
+
+
+REGISTER_FUNCTION(EllipticK,
+ evalf_func(EllipticK_evalf).
+ eval_func(EllipticK_eval).
+ derivative_func(EllipticK_deriv).
+ series_func(EllipticK_series).
+ print_func<print_latex>(EllipticK_print_latex).
+ do_not_evalf_params());
+
+
+static ex EllipticE_evalf(const ex& k)
+{
+ if ( !k.info(info_flags::numeric) ) {
+ return EllipticE(k).hold();
+ }
+
+ cln::cl_N kbar = sqrt(1-square(ex_to<numeric>(k).to_cl_N()));
+
+ ex result = Pi/2 * numeric( agm_helper_second_kind(1,kbar,ex_to<numeric>(k).to_cl_N()) / arithmetic_geometric_mean(1,kbar) );
+
+ return result.evalf();
+}
+
+
+static ex EllipticE_eval(const ex& k)
+{
+ if (k == _ex0) {
+ return Pi/2;
+ }
+
+ if ( (k == _ex1) || (k == _ex_1) ) {
+ return 1;
+ }
+
+ if ( k.info(info_flags::numeric) && !k.info(info_flags::crational) ) {
+ return EllipticE(k).evalf();
+ }
+
+ return EllipticE(k).hold();
+}
+
+
+static ex EllipticE_deriv(const ex& k, unsigned deriv_param)
+{
+ return -EllipticK(k)/k + EllipticE(k)/k;
+}
+
+
+static ex EllipticE_series(const ex& k, const relational& rel, int order, unsigned options)
+{
+ const ex k_pt = k.subs(rel, subs_options::no_pattern);
+
+ if (k_pt == _ex0) {
+ const symbol s;
+ ex ser;
+ // manually construct the primitive expansion
+ for (int i=0; i<(order+1)/2; ++i)
+ {
+ ser -= Pi/2 * numeric(cln::square(cln::binomial(2*i,i)))/(2*i-1) * pow(s/4,2*i);
+ }
+ // substitute the argument's series expansion
+ ser = ser.subs(s==k.series(rel, order), subs_options::no_pattern);
+ // maybe that was terminating, so add a proper order term
+ epvector nseq { expair(Order(_ex1), order) };
+ ser += pseries(rel, std::move(nseq));
+ // reexpanding it will collapse the series again
+ return ser.series(rel, order);
+ }
+
+ if ( (k_pt == _ex1) || (k_pt == _ex_1) ) {
+ throw std::runtime_error("EllipticE_series: don't know how to do the series expansion at this point!");
+ }
+
+ // all other cases
+ throw do_taylor();
+}
+
+static void EllipticE_print_latex(const ex& k, const print_context& c)
+{
+ c.s << "\\mathrm{K}(";
+ k.print(c);
+ c.s << ")";
+}
+
+
+REGISTER_FUNCTION(EllipticE,
+ evalf_func(EllipticE_evalf).
+ eval_func(EllipticE_eval).
+ derivative_func(EllipticE_deriv).
+ series_func(EllipticE_series).
+ print_func<print_latex>(EllipticE_print_latex).
+ do_not_evalf_params());
+
+
+//////////////////////////////////////////////////////////////////////
+//
+// Iterated integrals
+//
+// helper functions
+//
+//////////////////////////////////////////////////////////////////////
+
+// anonymous namespace for helper function
+namespace {
+
+// performs the actual series summation for an iterated integral
+cln::cl_N iterated_integral_do_sum(const std::vector<int> & m, const std::vector<const integration_kernel *> & kernel, const cln::cl_N & lambda, int N_trunc)
+{
+ if ( cln::zerop(lambda) ) {
+ return 0;
+ }
+
+ cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
+
+ const int depth = m.size();
+
+ cln::cl_N res = 0;
+ cln::cl_N resbuf;
+ cln::cl_N subexpr;
+
+ if ( N_trunc == 0 ) {
+ // sum until precision is reached
+ bool flag_accidental_zero = false;
+
+ int N = 1;
+
+ do {
+ resbuf = res;
+
+ if ( depth > 1 ) {
+ subexpr = 0;
+ multi_iterator_ordered_eq<int> i_multi(1,N+1,depth-1);
+ for( i_multi.init(); !i_multi.overflow(); i_multi++) {
+ cln::cl_N tt = one;
+ for (int j=1; j<depth; j++) {
+ if ( j==1 ) {
+ tt = tt * kernel[0]->series_coeff(N-i_multi[depth-2]) / cln::expt(cln::cl_I(i_multi[depth-2]),m[1]);
+ }
+ else {
+ tt = tt * kernel[j-1]->series_coeff(i_multi[depth-j]-i_multi[depth-j-1]) / cln::expt(cln::cl_I(i_multi[depth-j-1]),m[j]);
+ }
+ }
+ tt = tt * kernel[depth-1]->series_coeff(i_multi[0]);
+ subexpr += tt;
+ }
+ }
+ else {
+ // depth == 1
+ subexpr = kernel[0]->series_coeff(N) * one;
+ }
+ flag_accidental_zero = cln::zerop(subexpr);
+ res += cln::expt(lambda, N) / cln::expt(cln::cl_I(N),m[0]) * subexpr;
+ N++;
+
+ } while ( (res != resbuf) || flag_accidental_zero );
+ }
+ else {
+ // N_trunc > 0, sum up the first N_trunc terms
+ for (int N=1; N<=N_trunc; N++) {
+ if ( depth > 1 ) {
+ subexpr = 0;
+ multi_iterator_ordered_eq<int> i_multi(1,N+1,depth-1);
+ for( i_multi.init(); !i_multi.overflow(); i_multi++) {
+ cln::cl_N tt = one;
+ for (int j=1; j<depth; j++) {
+ if ( j==1 ) {
+ tt = tt * kernel[0]->series_coeff(N-i_multi[depth-2]) / cln::expt(cln::cl_I(i_multi[depth-2]),m[1]);
+ }
+ else {
+ tt = tt * kernel[j-1]->series_coeff(i_multi[depth-j]-i_multi[depth-j-1]) / cln::expt(cln::cl_I(i_multi[depth-j-1]),m[j]);
+ }
+ }
+ tt = tt * kernel[depth-1]->series_coeff(i_multi[0]);
+ subexpr += tt;
+ }
+ }
+ else {
+ // depth == 1
+ subexpr = kernel[0]->series_coeff(N) * one;
+ }
+ res += cln::expt(lambda, N) / cln::expt(cln::cl_I(N),m[0]) * subexpr;
+ }
+ }
+
+ return res;
+}
+
+// figure out the number of basic_log_kernels before a non-basic_log_kernel
+cln::cl_N iterated_integral_prepare_m_lst(const std::vector<const integration_kernel *> & kernel_in, const cln::cl_N & lambda, int N_trunc)
+{
+ size_t depth = kernel_in.size();
+
+ std::vector<int> m;
+ m.reserve(depth);
+
+ std::vector<const integration_kernel *> kernel;
+ kernel.reserve(depth);
+
+ int n = 1;
+
+ for (const auto & it : kernel_in) {
+ if ( is_a<basic_log_kernel>(*it) ) {
+ n++;
+ }
+ else {
+ m.push_back(n);
+ kernel.push_back( &ex_to<integration_kernel>(*it) );
+ n = 1;
+ }
+ }
+
+ cln::cl_N result = iterated_integral_do_sum(m, kernel, lambda, N_trunc);
+
+ return result;
+}
+
+// shuffle to remove trailing zeros,
+// integration kernels, which are not basic_log_kernels, are treated as regularised kernels
+cln::cl_N iterated_integral_shuffle(const std::vector<const integration_kernel *> & kernel, const cln::cl_N & lambda, int N_trunc)
+{
+ cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
+
+ const size_t depth = kernel.size();
+
+ size_t i_trailing = 0;
+ for (size_t i=0; i<depth; i++) {
+ if ( !(is_a<basic_log_kernel>(*(kernel[i]))) ) {
+ i_trailing = i+1;
+ }
+ }
+
+ if ( i_trailing == 0 ) {
+ return cln::expt(cln::log(lambda), depth) / cln::factorial(depth) * one;
+ }
+
+ if ( i_trailing == depth ) {
+ return iterated_integral_prepare_m_lst(kernel, lambda, N_trunc);
+ }
+
+ // shuffle
+ std::vector<const integration_kernel *> a,b;
+ for (size_t i=0; i<i_trailing; i++) {
+ a.push_back(kernel[i]);
+ }
+ for (size_t i=i_trailing; i<depth; i++) {
+ b.push_back(kernel[i]);
+ }
+
+ cln::cl_N result = iterated_integral_prepare_m_lst(a, lambda, N_trunc) * cln::expt(cln::log(lambda), depth-i_trailing) / cln::factorial(depth-i_trailing);
+ multi_iterator_shuffle_prime<const integration_kernel *> i_shuffle(a,b);
+ for( i_shuffle.init(); !i_shuffle.overflow(); i_shuffle++) {
+ std::vector<const integration_kernel *> new_kernel;
+ new_kernel.reserve(depth);
+ for (size_t i=0; i<depth; i++) {
+ new_kernel.push_back(i_shuffle[i]);
+ }
+
+ result -= iterated_integral_shuffle(new_kernel, lambda, N_trunc);
+ }
+
+ return result;
+}
+
+} // end of anonymous namespace
+
+//////////////////////////////////////////////////////////////////////
+//
+// Iterated integrals
+//
+// GiNaC function
+//
+//////////////////////////////////////////////////////////////////////
+
+static ex iterated_integral_evalf_impl(const ex& kernel_lst, const ex& lambda, const ex& N_trunc)
+{
+ // sanity check
+ if ((!kernel_lst.info(info_flags::list)) || (!lambda.evalf().info(info_flags::numeric)) || (!N_trunc.info(info_flags::nonnegint))) {
+ return iterated_integral(kernel_lst,lambda,N_trunc).hold();
+ }
+
+ lst k_lst = ex_to<lst>(kernel_lst);
+
+ bool flag_not_numeric = false;
+ for (const auto & it : k_lst) {
+ if ( !is_a<integration_kernel>(it) ) {
+ flag_not_numeric = true;
+ }
+ }
+ if ( flag_not_numeric) {
+ return iterated_integral(kernel_lst,lambda,N_trunc).hold();
+ }
+
+ for (const auto & it : k_lst) {
+ if ( !(ex_to<integration_kernel>(it).is_numeric()) ) {
+ flag_not_numeric = true;
+ }
+ }
+ if ( flag_not_numeric) {
+ return iterated_integral(kernel_lst,lambda,N_trunc).hold();
+ }
+
+ // now we know that iterated_integral gives a number
+
+ int N_trunc_int = ex_to<numeric>(N_trunc).to_int();
+
+ // check trailing zeros
+ const size_t depth = k_lst.nops();
+
+ std::vector<const integration_kernel *> kernel_vec;
+ kernel_vec.reserve(depth);
+
+ for (const auto & it : k_lst) {
+ kernel_vec.push_back( &ex_to<integration_kernel>(it) );
+ }
+
+ size_t i_trailing = 0;
+ for (size_t i=0; i<depth; i++) {
+ if ( !(kernel_vec[i]->has_trailing_zero()) ) {
+ i_trailing = i+1;
+ }
+ }
+
+ // split integral into regularised integrals and trailing basic_log_kernels
+ // non-basic_log_kernels are treated as regularised kernels in call to iterated_integral_shuffle
+ cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
+ cln::cl_N lambda_cln = ex_to<numeric>(lambda.evalf()).to_cl_N();
+ basic_log_kernel L0 = basic_log_kernel();
+
+ cln::cl_N result;
+ if ( is_a<basic_log_kernel>(*(kernel_vec[depth-1])) ) {
+ result = 0;
+ }
+ else {
+ result = iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
+ }
+
+ cln::cl_N coeff = one;
+ for (size_t i_plus=depth; i_plus>i_trailing; i_plus--) {
+ coeff = coeff * kernel_vec[i_plus-1]->series_coeff(0);
+ kernel_vec[i_plus-1] = &L0;
+ if ( i_plus==i_trailing+1 ) {
+ result += coeff * iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
+ }
+ else {
+ if ( !(is_a<basic_log_kernel>(*(kernel_vec[i_plus-2]))) ) {
+ result += coeff * iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
+ }
+ }
+ }
+
+ return numeric(result);
+}
+
+static ex iterated_integral2_evalf(const ex& kernel_lst, const ex& lambda)
+{
+ return iterated_integral_evalf_impl(kernel_lst,lambda,0);
+}
+
+static ex iterated_integral3_evalf(const ex& kernel_lst, const ex& lambda, const ex& N_trunc)
+{
+ return iterated_integral_evalf_impl(kernel_lst,lambda,N_trunc);
+}
+
+static ex iterated_integral2_eval(const ex& kernel_lst, const ex& lambda)
+{
+ if ( lambda.info(info_flags::numeric) && !lambda.info(info_flags::crational) ) {
+ return iterated_integral(kernel_lst,lambda).evalf();
+ }
+
+ return iterated_integral(kernel_lst,lambda).hold();
+}
+
+static ex iterated_integral3_eval(const ex& kernel_lst, const ex& lambda, const ex& N_trunc)
+{
+ if ( lambda.info(info_flags::numeric) && !lambda.info(info_flags::crational) ) {
+ return iterated_integral(kernel_lst,lambda,N_trunc).evalf();
+ }
+
+ return iterated_integral(kernel_lst,lambda,N_trunc).hold();
+}
+
+unsigned iterated_integral2_SERIAL::serial =
+ function::register_new(function_options("iterated_integral", 2).
+ eval_func(iterated_integral2_eval).
+ evalf_func(iterated_integral2_evalf).
+ do_not_evalf_params().
+ overloaded(2));
+
+unsigned iterated_integral3_SERIAL::serial =
+ function::register_new(function_options("iterated_integral", 3).
+ eval_func(iterated_integral3_eval).
+ evalf_func(iterated_integral3_evalf).
+ do_not_evalf_params().
+ overloaded(2));
+
+} // namespace GiNaC
+
--- /dev/null
+/** @file integration_kernel.cpp
+ *
+ * Implementation of GiNaC's integration kernels for iterated integrals. */
+
+/*
+ * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#include "integration_kernel.h"
+#include "add.h"
+#include "mul.h"
+#include "operators.h"
+#include "power.h"
+#include "relational.h"
+#include "symbol.h"
+#include "constant.h"
+#include "numeric.h"
+#include "function.h"
+#include "pseries.h"
+#include "utils.h"
+#include "inifcns.h"
+
+#include <iostream>
+#include <stdexcept>
+#include <cln/cln.h>
+
+
+namespace GiNaC {
+
+// anonymous namespace for helper function
+namespace {
+
+/**
+ *
+ * Returns the Kronecker symbol in the case where
+ * a: integer
+ * n: unit or prime number
+ *
+ * If n is an odd prime number, the routine returns the Legendre symbol.
+ *
+ * If n is a unit (e.g n equals 1 or -1) or if n is an even prime number (the only case is n=2)
+ * the routine returns the special values defined below.
+ *
+ * Ref.: Toshitsune Miyake, Modular Forms, Chapter 3.1
+ *
+ */
+numeric kronecker_symbol_prime(const numeric & a, const numeric & n)
+{
+ if ( n == 1 ) {
+ return 1;
+ }
+ else if ( n == -1 ) {
+ if ( a >= 0 ) {
+ return 1;
+ }
+ else {
+ return -1;
+ }
+ }
+ else if ( n == 2 ) {
+ if ( GiNaC::smod(a,8) == 1 ) {
+ return 1;
+ }
+ else if ( GiNaC::smod(a,8) == -1 ) {
+ return 1;
+ }
+ else if ( GiNaC::smod(a,8) == 3 ) {
+ return -1;
+ }
+ else if ( GiNaC::smod(a,8) == -3 ) {
+ return -1;
+ }
+ else {
+ return 0;
+ }
+ }
+
+ // n is an odd prime number
+ return GiNaC::smod( pow(a,(n-1)/2), n);
+}
+
+/**
+ *
+ * n: positive integer
+ *
+ * a: discriminant of a quadratic field, defines primitive character phi
+ * b: discriminant of a quadratic field, defines primitive character psi
+ * L=|a|: conductor of primitive character phi
+ * M=|b|: conductor of primitive character psi
+ *
+ * k: modular weight
+ *
+ * This function computes
+ * \f[
+ * \sum\limits_{d | n} \psi(d) \phi(n/d) d^{k-1}
+ * \f]
+ *
+ * Ref.: Eq.(5.3.1), William A. Stein, Modular Forms, A computational Approach;
+ *
+ * Eq.(32), arXiv:1704.08895
+ *
+ */
+numeric divisor_function(const numeric & n, const numeric & a, const numeric & b, const numeric & k)
+{
+ ex res = 0;
+
+ for (numeric i1=1; i1<=n; i1++) {
+ if ( irem(n,i1) == 0 ) {
+ numeric ratio = n/i1;
+ res += primitive_dirichlet_character(ratio,a) * primitive_dirichlet_character(i1,b) * pow(i1,k-1);
+ }
+ }
+
+ return ex_to<numeric>(res);
+}
+
+/**
+ *
+ * k: modular weight
+ *
+ * a: discriminant of a quadratic field, defines primitive character phi
+ * b: discriminant of a quadratic field, defines primitive character psi
+ * L=|a|: conductor of primitive character phi
+ * M=|b|: conductor of primitive character psi
+ *
+ * This function computes the constant term of the q-expansion of an Eisenstein series.
+ *
+ * The coefficient is given by the equation below eq.(5.3.1) in William A. Stein, Modular Forms, A computational Approach;
+ *
+ * or by eq.(33), arXiv:1704.08895
+ *
+ */
+numeric coefficient_a0(const numeric & k, const numeric & a, const numeric & b)
+{
+ ex conductor = abs(a);
+
+ numeric a0;
+ if ( conductor == 1 ) {
+ a0 = -numeric(1,2)/k*generalised_Bernoulli_number(k,b);
+ }
+ else {
+ a0 = 0;
+ }
+
+ return a0;
+}
+
+/**
+ *
+ * k: modular weight
+ * q: exp(2 Pi I tau/M)
+ *
+ * a: discriminant of a quadratic field, defines primitive character phi
+ * b: discriminant of a quadratic field, defines primitive character psi
+ * L=|a|: conductor of primitive character phi
+ * M=|b|: conductor of primitive character psi
+ *
+ * N: truncation order
+ *
+ * Returns the q-expansion of an Eisenstein to order N (the q^(N-1)-term is included, q^N is neglected).
+ *
+ * Ref.: Eq.(5.3.1), William A. Stein, Modular Forms, A computational Approach;
+ *
+ * Eq.(32), arXiv:1704.08895
+ *
+ */
+ex eisenstein_series(const numeric & k, const ex & q, const numeric & a, const numeric & b, const numeric & N)
+{
+ ex res = coefficient_a0(k,a,b);
+
+ for (numeric i1=1; i1<N; i1++) {
+ res += divisor_function(i1,a,b,k) * pow(q,i1);
+ }
+
+ return res;
+}
+
+/**
+ *
+ * Returns the q_N-expansion of the Eisenstein series E_{k,a,b}(K tau_N)
+ *
+ */
+ex E_eisenstein_series(const ex & q, const numeric & k, const numeric & N_level, const numeric & a, const numeric & b, const numeric & K, const numeric & N_order)
+{
+ int N_order_int = N_order.to_int();
+
+ ex res = eisenstein_series(k,pow(q,K),a,b,iquo(N_order,K));
+
+ res += Order(pow(q,N_order_int));
+ res = res.series(q,N_order_int);
+
+ return res;
+}
+
+/**
+ *
+ * In weight 2 we have a special case for trivial Dirichlet characters:
+ * Returns the q_N-expansion of the Eisenstein series E_{2,1,1}(tau_N) - K E_{2,1,1}(K tau_N).
+ *
+ */
+ex B_eisenstein_series(const ex & q, const numeric & N_level, const numeric & K, const numeric & N_order)
+{
+ int N_order_int = N_order.to_int();
+
+ ex res = eisenstein_series(2,q,1,1,N_order) - K*eisenstein_series(2,pow(q,K),1,1,iquo(N_order,K));
+
+ res += Order(pow(q,N_order_int));
+ res = res.series(q,N_order_int);
+
+ return res;
+}
+
+/**
+ *
+ * A helper function to expand polynomials in Eisenstein series.
+ *
+ */
+struct subs_q_expansion : public map_function
+{
+ subs_q_expansion(const ex & arg_qbar, int arg_order) : qbar(arg_qbar), order(arg_order)
+ {}
+
+ ex operator()(const ex & e)
+ {
+ if ( is_a<Eisenstein_kernel>(e) || is_a<Eisenstein_h_kernel>(e) ) {
+ return series_to_poly(e.series(qbar,order));
+ }
+ else {
+ return e.map(*this);
+ }
+ }
+
+ ex qbar;
+ int order;
+};
+
+/**
+ *
+ * \f[
+ * Li_{-n}(x), n>=0
+ * \f]
+ *
+ * This is a rational function in x.
+ *
+ * To speed things up, we cache it.
+ *
+ */
+class Li_negative
+{
+ // ctors
+public:
+ Li_negative();
+
+ // non-virtual functions
+public:
+ ex get_symbolic_value(int n, const ex & x_val);
+ ex get_numerical_value(int n, const ex & x_val);
+
+ // member variables :
+protected:
+ static std::vector<ex> cache_vec;
+ static symbol x;
+};
+
+
+Li_negative::Li_negative() {}
+
+ex Li_negative::get_symbolic_value(int n, const ex & x_val)
+{
+ int n_cache = cache_vec.size();
+
+ if ( n >= n_cache ) {
+ for (int j=n_cache; j<=n; j++) {
+ ex f;
+ if ( j == 0 ) {
+ f = x/(1-x);
+ }
+ else {
+ f = normal( x*diff(cache_vec[j-1],x));
+ }
+ cache_vec.push_back( f );
+ }
+ }
+
+ return cache_vec[n].subs(x==x_val);
+}
+
+ex Li_negative::get_numerical_value(int n, const ex & x_val)
+{
+ symbol x_symb("x_symb");
+
+ ex f = this->get_symbolic_value(n,x_symb);
+
+ ex res = f.subs(x_symb==x_val).evalf();
+
+ return res;
+}
+
+// initialise static data members
+std::vector<ex> Li_negative::cache_vec;
+symbol Li_negative::x = symbol("x");
+
+
+} // end of anonymous namespace
+
+/**
+ *
+ * Returns the decomposition of the positive integer n into prime numbers
+ * in the form
+ * lst( lst(p1,...,pr), lst(a1,...,ar) )
+ * such that
+ * n = p1^a1 ... pr^ar.
+ *
+ */
+ex ifactor(const numeric & n)
+{
+ if ( !n.is_pos_integer() ) throw (std::runtime_error("ifactor(): argument not a positive integer"));
+
+ lst p_lst, exp_lst;
+
+ // implementation for small integers
+ numeric n_temp = n;
+ for (numeric p=2; p<=n; p++) {
+ if ( p.info(info_flags::prime) ) {
+ numeric exp_temp = 0;
+ while ( irem(n_temp, p) == 0 ) {
+ n_temp = n_temp/p;
+ exp_temp++;
+ }
+ if ( exp_temp>0 ) {
+ p_lst.append(p);
+ exp_lst.append(exp_temp);
+ }
+ }
+ if ( n_temp == 1 ) break;
+ }
+
+ if ( n_temp != 1 ) throw (std::runtime_error("ifactor(): probabilistic primality test failed"));
+
+ lst res = {p_lst,exp_lst};
+
+ return res;
+}
+
+/**
+ *
+ * Returns true if the integer n is either one or the discriminant of a quadratic number field.
+ *
+ * Returns false otherwise.
+ *
+ * Ref.: Toshitsune Miyake, Modular Forms, Chapter 3.1
+ *
+ */
+bool is_discriminant_of_quadratic_number_field(const numeric & n)
+{
+ if ( n == 0 ) {
+ return false;
+ }
+
+ if ( n == 1 ) {
+ return true;
+ }
+
+ lst prime_factorisation = ex_to<lst>(ifactor(abs(n)));
+ lst p_lst = ex_to<lst>(prime_factorisation.op(0));
+ lst e_lst = ex_to<lst>(prime_factorisation.op(1));
+
+ size_t n_primes = p_lst.nops();
+
+ if ( n_primes > 0 ) {
+ // take the last prime
+ numeric p = ex_to<numeric>(p_lst.op(n_primes-1));
+
+ if ( p.is_odd() ) {
+ if ( e_lst.op(n_primes-1) != 1 ) {
+ return false;
+ }
+
+ numeric pstar = p;
+ if ( mod(p,4) == 3 ) {
+ pstar = -p;
+ }
+ return is_discriminant_of_quadratic_number_field(n/pstar);
+ }
+ }
+ // power of two now
+ if ( (n==-4) || (n==-8) || (n==8) || (n==-32) || (n==32) || (n==-64) || (n==128) ) {
+ return true;
+ }
+
+ return false;
+}
+
+/**
+ *
+ * Returns the Kronecker symbol
+ * a: integer
+ * n: integer
+ *
+ * This routine defines
+ * kronecker_symbol(1,0) = 1
+ * kronecker_symbol(-1,0) = 1
+ * kronecker_symbol(a,0) = 0, a != 1,-1
+ *
+ * In particular
+ * kronecker_symbol(-1,0) = 1 (in agreement with Sage)
+ *
+ * Ref.: Toshitsune Miyake, Modular Forms, Chapter 3.1
+ *
+ */
+numeric kronecker_symbol(const numeric & a, const numeric & n)
+{
+ // case n=0 first, include kronecker_symbol(0,0)=0
+ if ( n == 0 ) {
+ if ( (a == 1) || (a == -1) ) {
+ return 1;
+ }
+ else {
+ return 0;
+ }
+ }
+
+ numeric unit = 1;
+ numeric n_pos = n;
+ if ( n_pos<0 ) {
+ unit = -1;
+ n_pos = -n;
+ }
+
+ ex res = kronecker_symbol_prime(a,unit);
+
+ numeric n_odd = n_pos;
+ numeric alpha = 0;
+ while ( n_odd.is_even() ) {
+ alpha++;
+ n_odd = n_odd/2;
+ }
+ if ( alpha>0 ) {
+ res *= pow(kronecker_symbol_prime(a,2),alpha);
+ }
+
+ lst temp_lst = ex_to<lst>(ifactor(n_odd));
+ lst prime_lst = ex_to<lst>(temp_lst.op(0));
+ lst expo_lst = ex_to<lst>(temp_lst.op(1));
+
+ for (auto it_p = prime_lst.begin(), it_e = expo_lst.begin(); it_p != prime_lst.end(); it_p++, it_e++) {
+ res *= pow(kronecker_symbol_prime(a,ex_to<numeric>(*it_p)),ex_to<numeric>(*it_e));
+ }
+
+ return ex_to<numeric>(res);
+}
+
+/**
+ *
+ * Defines a primitive Dirichlet character through the Kronecker symbol.
+ *
+ * n: integer
+ * a: discriminant of a quadratic field
+ * |a|: conductor
+ *
+ * The character takes the values -1,0,1.
+ *
+ */
+numeric primitive_dirichlet_character(const numeric & n, const numeric & a)
+{
+ return kronecker_symbol(a,n);
+}
+
+/**
+ *
+ * Defines a Dirichlet character through the Kronecker symbol.
+ *
+ * n: integer
+ * a: discriminant of a quadratic field
+ * |a|: conductor
+ * N: modulus, needs to be multiple of |a|
+ *
+ * The character takes the values -1,0,1.
+ *
+ */
+numeric dirichlet_character(const numeric & n, const numeric & a, const numeric & N)
+{
+ if ( gcd(n,N) == 1 ) {
+ return primitive_dirichlet_character(n,a);
+ }
+
+ return 0;
+}
+
+/**
+ *
+ * The generalised Bernoulli number.
+ *
+ * k: index / modular weight
+ *
+ * b: discriminant of a quadratic field, defines primitive character psi
+ * M=|b|: conductor of primitive character psi
+ *
+ * The generalised Bernoulli number is computed from the series expansion of the generating function.
+ * The generating function is given in eq.(34), arXiv:1704.08895
+ *
+ */
+numeric generalised_Bernoulli_number(const numeric & k, const numeric & b)
+{
+ int k_int = k.to_int();
+
+ symbol x("x");
+
+ numeric conductor = abs(b);
+
+ ex gen_fct = 0;
+ for (numeric i1=1; i1<=conductor; i1++) {
+ gen_fct += primitive_dirichlet_character(i1,b) * x*exp(i1*x)/(exp(conductor*x)-1);
+ }
+
+ gen_fct = series_to_poly(gen_fct.series(x,k_int+1));
+
+ ex B = factorial(k) * gen_fct.coeff(x,k_int);
+
+ return ex_to<numeric>(B);
+}
+
+/**
+ *
+ * The Bernoulli polynomials
+ *
+ */
+ex Bernoulli_polynomial(const numeric & k, const ex & x)
+{
+ int k_int = k.to_int();
+
+ symbol t("t");
+
+ ex gen_fct = t*exp(x*t)/(exp(t)-1);
+
+ gen_fct = series_to_poly(gen_fct.series(t,k_int+1));
+
+ ex B = factorial(k) * gen_fct.coeff(t,k_int);
+
+ return B;
+}
+
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integration_kernel, basic,
+ print_func<print_context>(&integration_kernel::do_print))
+
+integration_kernel::integration_kernel() : inherited(), cache_step_size(100), series_vec()
+{
+}
+
+int integration_kernel::compare_same_type(const basic &other) const
+{
+ return 0;
+}
+
+ex integration_kernel::series(const relational & r, int order, unsigned options) const
+{
+ if ( r.rhs() != 0 ) {
+ throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
+ }
+
+ return Laurent_series(r.lhs(),order);
+}
+
+/**
+ *
+ * This routine returns true, if the integration kernel has a trailing zero.
+ *
+ */
+bool integration_kernel::has_trailing_zero(void) const
+{
+ if ( cln::zerop( series_coeff(0) ) ) {
+ return false;
+ }
+
+ return true;
+}
+
+/**
+ *
+ * This routine returns true, if the integration kernel can be evaluated numerically.
+ *
+ */
+bool integration_kernel::is_numeric(void) const
+{
+ return true;
+}
+
+/**
+ *
+ * Subclasses have either to implement series_coeff_impl
+ * or the two methods Laurent_series and uses_Laurent_series.
+ *
+ * The method series_coeff_impl can be used, if a single coefficient can be computed
+ * independently of the others.
+ *
+ * The method Laurent_series can be used, if it is more efficient to compute a Laurent series
+ * in one shot and to determine a range of coefficients from this Laurent series.
+ *
+ */
+cln::cl_N integration_kernel::series_coeff(int i) const
+{
+ int n_vec = series_vec.size();
+
+ if ( i >= n_vec ) {
+ int N = cache_step_size*(i/cache_step_size+1);
+
+ if ( uses_Laurent_series() ) {
+ symbol x("x");
+ // series_vec[0] gives coefficient of 1/z, series_vec[N-1] coefficient of z^(N-2),
+ // thus expansion up to order (N-1) is required
+ ex temp = Laurent_series(x, N-1);
+ for (int j=n_vec; j<N; j++) {
+ series_vec.push_back( ex_to<numeric>(temp.coeff(x,j-1).evalf()).to_cl_N() );
+ }
+ }
+ else {
+ for (int j=n_vec; j<N; j++) {
+ series_vec.push_back( series_coeff_impl(j) );
+ }
+ }
+ }
+
+ return series_vec[i];
+}
+
+/**
+ *
+ * For \f$ \omega = d\lambda \f$ only the coefficient of \f$ \lambda^0 \f$ is non-zero.
+ *
+ * The i-th coefficient corresponds to the power \f$ \lambda^{i-1} \f$.
+ *
+ */
+cln::cl_N integration_kernel::series_coeff_impl(int i) const
+{
+ if ( i == 1 ) {
+ return 1;
+ }
+
+ return 0;
+}
+
+/**
+ *
+ * Returns the Laurent series, starting possibly with the pole term.
+ * Neglected terms are of order \f$ O(x^order) \f$.
+ *
+ */
+ex integration_kernel::Laurent_series(const ex & x, int order) const
+{
+ ex res = 0;
+ for (int n=-1; n<order; n++) {
+ res += numeric(series_coeff(n+1)) * pow(x,n);
+ }
+ res += Order(pow(x,order));
+ res = res.series(x,order);
+
+ return res;
+}
+
+/**
+ *
+ * Evaluates the integrand at lambda.
+ *
+ */
+ex integration_kernel::get_numerical_value(const ex & lambda, int N_trunc) const
+{
+ return get_numerical_value_impl(lambda, 1, 0, N_trunc);
+}
+
+/**
+ *
+ * Returns true, if the coefficients are computed from the Laurent series
+ * (in which case the method Laurent_series needs to be implemented).
+ *
+ * Returns false if this is not the case
+ * (and the class has an implementation of series_coeff_impl).
+ *
+ */
+bool integration_kernel::uses_Laurent_series() const
+{
+ return false;
+}
+
+/**
+ *
+ * Returns the current size of the cache.
+ *
+ */
+size_t integration_kernel::get_cache_size(void) const
+{
+ return series_vec.size();
+}
+
+/**
+ *
+ * Sets the step size by which the cache is increased.
+ *
+ */
+void integration_kernel::set_cache_step(int cache_steps) const
+{
+ cache_step_size = cache_steps;
+}
+
+/**
+ *
+ * Wrapper around series_coeff(i), converts cl_N to numeric.
+ *
+ */
+ex integration_kernel::get_series_coeff(int i) const
+{
+ return numeric(series_coeff(i));
+}
+
+/**
+ *
+ * The actual implementation for computing a numerical value for the integrand.
+ *
+ */
+ex integration_kernel::get_numerical_value_impl(const ex & lambda, const ex & pre, int shift, int N_trunc) const
+{
+ cln::cl_N lambda_cln = ex_to<numeric>(lambda.evalf()).to_cl_N();
+ cln::cl_N pre_cln = ex_to<numeric>(pre.evalf()).to_cl_N();
+
+ cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
+
+ cln::cl_N res = 0;
+ cln::cl_N resbuf;
+ cln::cl_N subexpr;
+
+ if ( N_trunc == 0 ) {
+ // sum until precision is reached
+ bool flag_accidental_zero = false;
+
+ int N = 0;
+
+ do {
+ resbuf = res;
+
+ subexpr = series_coeff(N);
+
+ res += pre_cln * subexpr * cln::expt(lambda_cln,N-1+shift);
+
+ flag_accidental_zero = cln::zerop(subexpr);
+
+ N++;
+ } while ( (res != resbuf) || flag_accidental_zero );
+ }
+ else {
+ // N_trunc > 0, sum up the first N_trunc terms
+ for (int N=0; N<N_trunc; N++) {
+ subexpr = series_coeff(N);
+
+ res += pre_cln * subexpr * cln::expt(lambda_cln,N-1+shift);
+ }
+ }
+
+ return numeric(res);
+}
+
+void integration_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "integration_kernel()";
+}
+
+GINAC_BIND_UNARCHIVER(integration_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic_log_kernel, integration_kernel,
+ print_func<print_context>(&basic_log_kernel::do_print))
+
+basic_log_kernel::basic_log_kernel() : inherited()
+{
+}
+
+int basic_log_kernel::compare_same_type(const basic &other) const
+{
+ return 0;
+}
+
+cln::cl_N basic_log_kernel::series_coeff_impl(int i) const
+{
+ if ( i == 0 ) {
+ return 1;
+ }
+
+ return 0;
+}
+
+void basic_log_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "basic_log_kernel()";
+}
+
+GINAC_BIND_UNARCHIVER(basic_log_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(multiple_polylog_kernel, integration_kernel,
+ print_func<print_context>(&multiple_polylog_kernel::do_print))
+
+multiple_polylog_kernel::multiple_polylog_kernel() : inherited(), z(_ex1)
+{
+}
+
+multiple_polylog_kernel::multiple_polylog_kernel(const ex & arg_z) : inherited(), z(arg_z)
+{
+}
+
+int multiple_polylog_kernel::compare_same_type(const basic &other) const
+{
+ const multiple_polylog_kernel &o = static_cast<const multiple_polylog_kernel &>(other);
+
+ return z.compare(o.z);
+}
+
+size_t multiple_polylog_kernel::nops() const
+{
+ return 1;
+}
+
+ex multiple_polylog_kernel::op(size_t i) const
+{
+ if ( i != 0 ) {
+ throw(std::range_error("multiple_polylog_kernel::op(): out of range"));
+ }
+
+ return z;
+}
+
+ex & multiple_polylog_kernel::let_op(size_t i)
+{
+ ensure_if_modifiable();
+
+ if ( i != 0 ) {
+ throw(std::range_error("multiple_polylog_kernel::let_op(): out of range"));
+ }
+
+ return z;
+}
+
+bool multiple_polylog_kernel::is_numeric(void) const
+{
+ return z.evalf().info(info_flags::numeric);
+}
+
+cln::cl_N multiple_polylog_kernel::series_coeff_impl(int i) const
+{
+ if ( i == 0 ) {
+ return 0;
+ }
+
+ return -cln::expt(ex_to<numeric>(z.evalf()).to_cl_N(),-i);
+}
+
+void multiple_polylog_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "multiple_polylog_kernel(";
+ z.print(c);
+ c.s << ")";
+}
+
+GINAC_BIND_UNARCHIVER(multiple_polylog_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(ELi_kernel, integration_kernel,
+ print_func<print_context>(&ELi_kernel::do_print))
+
+ELi_kernel::ELi_kernel() : inherited(), n(_ex0), m(_ex0), x(_ex0), y(_ex0)
+{
+}
+
+ELi_kernel::ELi_kernel(const ex & arg_n, const ex & arg_m, const ex & arg_x, const ex & arg_y) : inherited(), n(arg_n), m(arg_m), x(arg_x), y(arg_y)
+{
+}
+
+int ELi_kernel::compare_same_type(const basic &other) const
+{
+ const ELi_kernel &o = static_cast<const ELi_kernel &>(other);
+ int cmpval;
+
+ cmpval = n.compare(o.n);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = m.compare(o.m);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = x.compare(o.x);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ return y.compare(o.y);
+}
+
+size_t ELi_kernel::nops() const
+{
+ return 4;
+}
+
+ex ELi_kernel::op(size_t i) const
+{
+ switch (i) {
+ case 0:
+ return n;
+ case 1:
+ return m;
+ case 2:
+ return x;
+ case 3:
+ return y;
+ default:
+ throw (std::out_of_range("ELi_kernel::op() out of range"));
+ }
+}
+
+ex & ELi_kernel::let_op(size_t i)
+{
+ ensure_if_modifiable();
+
+ switch (i) {
+ case 0:
+ return n;
+ case 1:
+ return m;
+ case 2:
+ return x;
+ case 3:
+ return y;
+ default:
+ throw (std::out_of_range("ELi_kernel::let_op() out of range"));
+ }
+}
+
+bool ELi_kernel::is_numeric(void) const
+{
+ return (n.info(info_flags::nonnegint) && m.info(info_flags::numeric) && x.evalf().info(info_flags::numeric) && y.evalf().info(info_flags::numeric));
+}
+
+cln::cl_N ELi_kernel::series_coeff_impl(int i) const
+{
+ if ( i == 0 ) {
+ return 0;
+ }
+
+ int n_int = ex_to<numeric>(n).to_int();
+ int m_int = ex_to<numeric>(m).to_int();
+
+ cln::cl_N x_cln = ex_to<numeric>(x.evalf()).to_cl_N();
+ cln::cl_N y_cln = ex_to<numeric>(y.evalf()).to_cl_N();
+
+ cln::cl_N res_cln = 0;
+
+ for (int j=1; j<=i; j++) {
+ if ( (i % j) == 0 ) {
+ int k = i/j;
+
+ res_cln += cln::expt(x_cln,j)/cln::expt(cln::cl_I(j),n_int) * cln::expt(y_cln,k)/cln::expt(cln::cl_I(k),m_int);
+ }
+ }
+
+ return res_cln;
+}
+
+/**
+ *
+ * Returns the value of ELi_{n,m}(x,y,qbar)
+ *
+ */
+ex ELi_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
+{
+ return get_numerical_value_impl(qbar, 1, 1, N_trunc);
+}
+
+void ELi_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "ELi_kernel(";
+ n.print(c);
+ c.s << ",";
+ m.print(c);
+ c.s << ",";
+ x.print(c);
+ c.s << ",";
+ y.print(c);
+ c.s << ")";
+}
+
+GINAC_BIND_UNARCHIVER(ELi_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Ebar_kernel, integration_kernel,
+ print_func<print_context>(&Ebar_kernel::do_print))
+
+Ebar_kernel::Ebar_kernel() : inherited(), n(_ex0), m(_ex0), x(_ex0), y(_ex0)
+{
+}
+
+Ebar_kernel::Ebar_kernel(const ex & arg_n, const ex & arg_m, const ex & arg_x, const ex & arg_y) : inherited(), n(arg_n), m(arg_m), x(arg_x), y(arg_y)
+{
+}
+
+int Ebar_kernel::compare_same_type(const basic &other) const
+{
+ const Ebar_kernel &o = static_cast<const Ebar_kernel &>(other);
+ int cmpval;
+
+ cmpval = n.compare(o.n);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = m.compare(o.m);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = x.compare(o.x);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ return y.compare(o.y);
+}
+
+size_t Ebar_kernel::nops() const
+{
+ return 4;
+}
+
+ex Ebar_kernel::op(size_t i) const
+{
+ switch (i) {
+ case 0:
+ return n;
+ case 1:
+ return m;
+ case 2:
+ return x;
+ case 3:
+ return y;
+ default:
+ throw (std::out_of_range("Ebar_kernel::op() out of range"));
+ }
+}
+
+ex & Ebar_kernel::let_op(size_t i)
+{
+ ensure_if_modifiable();
+
+ switch (i) {
+ case 0:
+ return n;
+ case 1:
+ return m;
+ case 2:
+ return x;
+ case 3:
+ return y;
+ default:
+ throw (std::out_of_range("Ebar_kernel::let_op() out of range"));
+ }
+}
+
+bool Ebar_kernel::is_numeric(void) const
+{
+ return (n.info(info_flags::nonnegint) && m.info(info_flags::numeric) && x.evalf().info(info_flags::numeric) && y.evalf().info(info_flags::numeric));
+}
+
+cln::cl_N Ebar_kernel::series_coeff_impl(int i) const
+{
+ if ( i == 0 ) {
+ return 0;
+ }
+
+ int n_int = ex_to<numeric>(n).to_int();
+ int m_int = ex_to<numeric>(m).to_int();
+
+ cln::cl_N x_cln = ex_to<numeric>(x.evalf()).to_cl_N();
+ cln::cl_N y_cln = ex_to<numeric>(y.evalf()).to_cl_N();
+
+ cln::cl_N res_cln = 0;
+
+ for (int j=1; j<=i; j++) {
+ if ( (i % j) == 0 ) {
+ int k = i/j;
+
+ res_cln += (cln::expt(x_cln,j)*cln::expt(y_cln,k)-cln::expt(cln::cl_I(-1),n_int+m_int)*cln::expt(x_cln,-j)*cln::expt(y_cln,-k))/cln::expt(cln::cl_I(j),n_int)/cln::expt(cln::cl_I(k),m_int);
+ }
+ }
+
+ return res_cln;
+}
+
+/**
+ *
+ * Returns the value of Ebar_{n,m}(x,y,qbar)
+ *
+ */
+ex Ebar_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
+{
+ return get_numerical_value_impl(qbar, 1, 1, N_trunc);
+}
+
+void Ebar_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "Ebar_kernel(";
+ n.print(c);
+ c.s << ",";
+ m.print(c);
+ c.s << ",";
+ x.print(c);
+ c.s << ",";
+ y.print(c);
+ c.s << ")";
+}
+
+GINAC_BIND_UNARCHIVER(Ebar_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Kronecker_dtau_kernel, integration_kernel,
+ print_func<print_context>(&Kronecker_dtau_kernel::do_print))
+
+Kronecker_dtau_kernel::Kronecker_dtau_kernel() : inherited(), n(_ex0), z(_ex0), K(_ex1), C_norm(_ex1)
+{
+}
+
+Kronecker_dtau_kernel::Kronecker_dtau_kernel(const ex & arg_n, const ex & arg_z, const ex & arg_K, const ex & arg_C_norm) : inherited(), n(arg_n), z(arg_z), K(arg_K), C_norm(arg_C_norm)
+{
+}
+
+int Kronecker_dtau_kernel::compare_same_type(const basic &other) const
+{
+ const Kronecker_dtau_kernel &o = static_cast<const Kronecker_dtau_kernel &>(other);
+ int cmpval;
+
+ cmpval = n.compare(o.n);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = z.compare(o.z);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = K.compare(o.K);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ return C_norm.compare(o.C_norm);
+}
+
+size_t Kronecker_dtau_kernel::nops() const
+{
+ return 4;
+}
+
+ex Kronecker_dtau_kernel::op(size_t i) const
+{
+ switch (i) {
+ case 0:
+ return n;
+ case 1:
+ return z;
+ case 2:
+ return K;
+ case 3:
+ return C_norm;
+ default:
+ throw (std::out_of_range("Kronecker_dtau_kernel::op() out of range"));
+ }
+}
+
+ex & Kronecker_dtau_kernel::let_op(size_t i)
+{
+ ensure_if_modifiable();
+
+ switch (i) {
+ case 0:
+ return n;
+ case 1:
+ return z;
+ case 2:
+ return K;
+ case 3:
+ return C_norm;
+ default:
+ throw (std::out_of_range("Kronecker_dtau_kernel::let_op() out of range"));
+ }
+}
+
+bool Kronecker_dtau_kernel::is_numeric(void) const
+{
+ return (n.info(info_flags::nonnegint) && z.evalf().info(info_flags::numeric) && K.info(info_flags::posint) && C_norm.evalf().info(info_flags::numeric));
+}
+
+cln::cl_N Kronecker_dtau_kernel::series_coeff_impl(int i) const
+{
+ numeric n_num = ex_to<numeric>(n);
+ int n_int = n_num.to_int();
+
+ // case n=0
+ if ( n_num == 0 ) {
+ if ( i == 0 ) {
+ ex res = -C_norm*K;
+
+ return ex_to<numeric>(res.evalf()).to_cl_N();
+ }
+
+ return 0;
+ }
+
+ // case n=1
+ if ( n_num == 1 ) {
+ return 0;
+ }
+
+ // case n>1
+ if ( i == 0 ) {
+ ex res = C_norm*K / factorial(n_num-2) * bernoulli(n_num)/n_num;
+
+ return ex_to<numeric>(res.evalf()).to_cl_N();
+ }
+
+ // n>1, i>0
+
+ // if K>1 the variable i needs to be a multiple of K
+ int K_int = ex_to<numeric>(K).to_int();
+
+ if ( (i % K_int) != 0 ) {
+ return 0;
+ }
+ int i_local = i/K_int;
+
+ ex w = exp(ex_to<numeric>((2*Pi*I*z).evalf()));
+ cln::cl_N w_cln = ex_to<numeric>(w).to_cl_N();
+ cln::cl_N res_cln = 0;
+ for (int j=1; j<=i_local; j++) {
+ if ( (i_local % j) == 0 ) {
+ res_cln += (cln::expt(w_cln,j)+cln::expt(cln::cl_I(-1),n_int)*cln::expt(w_cln,-j)) * cln::expt(cln::cl_I(i_local/j),n_int-1);
+ }
+ }
+ ex pre = -C_norm*K/factorial(n_num-2);
+
+ return ex_to<numeric>(pre.evalf()).to_cl_N() * res_cln;
+}
+
+/**
+ *
+ * Returns the value of the g^(n)(z,K*tau), where tau is given by qbar.
+ *
+ */
+ex Kronecker_dtau_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
+{
+ numeric n_num = ex_to<numeric>(n);
+
+ if ( n_num == 0 ) {
+ return 1;
+ }
+
+ // use the direct formula here
+ if ( n_num == 1 ) {
+ ex wbar = exp(ex_to<numeric>((2*Pi*I*z).evalf()));
+ ex res = -2*Pi*I*( numeric(1,2)*(1+wbar)/(1-wbar) + Ebar_kernel(0,0,wbar,1).get_numerical_value(pow(qbar,K),N_trunc));
+
+ return ex_to<numeric>(res.evalf());
+ }
+
+ ex pre = pow(2*Pi*I,n_num)/C_norm/K/(n_num-1);
+
+ return get_numerical_value_impl(qbar, pre, 1, N_trunc);
+}
+
+void Kronecker_dtau_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "Kronecker_dtau_kernel(";
+ n.print(c);
+ c.s << ",";
+ z.print(c);
+ c.s << ",";
+ K.print(c);
+ c.s << ",";
+ C_norm.print(c);
+ c.s << ")";
+}
+
+GINAC_BIND_UNARCHIVER(Kronecker_dtau_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Kronecker_dz_kernel, integration_kernel,
+ print_func<print_context>(&Kronecker_dz_kernel::do_print))
+
+Kronecker_dz_kernel::Kronecker_dz_kernel() : inherited(), n(_ex0), z_j(_ex0), tau(_ex0), K(_ex1), C_norm(_ex1)
+{
+}
+
+Kronecker_dz_kernel::Kronecker_dz_kernel(const ex & arg_n, const ex & arg_z_j, const ex & arg_tau, const ex & arg_K, const ex & arg_C_norm) : inherited(), n(arg_n), z_j(arg_z_j), tau(arg_tau), K(arg_K), C_norm(arg_C_norm)
+{
+}
+
+int Kronecker_dz_kernel::compare_same_type(const basic &other) const
+{
+ const Kronecker_dz_kernel &o = static_cast<const Kronecker_dz_kernel &>(other);
+ int cmpval;
+
+ cmpval = n.compare(o.n);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = z_j.compare(o.z_j);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = tau.compare(o.tau);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = K.compare(o.K);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ return C_norm.compare(o.C_norm);
+}
+
+size_t Kronecker_dz_kernel::nops() const
+{
+ return 5;
+}
+
+ex Kronecker_dz_kernel::op(size_t i) const
+{
+ switch (i) {
+ case 0:
+ return n;
+ case 1:
+ return z_j;
+ case 2:
+ return tau;
+ case 3:
+ return K;
+ case 4:
+ return C_norm;
+ default:
+ throw (std::out_of_range("Kronecker_dz_kernel::op() out of range"));
+ }
+}
+
+ex & Kronecker_dz_kernel::let_op(size_t i)
+{
+ ensure_if_modifiable();
+
+ switch (i) {
+ case 0:
+ return n;
+ case 1:
+ return z_j;
+ case 2:
+ return tau;
+ case 3:
+ return K;
+ case 4:
+ return C_norm;
+ default:
+ throw (std::out_of_range("Kronecker_dz_kernel::let_op() out of range"));
+ }
+}
+
+bool Kronecker_dz_kernel::is_numeric(void) const
+{
+ return (n.info(info_flags::nonnegint) && z_j.evalf().info(info_flags::numeric) && tau.evalf().info(info_flags::numeric) && K.info(info_flags::posint) && C_norm.evalf().info(info_flags::numeric));
+}
+
+cln::cl_N Kronecker_dz_kernel::series_coeff_impl(int i) const
+{
+ numeric n_num = ex_to<numeric>(n);
+ int n_int = n_num.to_int();
+
+ ex w_j_inv = exp(ex_to<numeric>((-2*Pi*I*z_j).evalf()));
+ cln::cl_N w_j_inv_cln = ex_to<numeric>(w_j_inv).to_cl_N();
+
+ ex qbar = exp(ex_to<numeric>((2*Pi*I*K*tau).evalf()));
+
+ // case n=0
+ if ( n_num == 0 ) {
+ return 0;
+ }
+
+ // case n=1
+ if ( n_num == 1 ) {
+ if ( i == 1 ) {
+ return ex_to<numeric>((C_norm * 2*Pi*I).evalf()).to_cl_N();
+ }
+
+ return 0;
+ }
+
+ // case n=2
+ if ( n_num == 2 ) {
+ if ( ex_to<numeric>(z_j.evalf()).is_zero() ) {
+ if ( i == 0 ) {
+ return ex_to<numeric>((C_norm).evalf()).to_cl_N();
+ }
+ else if ( i == 1 ) {
+ return 0;
+ }
+ else {
+ ex res = -bernoulli(i)/numeric(i);
+ if ( numeric(i).is_even() ) {
+ Ebar_kernel Ebar = Ebar_kernel( 1-i, 0, numeric(1), numeric(1) );
+ res += Ebar.get_numerical_value(qbar);
+ }
+
+ res *= -pow(2*Pi*I,i)*C_norm/factorial(i-1);
+
+ return ex_to<numeric>(res.evalf()).to_cl_N();
+ }
+ }
+ else {
+ // z_j is not zero
+ if ( i == 0 ) {
+ return 0;
+ }
+ else {
+ Li_negative my_Li_negative;
+
+ ex res = 0;
+ if ( i == 1 ) {
+ res = numeric(1,2);
+ }
+
+ Ebar_kernel Ebar = Ebar_kernel( 1-i, 0, w_j_inv, numeric(1) );
+
+ res += my_Li_negative.get_numerical_value(i-1,w_j_inv) + Ebar.get_numerical_value(qbar);
+
+ res *= -pow(2*Pi*I,i)*C_norm/factorial(i-1);
+
+ return ex_to<numeric>(res.evalf()).to_cl_N();
+ }
+ }
+ }
+
+ // case n>2
+ ex res = 0;
+ if ( i == 1 ) {
+ res += - bernoulli(n_num-1)/(n_num-1);
+ }
+ if ( i > 0 ) {
+ if ( ex_to<numeric>(z_j.evalf()).is_zero() ) {
+ if ( (numeric(i)+n_num).is_even() ) {
+ Ebar_kernel Ebar = Ebar_kernel( 1-i, 2-n_num, numeric(1), numeric(1) );
+
+ res += pow(2*Pi*I,i-1)/factorial(i-1) * Ebar.get_numerical_value(qbar);
+ }
+ }
+ else {
+ // z_j is not zero
+ Ebar_kernel Ebar = Ebar_kernel( 1-i, 2-n_num, w_j_inv, numeric(1) );
+
+ res += pow(2*Pi*I,i-1)/factorial(i-1) * Ebar.get_numerical_value(qbar);
+ }
+ }
+
+ res *= - C_norm * 2*Pi*I/factorial(n_num-2);
+
+ return ex_to<numeric>(res.evalf()).to_cl_N();
+}
+
+/**
+ *
+ * Returns the value of the g^(n-1)(z-z_j,K*tau).
+ *
+ */
+ex Kronecker_dz_kernel::get_numerical_value(const ex & z, int N_trunc) const
+{
+ numeric n_num = ex_to<numeric>(n);
+
+ if ( n_num == 1 ) {
+ return 1;
+ }
+
+ ex pre = pow(2*Pi*I,n-2)/C_norm;
+
+ return get_numerical_value_impl(z, pre, 0, N_trunc);
+}
+
+void Kronecker_dz_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "Kronecker_dz_kernel(";
+ n.print(c);
+ c.s << ",";
+ z_j.print(c);
+ c.s << ",";
+ tau.print(c);
+ c.s << ",";
+ K.print(c);
+ c.s << ",";
+ C_norm.print(c);
+ c.s << ")";
+}
+
+GINAC_BIND_UNARCHIVER(Kronecker_dz_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eisenstein_kernel, integration_kernel,
+ print_func<print_context>(&Eisenstein_kernel::do_print))
+
+Eisenstein_kernel::Eisenstein_kernel() : inherited(), k(_ex0), N(_ex0), a(_ex0), b(_ex0), K(_ex0), C_norm(_ex1)
+{
+}
+
+Eisenstein_kernel::Eisenstein_kernel(const ex & arg_k, const ex & arg_N, const ex & arg_a, const ex & arg_b, const ex & arg_K, const ex & arg_C_norm) : inherited(), k(arg_k), N(arg_N), a(arg_a), b(arg_b), K(arg_K), C_norm(arg_C_norm)
+{
+}
+
+int Eisenstein_kernel::compare_same_type(const basic &other) const
+{
+ const Eisenstein_kernel &o = static_cast<const Eisenstein_kernel &>(other);
+ int cmpval;
+
+ cmpval = k.compare(o.k);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = N.compare(o.N);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = a.compare(o.a);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = b.compare(o.b);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = K.compare(o.K);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ return C_norm.compare(o.C_norm);
+}
+
+/**
+ *
+ * The series method for this class returns the qbar-expansion of the modular form,
+ * without an additional factor of C_norm/qbar.
+ *
+ * This allows for easy use in the class modular_form_kernel.
+ *
+ */
+ex Eisenstein_kernel::series(const relational & r, int order, unsigned options) const
+{
+ if ( r.rhs() != 0 ) {
+ throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
+ }
+
+ ex qbar = r.lhs();
+ ex res = q_expansion_modular_form(qbar, order);
+ res = res.series(qbar,order);
+
+ return res;
+}
+
+size_t Eisenstein_kernel::nops() const
+{
+ return 6;
+}
+
+ex Eisenstein_kernel::op(size_t i) const
+{
+ switch (i) {
+ case 0:
+ return k;
+ case 1:
+ return N;
+ case 2:
+ return a;
+ case 3:
+ return b;
+ case 4:
+ return K;
+ case 5:
+ return C_norm;
+ default:
+ throw (std::out_of_range("Eisenstein_kernel::op() out of range"));
+ }
+}
+
+ex & Eisenstein_kernel::let_op(size_t i)
+{
+ ensure_if_modifiable();
+
+ switch (i) {
+ case 0:
+ return k;
+ case 1:
+ return N;
+ case 2:
+ return a;
+ case 3:
+ return b;
+ case 4:
+ return K;
+ case 5:
+ return C_norm;
+ default:
+ throw (std::out_of_range("Eisenstein_kernel::let_op() out of range"));
+ }
+}
+
+bool Eisenstein_kernel::is_numeric(void) const
+{
+ return (k.info(info_flags::nonnegint) && N.info(info_flags::posint) && a.info(info_flags::integer) && b.info(info_flags::integer) && K.info(info_flags::posint) && C_norm.evalf().info(info_flags::numeric));
+}
+
+ex Eisenstein_kernel::Laurent_series(const ex & x, int order) const
+{
+ ex res = C_norm * q_expansion_modular_form(x, order+1)/x;
+ res = res.series(x,order);
+
+ return res;
+}
+
+/**
+ *
+ * Returns the value of the modular form.
+ *
+ */
+ex Eisenstein_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
+{
+ ex pre = numeric(1)/C_norm;
+
+ return get_numerical_value_impl(qbar, pre, 1, N_trunc);
+}
+
+bool Eisenstein_kernel::uses_Laurent_series() const
+{
+ return true;
+}
+
+ex Eisenstein_kernel::q_expansion_modular_form(const ex & q, int order) const
+{
+ numeric k_num = ex_to<numeric>(k);
+ numeric N_num = ex_to<numeric>(N);
+ numeric a_num = ex_to<numeric>(a);
+ numeric b_num = ex_to<numeric>(b);
+ numeric K_num = ex_to<numeric>(K);
+
+ if ( (k==2) && (a==1) && (b==1) ) {
+ return B_eisenstein_series(q, N_num, K_num, order);
+ }
+
+ return E_eisenstein_series(q, k_num, N_num, a_num, b_num, K_num, order);
+}
+
+void Eisenstein_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "Eisenstein_kernel(";
+ k.print(c);
+ c.s << ",";
+ N.print(c);
+ c.s << ",";
+ a.print(c);
+ c.s << ",";
+ b.print(c);
+ c.s << ",";
+ K.print(c);
+ c.s << ",";
+ C_norm.print(c);
+ c.s << ")";
+}
+
+GINAC_BIND_UNARCHIVER(Eisenstein_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eisenstein_h_kernel, integration_kernel,
+ print_func<print_context>(&Eisenstein_h_kernel::do_print))
+
+Eisenstein_h_kernel::Eisenstein_h_kernel() : inherited(), k(_ex0), N(_ex0), r(_ex0), s(_ex0), C_norm(_ex1)
+{
+}
+
+Eisenstein_h_kernel::Eisenstein_h_kernel(const ex & arg_k, const ex & arg_N, const ex & arg_r, const ex & arg_s, const ex & arg_C_norm) : inherited(), k(arg_k), N(arg_N), r(arg_r), s(arg_s), C_norm(arg_C_norm)
+{
+}
+
+int Eisenstein_h_kernel::compare_same_type(const basic &other) const
+{
+ const Eisenstein_h_kernel &o = static_cast<const Eisenstein_h_kernel &>(other);
+ int cmpval;
+
+ cmpval = k.compare(o.k);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = N.compare(o.N);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = r.compare(o.r);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = s.compare(o.s);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ return C_norm.compare(o.C_norm);
+}
+
+/**
+ *
+ * The series method for this class returns the qbar-expansion of the modular form,
+ * without an additional factor of C_norm/qbar.
+ *
+ * This allows for easy use in the class modular_form_kernel.
+ *
+ */
+ex Eisenstein_h_kernel::series(const relational & r, int order, unsigned options) const
+{
+ if ( r.rhs() != 0 ) {
+ throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
+ }
+
+ ex qbar = r.lhs();
+ ex res = q_expansion_modular_form(qbar, order);
+ res = res.series(qbar,order);
+
+ return res;
+}
+
+size_t Eisenstein_h_kernel::nops() const
+{
+ return 5;
+}
+
+ex Eisenstein_h_kernel::op(size_t i) const
+{
+ switch (i) {
+ case 0:
+ return k;
+ case 1:
+ return N;
+ case 2:
+ return r;
+ case 3:
+ return s;
+ case 4:
+ return C_norm;
+ default:
+ throw (std::out_of_range("Eisenstein_h_kernel::op() out of range"));
+ }
+}
+
+ex & Eisenstein_h_kernel::let_op(size_t i)
+{
+ ensure_if_modifiable();
+
+ switch (i) {
+ case 0:
+ return k;
+ case 1:
+ return N;
+ case 2:
+ return r;
+ case 3:
+ return s;
+ case 4:
+ return C_norm;
+ default:
+ throw (std::out_of_range("Eisenstein_h_kernel::let_op() out of range"));
+ }
+}
+
+bool Eisenstein_h_kernel::is_numeric(void) const
+{
+ return (k.info(info_flags::nonnegint) && N.info(info_flags::posint) && r.info(info_flags::integer) && s.info(info_flags::integer) && C_norm.evalf().info(info_flags::numeric));
+}
+
+ex Eisenstein_h_kernel::Laurent_series(const ex & x, int order) const
+{
+ ex res = C_norm * q_expansion_modular_form(x, order+1)/x;
+ res = res.series(x,order);
+
+ return res;
+}
+
+/**
+ *
+ * Returns the value of the modular form.
+ *
+ */
+ex Eisenstein_h_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
+{
+ ex pre = numeric(1)/C_norm;
+
+ return get_numerical_value_impl(qbar, pre, 1, N_trunc);
+}
+
+bool Eisenstein_h_kernel::uses_Laurent_series() const
+{
+ return true;
+}
+
+/**
+ *
+ * The constant coefficient in the Fourier expansion.
+ *
+ */
+ex Eisenstein_h_kernel::coefficient_a0(const numeric & k, const numeric & r, const numeric & s, const numeric & N) const
+{
+ if ( k == 1 ) {
+ if ( irem(s,N) != 0 ) {
+ return numeric(1,4) - mod(s,N)/numeric(2)/N;
+ }
+ else if ( (irem(r,N)==0) && (irem(s,N)==0) ) {
+ return 0;
+ }
+ else {
+ return I*numeric(1,4)*cos(Pi*mod(r,N)/N)/sin(Pi*mod(r,N)/N);
+ }
+ }
+
+ // case k > 1
+ return -Bernoulli_polynomial(k,mod(s,N)/N)/numeric(2)/k;
+}
+
+/**
+ *
+ * The higher coefficients in the Fourier expansion.
+ *
+ */
+ex Eisenstein_h_kernel::coefficient_an(const numeric & n, const numeric & k, const numeric & r, const numeric & s, const numeric & N) const
+{
+ ex res = 0;
+
+ for (numeric m=1; m<=n; m++) {
+ if ( irem(n,m) == 0 ) {
+ for (numeric c1=0; c1<N; c1++) {
+ numeric c2 = n/m;
+ res += pow(m,k-1)*exp(2*Pi*I/N*mod(r*c2-(s-m)*c1,N)) - pow(-m,k-1)*exp(2*Pi*I/N*mod(-r*c2+(s+m)*c1,N));
+ }
+ }
+ }
+
+ return res/numeric(2)/pow(N,k);
+}
+
+ex Eisenstein_h_kernel::q_expansion_modular_form(const ex & q, int N_order) const
+{
+ numeric N_order_num = numeric(N_order);
+
+ numeric k_num = ex_to<numeric>(k);
+ numeric r_num = ex_to<numeric>(r);
+ numeric s_num = ex_to<numeric>(s);
+ numeric N_num = ex_to<numeric>(N);
+
+ ex res = coefficient_a0(k_num,r_num,s_num,N_num);
+
+ for (numeric i1=1; i1<N_order_num; i1++) {
+ res += coefficient_an(i1,k_num,r_num,s_num,N_num) * pow(q,i1);
+ }
+
+ res += Order(pow(q,N_order));
+ res = res.series(q,N_order);
+
+ return res;
+}
+
+void Eisenstein_h_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "Eisenstein_h_kernel(";
+ k.print(c);
+ c.s << ",";
+ N.print(c);
+ c.s << ",";
+ r.print(c);
+ c.s << ",";
+ s.print(c);
+ c.s << ",";
+ C_norm.print(c);
+ c.s << ")";
+}
+
+GINAC_BIND_UNARCHIVER(Eisenstein_h_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(modular_form_kernel, integration_kernel,
+ print_func<print_context>(&modular_form_kernel::do_print))
+
+modular_form_kernel::modular_form_kernel() : inherited(), k(_ex0), P(_ex0), C_norm(_ex1)
+{
+}
+
+modular_form_kernel::modular_form_kernel(const ex & arg_k, const ex & arg_P, const ex & arg_C_norm) : inherited(), k(arg_k), P(arg_P), C_norm(arg_C_norm)
+{
+}
+
+int modular_form_kernel::compare_same_type(const basic &other) const
+{
+ const modular_form_kernel &o = static_cast<const modular_form_kernel &>(other);
+ int cmpval;
+
+ cmpval = k.compare(o.k);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ cmpval = P.compare(o.P);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ return C_norm.compare(o.C_norm);
+}
+
+/**
+ *
+ * The series method for this class returns the qbar-expansion of the modular form,
+ * without an additional factor of C_norm/qbar.
+ *
+ */
+ex modular_form_kernel::series(const relational & r, int order, unsigned options) const
+{
+ if ( r.rhs() != 0 ) {
+ throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
+ }
+
+ ex qbar = r.lhs();
+
+ subs_q_expansion do_subs_q_expansion(qbar, order);
+
+ ex res = do_subs_q_expansion(P).series(qbar,order);
+
+ return res;
+}
+
+size_t modular_form_kernel::nops() const
+{
+ return 3;
+}
+
+ex modular_form_kernel::op(size_t i) const
+{
+ switch (i) {
+ case 0:
+ return k;
+ case 1:
+ return P;
+ case 2:
+ return C_norm;
+ default:
+ throw (std::out_of_range("modular_form_kernel::op() out of range"));
+ }
+}
+
+ex & modular_form_kernel::let_op(size_t i)
+{
+ ensure_if_modifiable();
+
+ switch (i) {
+ case 0:
+ return k;
+ case 1:
+ return P;
+ case 2:
+ return C_norm;
+ default:
+ throw (std::out_of_range("modular_form_kernel::let_op() out of range"));
+ }
+}
+
+bool modular_form_kernel::is_numeric(void) const
+{
+ bool flag = (k.info(info_flags::nonnegint) && C_norm.evalf().info(info_flags::numeric));
+ if ( !flag ) {
+ return false;
+ }
+
+ symbol qbar("qbar");
+
+ // test with a random number and random expansion
+ return series_to_poly(P.series(qbar,18)).subs(qbar==numeric(1,937)).evalf().info(info_flags::numeric);
+}
+
+ex modular_form_kernel::Laurent_series(const ex & qbar, int order) const
+{
+ ex res = series_to_poly(P.series(qbar,order+1));
+ res = C_norm * res/qbar;
+ res = res.series(qbar,order);
+ return res;
+}
+
+/**
+ *
+ * Returns the value of the modular form.
+ *
+ */
+ex modular_form_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
+{
+ ex pre = numeric(1)/C_norm;
+
+ return get_numerical_value_impl(qbar, pre, 1, N_trunc);
+}
+
+bool modular_form_kernel::uses_Laurent_series() const
+{
+ return true;
+}
+
+ex modular_form_kernel::q_expansion_modular_form(const ex & q, int N_order) const
+{
+ return this->series(q==0,N_order);
+}
+
+void modular_form_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "modular_form_kernel(";
+ k.print(c);
+ c.s << ",";
+ P.print(c);
+ c.s << ",";
+ C_norm.print(c);
+ c.s << ")";
+}
+
+GINAC_BIND_UNARCHIVER(modular_form_kernel);
+
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(user_defined_kernel, integration_kernel,
+ print_func<print_context>(&user_defined_kernel::do_print))
+
+user_defined_kernel::user_defined_kernel() : inherited(), f(_ex0), x(_ex0)
+{
+}
+
+user_defined_kernel::user_defined_kernel(const ex & arg_f, const ex & arg_x) : inherited(), f(arg_f), x(arg_x)
+{
+}
+
+int user_defined_kernel::compare_same_type(const basic &other) const
+{
+ const user_defined_kernel &o = static_cast<const user_defined_kernel &>(other);
+ int cmpval;
+
+ cmpval = f.compare(o.f);
+ if ( cmpval) {
+ return cmpval;
+ }
+
+ return x.compare(o.x);
+}
+
+size_t user_defined_kernel::nops() const
+{
+ return 2;
+}
+
+ex user_defined_kernel::op(size_t i) const
+{
+ switch (i) {
+ case 0:
+ return f;
+ case 1:
+ return x;
+ default:
+ throw (std::out_of_range("user_defined_kernel::op() out of range"));
+ }
+}
+
+ex & user_defined_kernel::let_op(size_t i)
+{
+ ensure_if_modifiable();
+
+ switch (i) {
+ case 0:
+ return f;
+ case 1:
+ return x;
+ default:
+ throw (std::out_of_range("user_defined_kernel::let_op() out of range"));
+ }
+}
+
+bool user_defined_kernel::is_numeric(void) const
+{
+ // test with a random number
+ return f.subs(x==numeric(1,937)).evalf().info(info_flags::numeric);
+}
+
+ex user_defined_kernel::Laurent_series(const ex & x_up, int order) const
+{
+ ex res = f.series(x,order).subs(x==x_up);
+
+ return res;
+}
+
+bool user_defined_kernel::uses_Laurent_series() const
+{
+ return true;
+}
+
+void user_defined_kernel::do_print(const print_context & c, unsigned level) const
+{
+ c.s << "user_defined_kernel(";
+ f.print(c);
+ c.s << ",";
+ x.print(c);
+ c.s << ")";
+}
+
+GINAC_BIND_UNARCHIVER(user_defined_kernel);
+
+} // namespace GiNaC
--- /dev/null
+/** @file integration_kernel.h
+ *
+ * Interface to GiNaC's integration kernels for iterated integrals. */
+
+/*
+ * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifndef GINAC_INTEGRATION_KERNEL_H
+#define GINAC_INTEGRATION_KERNEL_H
+
+#include "basic.h"
+#include "archive.h"
+#include "numeric.h"
+
+#include <cln/complex.h>
+#include <vector>
+
+namespace GiNaC {
+
+ex ifactor(const numeric & n);
+bool is_discriminant_of_quadratic_number_field(const numeric & n);
+numeric kronecker_symbol(const numeric & a, const numeric & n);
+numeric primitive_dirichlet_character(const numeric & n, const numeric & a);
+numeric dirichlet_character(const numeric & n, const numeric & a, const numeric & N);
+numeric generalised_Bernoulli_number(const numeric & k, const numeric & b);
+ex Bernoulli_polynomial(const numeric & k, const ex & x);
+
+/**
+ *
+ * The base class for integration kernels for iterated integrals.
+ *
+ * This class represents the differential one-form
+ * \f[
+ * \omega = d\lambda
+ * \f]
+ * The integration variable is a dummy variable and does not need to be specified.
+ *
+ */
+class integration_kernel : public basic
+{
+ GINAC_DECLARE_REGISTERED_CLASS(integration_kernel, basic)
+
+ // ctors
+public:
+
+ // functions overriding virtual functions from base classes
+public:
+ ex series(const relational & r, int order, unsigned options = 0) const override;
+
+protected:
+
+ // new virtual functions which can be overridden by derived classes
+public:
+ virtual bool has_trailing_zero(void) const;
+ virtual bool is_numeric(void) const;
+ virtual ex Laurent_series(const ex & x, int order) const;
+ virtual ex get_numerical_value(const ex & lambda, int N_trunc = 0) const;
+
+protected:
+ virtual bool uses_Laurent_series() const;
+ virtual cln::cl_N series_coeff_impl(int i) const;
+
+ // non-virtual functions
+public:
+ size_t get_cache_size(void) const;
+ void set_cache_step(int cache_steps) const;
+ ex get_series_coeff(int i) const;
+ cln::cl_N series_coeff(int i) const;
+
+protected:
+ ex get_numerical_value_impl(const ex & lambda, const ex & pre, int shift, int N_trunc) const;
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ // cache is increased by steps of cache_step_size
+ mutable int cache_step_size;
+ // cache already computed series coefficients
+ mutable std::vector<cln::cl_N> series_vec;
+
+};
+
+GINAC_DECLARE_UNARCHIVER(integration_kernel);
+
+/**
+ *
+ * The basic integration kernel with a logarithmic singularity at the origin.
+ *
+ * This class represents the differential one-form
+ * \f[
+ * L_0 = \frac{d\lambda}{\lambda}
+ * \f]
+ *
+ */
+class basic_log_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(basic_log_kernel, integration_kernel)
+
+ // ctors
+public:
+
+ // functions overriding virtual functions from base classes
+public:
+
+protected:
+ cln::cl_N series_coeff_impl(int i) const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+
+};
+
+GINAC_DECLARE_UNARCHIVER(basic_log_kernel);
+
+/**
+ *
+ * The integration kernel for multiple polylogarithms.
+ *
+ * This class represents the differential one-form
+ * \f[
+ * \omega^{\mathrm{mpl}}(z) = \frac{d\lambda}{\lambda-z}
+ * \f]
+ *
+ * For the case \f$ z=0 \f$ the class basic_log_kernel should be used.
+ *
+ */
+class multiple_polylog_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(multiple_polylog_kernel, integration_kernel)
+
+ // ctors
+public:
+ multiple_polylog_kernel(const ex & z);
+
+ // functions overriding virtual functions from base classes
+public:
+ size_t nops() const override;
+ ex op(size_t i) const override;
+ ex & let_op(size_t i) override;
+
+ bool is_numeric(void) const override;
+
+protected:
+ cln::cl_N series_coeff_impl(int i) const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ ex z;
+
+};
+
+GINAC_DECLARE_UNARCHIVER(multiple_polylog_kernel);
+
+/**
+ *
+ * The ELi-kernel.
+ *
+ * This class represents the differential one-form
+ * \f[
+ * \omega^{\mathrm{ELi}}_{n;m}(x;y) = \mathrm{ELi}_{n;m}(x;y;\bar{q}) \frac{d\bar{q}}{\bar{q}}
+ * \f]
+ *
+ */
+class ELi_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(ELi_kernel, integration_kernel)
+
+ // ctors
+public:
+ ELi_kernel(const ex & n, const ex & m, const ex & x, const ex & y);
+
+ // functions overriding virtual functions from base classes
+public:
+ size_t nops() const override;
+ ex op(size_t i) const override;
+ ex & let_op(size_t i) override;
+
+ bool is_numeric(void) const override;
+ ex get_numerical_value(const ex & qbar, int N_trunc = 0) const override;
+
+protected:
+ cln::cl_N series_coeff_impl(int i) const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ ex n;
+ ex m;
+ ex x;
+ ex y;
+
+};
+
+GINAC_DECLARE_UNARCHIVER(ELi_kernel);
+
+/**
+ *
+ * The Ebar-kernel
+ *
+ * This class represents the differential one-form
+ * \f[
+ * \omega^{\overline{\mathrm{E}}}_{n;m}(x;y) = \overline{\mathrm{E}}_{n;m}(x;y;\bar{q}) \frac{d\bar{q}}{\bar{q}}
+ * \f]
+ *
+ */
+class Ebar_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(Ebar_kernel, integration_kernel)
+
+ // ctors
+public:
+ Ebar_kernel(const ex & n, const ex & m, const ex & x, const ex & y);
+
+ // functions overriding virtual functions from base classes
+public:
+ size_t nops() const override;
+ ex op(size_t i) const override;
+ ex & let_op(size_t i) override;
+
+ bool is_numeric(void) const override;
+ ex get_numerical_value(const ex & qbar, int N_trunc = 0) const override;
+
+protected:
+ cln::cl_N series_coeff_impl(int i) const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ ex n;
+ ex m;
+ ex x;
+ ex y;
+
+};
+
+GINAC_DECLARE_UNARCHIVER(Ebar_kernel);
+
+/**
+ *
+ * The kernel corresponding to integrating the Kronecker coefficient function \f$ g^{(n)}(z_j,K \tau) \f$
+ * in \f$ \tau \f$ (or equivalently in \f$ \bar{q} \f$).
+ *
+ * This class represents the differential one-form
+ * \f[
+ * \omega^{\mathrm{Kronecker},\tau}_{n,K}(z_j) = \frac{C_n K (n-1)}{(2\pi i)^n} g^{(n)}(z_j,K \tau) \frac{d\bar{q}}{\bar{q}}
+ * \f]
+ *
+ */
+class Kronecker_dtau_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(Kronecker_dtau_kernel, integration_kernel)
+
+ // ctors
+public:
+ Kronecker_dtau_kernel(const ex & n, const ex & z, const ex & K = numeric(1), const ex & C_norm = numeric(1));
+
+ // functions overriding virtual functions from base classes
+public:
+ size_t nops() const override;
+ ex op(size_t i) const override;
+ ex & let_op(size_t i) override;
+
+ bool is_numeric(void) const override;
+ ex get_numerical_value(const ex & qbar, int N_trunc = 0) const override;
+
+protected:
+ cln::cl_N series_coeff_impl(int i) const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ ex n;
+ ex z;
+ ex K;
+ ex C_norm;
+};
+
+GINAC_DECLARE_UNARCHIVER(Kronecker_dtau_kernel);
+
+
+/**
+ *
+ * The kernel corresponding to integrating the Kronecker coefficient function \f$ g^{(n-1)}(z-z_j, K \tau) \f$
+ * in \f$ z \f$.
+ *
+ * This class represents the differential one-form
+ * \f[
+ * \omega^{\mathrm{Kronecker},z}_{n,K}(z_j,\tau) = C_n (2\pi i)^{2-n} g^{(n-1)}(z-z_j, K \tau) dz
+ * \f]
+ *
+ */
+class Kronecker_dz_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(Kronecker_dz_kernel, integration_kernel)
+
+ // ctors
+public:
+ Kronecker_dz_kernel(const ex & n, const ex & z_j, const ex & tau, const ex & K = numeric(1), const ex & C_norm = numeric(1));
+
+ // functions overriding virtual functions from base classes
+public:
+ size_t nops() const override;
+ ex op(size_t i) const override;
+ ex & let_op(size_t i) override;
+
+ bool is_numeric(void) const override;
+ ex get_numerical_value(const ex & z, int N_trunc = 0) const override;
+
+protected:
+ cln::cl_N series_coeff_impl(int i) const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ ex n;
+ ex z_j;
+ ex tau;
+ ex K;
+ ex C_norm;
+
+};
+
+GINAC_DECLARE_UNARCHIVER(Kronecker_dz_kernel);
+
+
+/**
+ *
+ * The kernel corresponding to the Eisenstein series \f$ E_{k,N,a,b,K}(\tau) \f$.
+ *
+ * This class represents the differential one-form
+ * \f[
+ * \omega^{\mathrm{Eisenstein}}_{k,N,a,b,K} = C_k E_{k,N,a,b,K}(\tau) \frac{d\bar{q}_N}{\bar{q}_N}
+ * \f]
+ *
+ * The integers a and b are either one or the discriminant of a quadratic number field.
+ * This class represents Eisenstein series, which can be defined by primitive Dirichlet characters from the Kronecker symbol.
+ * This implies that the characters take the values -1,0,1, i.e. no higher roots of unity occur.
+ * The \f[ \bar{q} \f]-expansion has then rational coefficients.
+ *
+ * Ref.: W. Stein, Modular Forms: A Computational Approach, Chapter 5
+ *
+ */
+class Eisenstein_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(Eisenstein_kernel, integration_kernel)
+
+ // ctors
+public:
+ Eisenstein_kernel(const ex & k, const ex & N, const ex & a, const ex & b, const ex & K, const ex & C_norm = numeric(1));
+
+ // functions overriding virtual functions from base classes
+public:
+ ex series(const relational & r, int order, unsigned options = 0) const override;
+
+ size_t nops() const override;
+ ex op(size_t i) const override;
+ ex & let_op(size_t i) override;
+
+ bool is_numeric(void) const override;
+ ex Laurent_series(const ex & x, int order) const override;
+ ex get_numerical_value(const ex & qbar, int N_trunc = 0) const override;
+
+protected:
+ bool uses_Laurent_series() const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+ ex q_expansion_modular_form(const ex & q, int order) const;
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ ex k;
+ ex N;
+ ex a;
+ ex b;
+ ex K;
+ ex C_norm;
+
+};
+
+GINAC_DECLARE_UNARCHIVER(Eisenstein_kernel);
+
+
+/**
+ *
+ * The kernel corresponding to the Eisenstein series \f$ h_{k,N,r,s}(\tau) \f$.
+ *
+ * This class represents the differential one-form
+ * \f[
+ * \omega^{\mathrm{Eisenstein,h}}_{k,N,r,s} = C_k h_{k,N,r,s}(\tau) \frac{d\bar{q}_N}{\bar{q}_N}
+ * \f]
+ *
+ */
+class Eisenstein_h_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(Eisenstein_h_kernel, integration_kernel)
+
+ // ctors
+public:
+ Eisenstein_h_kernel(const ex & k, const ex & N, const ex & r, const ex & s, const ex & C_norm = numeric(1));
+
+ // functions overriding virtual functions from base classes
+public:
+ ex series(const relational & r, int order, unsigned options = 0) const override;
+
+ size_t nops() const override;
+ ex op(size_t i) const override;
+ ex & let_op(size_t i) override;
+
+ bool is_numeric(void) const override;
+ ex Laurent_series(const ex & x, int order) const override;
+ ex get_numerical_value(const ex & qbar, int N_trunc = 0) const override;
+
+protected:
+ bool uses_Laurent_series() const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+ ex coefficient_a0(const numeric & k, const numeric & r, const numeric & s, const numeric & N) const;
+ ex coefficient_an(const numeric & n, const numeric & k, const numeric & r, const numeric & s, const numeric & N) const;
+ ex q_expansion_modular_form(const ex & q, int order) const;
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ ex k;
+ ex N;
+ ex r;
+ ex s;
+ ex C_norm;
+
+};
+
+GINAC_DECLARE_UNARCHIVER(Eisenstein_h_kernel);
+
+
+/**
+ *
+ * A kernel corresponding to a polynomial in Eisenstein series.
+ *
+ * This class represents the differential one-form
+ * \f[
+ * \omega^{\mathrm{modular}}(P_k(\eta^{(1)}_{k_1}, \dots, \eta^{(r)}_{k_r})) = C_k P_k(\eta^{(1)}_{k_1}, \dots, \eta^{(r)}_{k_r}) \frac{d\bar{q}_N}{\bar{q}_N}.
+ * \f]
+ *
+ */
+class modular_form_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(modular_form_kernel, integration_kernel)
+
+ // ctors
+public:
+ modular_form_kernel(const ex & k, const ex & P, const ex & C_norm = numeric(1));
+
+ // functions overriding virtual functions from base classes
+public:
+ ex series(const relational & r, int order, unsigned options = 0) const override;
+
+ size_t nops() const override;
+ ex op(size_t i) const override;
+ ex & let_op(size_t i) override;
+
+ bool is_numeric(void) const override;
+ ex Laurent_series(const ex & qbar, int order) const override;
+ ex get_numerical_value(const ex & qbar, int N_trunc = 0) const override;
+
+protected:
+ bool uses_Laurent_series() const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+ ex q_expansion_modular_form(const ex & q, int order) const;
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ ex k;
+ ex P;
+ ex C_norm;
+
+};
+
+GINAC_DECLARE_UNARCHIVER(modular_form_kernel);
+
+
+/**
+ *
+ * A user-defined integration kernel.
+ * The input is an expression \f$ f \f$, depending on a variable \f$ x \f$.
+ * It is assumed that \f$ f \f$ has a Laurent expansion around \f$ x=0 \f$ and
+ * maximally a simple pole at \f$ x=0 \f$.
+ *
+ */
+class user_defined_kernel : public integration_kernel
+{
+ GINAC_DECLARE_REGISTERED_CLASS(user_defined_kernel, integration_kernel)
+
+ // ctors
+public:
+ user_defined_kernel(const ex & f, const ex & x);
+
+ // functions overriding virtual functions from base classes
+public:
+ size_t nops() const override;
+ ex op(size_t i) const override;
+ ex & let_op(size_t i) override;
+
+ bool is_numeric(void) const override;
+ ex Laurent_series(const ex & x, int order) const override;
+
+protected:
+ bool uses_Laurent_series() const override;
+
+ // new virtual functions which can be overridden by derived classes
+public:
+
+protected:
+
+ // non-virtual functions
+public:
+
+protected:
+ void do_print(const print_context & c, unsigned level) const;
+
+ // friends :
+
+ // member variables :
+
+protected:
+ ex f;
+ ex x;
+
+};
+
+GINAC_DECLARE_UNARCHIVER(user_defined_kernel);
+
+} // namespace GiNaC
+
+#endif // ndef GINAC_INTEGRATION_KERNEL_H
--- /dev/null
+/** @file utils_multi_iterator.h
+ *
+ * Utilities for summing over multiple indices */
+
+/*
+ * GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
+ *
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
+ */
+
+#ifndef GINAC_UTILS_MULTI_ITERATOR_H
+#define GINAC_UTILS_MULTI_ITERATOR_H
+
+#include <cstddef>
+#include <vector>
+#include <ostream>
+#include <iterator>
+
+namespace GiNaC {
+
+/**
+ *
+ * SFINAE test for distance
+ *
+ */
+template <typename T> class has_distance {
+private:
+ typedef char yes_type[1];
+ typedef char no_type[2];
+
+ template <typename C> static yes_type & test( decltype(std::distance<C>) ) ;
+ template <typename C> static no_type & test(...);
+
+public:
+ enum { value = sizeof(test<T>(0)) == sizeof(yes_type) };
+};
+
+/**
+ *
+ * For printing a multi-index:
+ * If the templates are used, where T is an iterator, printing the address where the iterator points to is not meaningful.
+ * However, we may print the difference to the starting point.
+ *
+ */
+template<typename T> typename std::enable_if<has_distance<T>::value, typename std::iterator_traits<T>::difference_type>::type format_index_value(const T & a, const T & b) {
+ return std::distance(a,b);
+}
+
+/**
+ *
+ * For all other cases we simply print the value.
+ *
+ */
+template<typename T> typename std::enable_if<!has_distance<T>::value, T>::type format_index_value(const T & a, const T & b) {
+ return b;
+}
+
+/**
+ *
+ * basic_multi_iterator is a base class.
+ *
+ * The base class itself does not do anything useful.
+ * A typical use of a class derived from basic_multi_iterator is
+ *
+ * multi_iterator_ordered<int> k(0,4,2);
+ *
+ * for( k.init(); !k.overflow(); k++) {
+ * std::cout << k << std::endl;
+ * }
+ *
+ * which prints out
+ *
+ * multi_iterator_ordered(0,1)
+ * multi_iterator_ordered(0,2)
+ * multi_iterator_ordered(0,3)
+ * multi_iterator_ordered(1,2)
+ * multi_iterator_ordered(1,3)
+ * multi_iterator_ordered(2,3)
+ *
+ * Individual components of k can be accessed with k[i] or k(i).
+ *
+ * All classes derived from basic_multi_iterator follow the same syntax.
+ *
+ */
+template<class T> class basic_multi_iterator {
+
+ // ctors
+public :
+ basic_multi_iterator(void);
+ explicit basic_multi_iterator(T B, T N, size_t k);
+ explicit basic_multi_iterator(T B, T N, const std::vector<T> & vv);
+
+ // dtor
+ virtual ~basic_multi_iterator();
+
+ // functions
+public :
+ size_t size(void) const;
+ bool overflow(void) const;
+ const std::vector<T> & get_vector(void) const;
+
+ // subscripting
+public :
+ T operator[](size_t i) const;
+ T & operator[](size_t i);
+
+ T operator()(size_t i) const;
+ T & operator()(size_t i);
+
+ // virtual functions
+public :
+ // initialization
+ virtual basic_multi_iterator<T> & init(void);
+ // postfix increment
+ virtual basic_multi_iterator<T> & operator++ (int);
+
+ // I/O operators
+ template <class TT> friend std::ostream & operator<< (std::ostream & os, const basic_multi_iterator<TT> & v);
+
+ // member variables :
+protected :
+ T N;
+ T B;
+ std::vector<T> v;
+ bool flag_overflow;
+
+};
+
+/**
+ *
+ * The class multi_iterator_ordered defines a multi_iterator
+ * \f$(i_1,i_2,...,i_k)\f$, such that
+ * \f[
+ * B \le i_j < N
+ * \f]
+ * and
+ * \f[
+ * i_j < i_{j+1}.
+ * \f]
+ * It is assumed that \f$k>0\f$ and \f$ N-B \ge k \f$.
+ *
+ */
+template<class T> class multi_iterator_ordered : public basic_multi_iterator<T> {
+
+ // ctors
+public :
+ multi_iterator_ordered(void);
+ explicit multi_iterator_ordered(T B, T N, size_t k);
+ explicit multi_iterator_ordered(T B, T N, const std::vector<T> & vv);
+
+ // overriding virtual functions from base class
+public :
+ // initialization
+ basic_multi_iterator<T> & init(void);
+ // postfix increment
+ basic_multi_iterator<T> & operator++ (int);
+
+ // I/O operators
+ template <class TT> friend std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered<TT> & v);
+
+};
+
+/**
+ *
+ * The class multi_iterator_ordered_eq defines a multi_iterator
+ * \f$(i_1,i_2,...,i_k)\f$, such that
+ * \f[
+ * B \le i_j < N
+ * \f]
+ * and
+ * \f[
+ * i_j \le i_{j+1}.
+ * \f]
+ * It is assumed that \f$k>0\f$.
+ *
+ */
+template<class T> class multi_iterator_ordered_eq : public basic_multi_iterator<T> {
+
+ // ctors
+public :
+ multi_iterator_ordered_eq(void);
+ explicit multi_iterator_ordered_eq(T B, T N, size_t k);
+ explicit multi_iterator_ordered_eq(T B, T N, const std::vector<T> & vv);
+
+ // overriding virtual functions from base class
+public :
+ // initialization
+ basic_multi_iterator<T> & init(void);
+ // postfix increment
+ basic_multi_iterator<T> & operator++ (int);
+
+ // I/O operators
+ template <class TT> friend std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered_eq<TT> & v);
+
+};
+
+/**
+ *
+ * The class multi_iterator_ordered_eq_indv defines a multi_iterator
+ * \f$(i_1,i_2,...,i_k)\f$, such that
+ * \f[
+ * B \le i_j < N_j
+ * \f]
+ * and
+ * \f[
+ * i_j \le i_{j+1}.
+ * \f]
+ *
+ */
+template<class T> class multi_iterator_ordered_eq_indv : public basic_multi_iterator<T> {
+
+ // ctors
+public :
+ multi_iterator_ordered_eq_indv(void);
+ explicit multi_iterator_ordered_eq_indv(T B, const std::vector<T> & Nv, size_t k);
+ explicit multi_iterator_ordered_eq_indv(T B, const std::vector<T> & Nv, const std::vector<T> & vv);
+
+ // overriding virtual functions from base class
+public :
+ // initialization
+ basic_multi_iterator<T> & init(void);
+ // postfix increment
+ basic_multi_iterator<T> & operator++ (int);
+
+ // I/O operators
+ template <class TT> friend std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered_eq_indv<TT> & v);
+
+ // member variables :
+protected :
+ std::vector<T> Nv;
+};
+
+/**
+ *
+ * The class multi_iterator_counter defines a multi_iterator
+ * \f$(i_1,i_2,...,i_k)\f$, such that
+ * \f[
+ * B \le i_j < N
+ * \f]
+ *
+ */
+template<class T> class multi_iterator_counter : public basic_multi_iterator<T> {
+
+ // ctors
+public :
+ multi_iterator_counter(void);
+ explicit multi_iterator_counter(T B, T N, size_t k);
+ explicit multi_iterator_counter(T B, T N, const std::vector<T> & vv);
+
+ // overriding virtual functions from base class
+public :
+ // initialization
+ basic_multi_iterator<T> & init(void);
+ // postfix increment
+ basic_multi_iterator<T> & operator++ (int);
+
+ // I/O operators
+ template <class TT> friend std::ostream & operator<< (std::ostream & os, const multi_iterator_counter<TT> & v);
+
+};
+
+/**
+ *
+ * The class multi_iterator_counter_indv defines a multi_iterator
+ * \f$(i_1,i_2,...,i_k)\f$, such that
+ * \f[
+ * B \le i_j < N_j
+ * \f]
+ *
+ */
+template<class T> class multi_iterator_counter_indv : public basic_multi_iterator<T> {
+
+ // ctors
+public :
+ multi_iterator_counter_indv(void);
+ explicit multi_iterator_counter_indv(T B, const std::vector<T> & Nv, size_t k);
+ explicit multi_iterator_counter_indv(T B, const std::vector<T> & Nv, const std::vector<T> & vv);
+
+ // overriding virtual functions from base class
+public :
+ // initialization
+ basic_multi_iterator<T> & init(void);
+ // postfix increment
+ basic_multi_iterator<T> & operator++ (int);
+
+ // I/O operators
+ template <class TT> friend std::ostream & operator<< (std::ostream & os, const multi_iterator_counter_indv<TT> & v);
+
+ // member variables :
+protected :
+ std::vector<T> Nv;
+};
+
+/**
+ *
+ * The class multi_iterator_permutation defines a multi_iterator
+ * \f$(i_1,i_2,...,i_k)\f$, for which
+ * \f[
+ * B \le i_j < N
+ * \f]
+ * and
+ * \f[
+ * i_i \neq i_j
+ * \f]
+ * In particular, if \f$N-B=k\f$, multi_iterator_permutation loops over all
+ * permutations of \f$k\f$ elements.
+ *
+ */
+template<class T> class multi_iterator_permutation : public basic_multi_iterator<T> {
+
+ // ctors
+public :
+ multi_iterator_permutation(void);
+ explicit multi_iterator_permutation(T B, T N, size_t k);
+ explicit multi_iterator_permutation(T B, T N, const std::vector<T> & vv);
+
+ // overriding virtual functions from base class
+public :
+ // initialization
+ basic_multi_iterator<T> & init(void);
+ // postfix increment
+ basic_multi_iterator<T> & operator++ (int);
+
+ // new functions in this class
+ int get_sign(void) const;
+
+ // I/O operators
+ template <class TT> friend std::ostream & operator<< (std::ostream & os, const multi_iterator_permutation<TT> & v);
+
+};
+
+/**
+ *
+ * The class multi_iterator_shuffle defines a multi_iterator,
+ * which runs over all shuffles of a and b.
+ *
+ */
+template<class T> class multi_iterator_shuffle : public basic_multi_iterator<T> {
+
+ // ctors
+public :
+ multi_iterator_shuffle(void);
+ explicit multi_iterator_shuffle(const std::vector<T> & a, const std::vector<T> & b);
+
+ // overriding virtual functions from base class
+public :
+ // initialization
+ basic_multi_iterator<T> & init(void);
+ // postfix increment
+ basic_multi_iterator<T> & operator++ (int);
+
+ // I/O operators
+ template <class TT> friend std::ostream & operator<< (std::ostream & os, const multi_iterator_shuffle<TT> & v);
+
+ // member variables :
+protected :
+ size_t N_internal;
+ std::vector<size_t> v_internal;
+ std::vector<T> v_orig;
+};
+
+/**
+ *
+ * The class multi_iterator_shuffle_prime defines a multi_iterator,
+ * which runs over all shuffles of a and b, excluding the first one (a,b).
+ *
+ */
+template<class T> class multi_iterator_shuffle_prime : public multi_iterator_shuffle<T> {
+
+ // ctors
+public :
+ multi_iterator_shuffle_prime(void);
+ explicit multi_iterator_shuffle_prime(const std::vector<T> & a, const std::vector<T> & b);
+
+ // overriding virtual functions from base class
+public :
+ // initialization
+ basic_multi_iterator<T> & init(void);
+
+ // I/O operators
+ template <class TT> friend std::ostream & operator<< (std::ostream & os, const multi_iterator_shuffle_prime<TT> & v);
+};
+
+// ----------------------------------------------------------------------------------------------------------------
+
+// ctors
+
+/**
+ *
+ * Default constructor
+ *
+ */
+template<class T> inline basic_multi_iterator<T>::basic_multi_iterator(void) : N(), B(), v(), flag_overflow(false)
+{}
+
+/**
+ *
+ * Construct a multi_iterator with upper limit N, lower limit B and size k .
+ *
+ */
+template<class T> inline basic_multi_iterator<T>::basic_multi_iterator(T BB, T NN, size_t k) : N(NN), B(BB), v(k), flag_overflow(false)
+{}
+
+/**
+ *
+ * Construct from a vector.
+ *
+ */
+template<class T> inline basic_multi_iterator<T>::basic_multi_iterator(T BB, T NN, const std::vector<T> & vv) : N(NN), B(BB), v(vv), flag_overflow(false)
+{}
+
+/**
+ *
+ * Destructor
+ *
+ */
+template<class T> inline basic_multi_iterator<T>::~basic_multi_iterator()
+{}
+
+// functions
+
+/**
+ *
+ * Returns the size of a multi_iterator.
+ *
+ */
+template<class T> inline size_t basic_multi_iterator<T>::size(void) const
+{
+ return v.size();
+}
+
+/**
+ *
+ * Initialize the multi-index to
+ * \f[
+ * (n_1,n_2,n_3,...,n_k) = (B,B,...,B)
+ * \f]
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & basic_multi_iterator<T>::init(void)
+{
+ flag_overflow = false;
+
+ for ( size_t i=0; i<v.size(); i++) {
+ v[i] = B;
+ }
+ return *this;
+}
+
+/**
+ *
+ * Return the overflow flag.
+ *
+ */
+template<class T> inline bool basic_multi_iterator<T>::overflow(void) const
+{
+ return flag_overflow;
+}
+
+/**
+ *
+ * Returns a reference to the vector v.
+ *
+ */
+template<class T> inline const std::vector<T> & basic_multi_iterator<T>::get_vector(void) const
+{
+ return v;
+}
+
+// subscripting
+
+/**
+ *
+ * Subscription via []
+ *
+ */
+template<class T> inline T basic_multi_iterator<T>::operator[](size_t i) const
+{
+ return v[i];
+}
+
+/**
+ *
+ * Subscription via []
+ *
+ */
+template<class T> inline T & basic_multi_iterator<T>::operator[](size_t i)
+{
+ return v[i];
+}
+
+/**
+ *
+ * Subscription via ()
+ *
+ */
+template<class T> inline T basic_multi_iterator<T>::operator()(size_t i) const
+{
+ return v[i];
+}
+
+/**
+ *
+ * Subscription via ()
+ *
+ */
+template<class T> inline T & basic_multi_iterator<T>::operator()(size_t i)
+{
+ return v[i];
+}
+
+
+/**
+ *
+ * No effect for basic_multi_iterator
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & basic_multi_iterator<T>::operator++ (int)
+{
+ return *this;
+}
+
+// I/O operators
+
+/**
+ *
+ * Output operator. A multi_iterator prints out as
+ * basic_multi_iterator(\f$n_0,n_1,...\f$).
+ *
+ */
+template<class T> inline std::ostream & operator<< (std::ostream & os, const basic_multi_iterator<T> & v)
+{
+ os << "basic_multi_iterator(";
+ for ( size_t i=0; i<v.size(); i++) {
+ if (i>0) {
+ os << ",";
+ }
+ os << format_index_value(v.B,v(i));
+ }
+
+ return os << ")";
+}
+
+
+
+// ctors
+
+/**
+ *
+ * Default constructor
+ *
+ */
+template<class T> inline multi_iterator_ordered<T>::multi_iterator_ordered(void) : basic_multi_iterator<T>()
+{}
+
+/**
+ *
+ * Construct a multi_iterator with upper limit N and size k .
+ *
+ */
+template<class T> inline multi_iterator_ordered<T>::multi_iterator_ordered(T B, T N, size_t k) : basic_multi_iterator<T>(B,N,k)
+{}
+
+/**
+ *
+ * Construct from a vector.
+ *
+ */
+template<class T> inline multi_iterator_ordered<T>::multi_iterator_ordered(T B, T N, const std::vector<T> & v) : basic_multi_iterator<T>(B,N,v)
+{}
+
+// functions
+
+/**
+ *
+ * Initialize the multi-index to
+ * \f[
+ * (n_1,n_2,n_3,...,n_k) = (B+0,B+1,B+2,...,B+k-1)
+ * \f]
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_ordered<T>::init(void)
+{
+ this->flag_overflow = false;
+ T it = this->B;
+
+ for ( size_t i=0; i < this->v.size(); i++) {
+ this->v[i] = it;
+ it++;
+ }
+ return *this;
+}
+
+/**
+ *
+ * The postfix increment operator allows to
+ * write for a multi-index n++, which will
+ * update n to the next configuration.
+ *
+ * If n is in the last configuration and the
+ * increment operator ++ is applied to n,
+ * the overflow flag will be raised.
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_ordered<T>::operator++ (int)
+{
+ int k = this->size();
+ int j = k - 1;
+ T Upper_limit = this->N;
+
+ while ( j>0 ) {
+ this->v[j]++;
+ if ( this->v[j] == Upper_limit ) {
+ j--;
+ Upper_limit--;
+ }
+ else {
+ break;
+ }
+ }
+
+ if (j==0) {
+ this->v[j]++;
+ if (this->v[j] == Upper_limit) this->flag_overflow=true;
+ }
+
+ if ( j>= 0) {
+ for (int jj=j+1;jj<k;jj++) {
+ this->v[jj] = this->v[jj-1];
+ this->v[jj]++;
+ }
+ }
+
+ return *this;
+}
+
+// I/O operators
+
+/**
+ *
+ * Output operator. A multi_iterator_ordered prints out as
+ * multi_iterator_ordered(\f$n_0,n_1,...\f$).
+ *
+ */
+template<class T> inline std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered<T> & v)
+{
+ os << "multi_iterator_ordered(";
+ for ( size_t i=0; i<v.size(); i++) {
+ if (i>0) {
+ os << ",";
+ }
+ os << format_index_value(v.B,v(i));
+ }
+
+ return os << ")";
+}
+
+
+
+// ctors
+
+/**
+ *
+ * Default constructor
+ *
+ */
+template<class T> inline multi_iterator_ordered_eq<T>::multi_iterator_ordered_eq(void) : basic_multi_iterator<T>()
+{}
+
+/**
+ *
+ * Construct a multi_iterator with upper limit N and size k .
+ *
+ */
+template<class T> inline multi_iterator_ordered_eq<T>::multi_iterator_ordered_eq(T B, T N, size_t k) : basic_multi_iterator<T>(B,N,k)
+{}
+
+/**
+ *
+ * Construct from a vector.
+ *
+ */
+template<class T> inline multi_iterator_ordered_eq<T>::multi_iterator_ordered_eq(T B, T N, const std::vector<T> & v) : basic_multi_iterator<T>(B,N,v)
+{}
+
+// functions
+
+/**
+ *
+ * Initialize the multi-index to
+ * \f[
+ * (n_1,n_2,...,n_k) = (B,B,...,B)
+ * \f]
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_ordered_eq<T>::init(void)
+{
+ this->flag_overflow = false;
+
+ for ( size_t i=0; i < this->v.size(); i++) {
+ this->v[i] = this->B;
+ }
+ return *this;
+}
+
+/**
+ *
+ * The postfix increment operator allows to
+ * write for a multi-index n++, which will
+ * update n to the next configuration.
+ *
+ * If n is in the last configuration and the
+ * increment operator ++ is applied to n,
+ * the overflow flag will be raised.
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_ordered_eq<T>::operator++ (int)
+{
+ int k = this->size();
+ int j = k - 1;
+
+ while ( j>0 ) {
+ this->v[j]++;
+ if ( this->v[j] == this->N ) {
+ j--;
+ }
+ else {
+ break;
+ }
+ }
+
+ if (j==0) {
+ this->v[j]++;
+ if (this->v[j] == this->N) {
+ this->flag_overflow=true;
+ }
+ }
+
+ if ( j>= 0) {
+ for (int jj=j+1;jj<k;jj++) {
+ this->v[jj] = this->v[jj-1];
+ }
+ }
+
+ return *this;
+}
+
+// I/O operators
+
+/**
+ *
+ * Output operator. A multi_iterator_ordered_eq prints out as
+ * multi_iterator_ordered_eq(\f$n_0,n_1,...\f$).
+ *
+ */
+template<class T> inline std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered_eq<T> & v)
+{
+ os << "multi_iterator_ordered_eq(";
+ for ( size_t i=0; i<v.size(); i++) {
+ if (i>0) {
+ os << ",";
+ }
+ os << format_index_value(v.B,v(i));
+ }
+
+ return os << ")";
+}
+
+
+
+
+// ctors
+
+/**
+ *
+ * Default constructor
+ *
+ */
+template<class T> inline multi_iterator_ordered_eq_indv<T>::multi_iterator_ordered_eq_indv(void) : basic_multi_iterator<T>(), Nv()
+{}
+
+/**
+ *
+ * Construct a multi_iterator with upper limit N and size k .
+ *
+ */
+template<class T> inline multi_iterator_ordered_eq_indv<T>::multi_iterator_ordered_eq_indv(T B, const std::vector<T> & Nvv, size_t k) : basic_multi_iterator<T>(B,B,k), Nv(Nvv)
+{}
+
+/**
+ *
+ * Construct from a vector.
+ *
+ */
+template<class T> inline multi_iterator_ordered_eq_indv<T>::multi_iterator_ordered_eq_indv(T B, const std::vector<T> & Nvv, const std::vector<T> & v) : basic_multi_iterator<T>(B,B,v), Nv(Nvv)
+{}
+
+// functions
+
+/**
+ *
+ * Initialize the multi-index to
+ * \f[
+ * (n_1,n_2,n_3,...,n_k) = (B,B,B,...,B)
+ * \f]
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_ordered_eq_indv<T>::init(void)
+{
+ this->flag_overflow = false;
+
+ for ( size_t i=0; i < this->v.size(); i++) {
+ this->v[i] = this->B;
+ }
+ return *this;
+}
+
+/**
+ *
+ * The postfix increment operator allows to
+ * write for a multi-index n++, which will
+ * update n to the next configuration.
+ *
+ * If n is in the last configuration and the
+ * increment operator ++ is applied to n,
+ * the overflow flag will be raised.
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_ordered_eq_indv<T>::operator++ (int)
+{
+ int k = this->size();
+ int j = k - 1;
+
+ while ( j>0 ) {
+ this->v[j]++;
+ if ( this->v[j] == Nv[j] ) {
+ j--;
+ }
+ else {
+ break;
+ }
+ }
+
+ if (j==0) {
+ this->v[j]++;
+ if (this->v[j] == Nv[j]) {
+ this->flag_overflow=true;
+ }
+ }
+
+ if ( j>= 0) {
+ for (int jj=j+1;jj<k;jj++) {
+ this->v[jj] = this->v[jj-1];
+ }
+ }
+
+ return *this;
+}
+
+// I/O operators
+
+/**
+ *
+ * Output operator. A multi_iterator_ordered_eq_indv prints out as
+ * multi_iterator_ordered_eq_indv(\f$n_0,n_1,...\f$).
+ *
+ */
+template<class T> inline std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered_eq_indv<T> & v)
+{
+ os << "multi_iterator_ordered_eq_indv(";
+ for ( size_t i=0; i<v.size(); i++) {
+ if (i>0) {
+ os << ",";
+ }
+ os << format_index_value(v.B,v(i));
+ }
+
+ return os << ")";
+}
+
+
+
+
+// ctors
+
+/**
+ *
+ * Default constructor
+ *
+ */
+template<class T> inline multi_iterator_counter<T>::multi_iterator_counter(void) : basic_multi_iterator<T>()
+{}
+
+/**
+ *
+ * Construct a multi_iterator with upper limit N and size k .
+ *
+ */
+template<class T> inline multi_iterator_counter<T>::multi_iterator_counter(T B, T N, size_t k) : basic_multi_iterator<T>(B,N,k)
+{}
+
+/**
+ *
+ * Construct from a vector.
+ *
+ */
+template<class T> inline multi_iterator_counter<T>::multi_iterator_counter(T B, T N, const std::vector<T> & v) : basic_multi_iterator<T>(B,N,v)
+{}
+
+// functions
+
+/**
+ *
+ * Initialize the multi-index to
+ * \f[
+ * (n_1,n_2,n_3,...,n_k) = (B,B,...,B)
+ * \f]
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_counter<T>::init(void)
+{
+ this->flag_overflow = false;
+
+ for ( size_t i=0; i < this->v.size(); i++) {
+ this->v[i] = this->B;
+ }
+ return *this;
+}
+
+/**
+ *
+ * The postfix increment operator allows to
+ * write for a multi-index n++, which will
+ * update n to the next configuration.
+ *
+ * If n is in the last configuration and the
+ * increment operator ++ is applied to n,
+ * the overflow flag will be raised.
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_counter<T>::operator++ (int)
+{
+ int k = this->size();
+ int j = k - 1;
+
+ while ( j>0 ) {
+ this->v[j]++;
+ if ( this->v[j] == this->N ) {
+ this->v[j] = this->B;
+ j--;
+ }
+ else {
+ break;
+ }
+ }
+
+ if (j==0) {
+ this->v[j]++;
+ if (this->v[j] == this->N) {
+ this->v[j] = this->B;
+ this->flag_overflow=true;
+ }
+ }
+
+ return *this;
+}
+
+// I/O operators
+
+/**
+ *
+ * Output operator. A multi_iterator_counter prints out as
+ * multi_iterator_counter(\f$n_0,n_1,...\f$).
+ *
+ */
+template<class T> inline std::ostream & operator<< (std::ostream & os, const multi_iterator_counter<T> & v)
+{
+ os << "multi_iterator_counter(";
+ for ( size_t i=0; i<v.size(); i++) {
+ if (i>0) {
+ os << ",";
+ }
+ os << format_index_value(v.B,v(i));
+ }
+
+ return os << ")";
+}
+
+
+
+
+// ctors
+
+/**
+ *
+ * Default constructor
+ *
+ */
+template<class T> inline multi_iterator_counter_indv<T>::multi_iterator_counter_indv(void) : basic_multi_iterator<T>(), Nv()
+{}
+
+/**
+ *
+ * Construct a multi_iterator with upper limit N and size k .
+ *
+ */
+template<class T> inline multi_iterator_counter_indv<T>::multi_iterator_counter_indv(T B, const std::vector<T> & Nvv, size_t k) : basic_multi_iterator<T>(B,B,k), Nv(Nvv)
+{}
+
+/**
+ *
+ * Construct from a vector.
+ *
+ */
+template<class T> inline multi_iterator_counter_indv<T>::multi_iterator_counter_indv(T B, const std::vector<T> & Nvv, const std::vector<T> & v) : basic_multi_iterator<T>(B,B,v), Nv(Nvv)
+{}
+
+// functions
+
+/**
+ *
+ * Initialize the multi-index to
+ * \f[
+ * (n_1,n_2,n_3,...,n_k) = (B,B,...,B)
+ * \f]
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_counter_indv<T>::init(void)
+{
+ this->flag_overflow = false;
+
+ for ( size_t i=0; i < this->v.size(); i++) {
+ this->v[i] = this->B;
+ }
+ return *this;
+}
+
+/**
+ *
+ * The postfix increment operator allows to
+ * write for a multi-index n++, which will
+ * update n to the next configuration.
+ *
+ * If n is in the last configuration and the
+ * increment operator ++ is applied to n,
+ * the overflow flag will be raised.
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_counter_indv<T>::operator++ (int)
+{
+ int k = this->size();
+ int j = k - 1;
+
+ while ( j>0 ) {
+ this->v[j]++;
+ if ( this->v[j] == Nv[j] ) {
+ this->v[j] = this->B;
+ j--;
+ }
+ else {
+ break;
+ }
+ }
+
+ if (j==0) {
+ this->v[j]++;
+ if (this->v[j] == Nv[j]) {
+ this->v[j] = this->B;
+ this->flag_overflow=true;
+ }
+ }
+
+ return *this;
+}
+
+// I/O operators
+
+/**
+ *
+ * Output operator. A multi_iterator_counter_indv prints out as
+ * multi_iterator_counter_indv(\f$n_0,n_1,...\f$).
+ *
+ */
+template<class T> inline std::ostream & operator<< (std::ostream & os, const multi_iterator_counter_indv<T> & v)
+{
+ os << "multi_iterator_counter_indv(";
+ for ( size_t i=0; i<v.size(); i++) {
+ if (i>0) {
+ os << ",";
+ }
+ os << format_index_value(v.B,v(i));
+ }
+
+ return os << ")";
+}
+
+
+
+
+// ctors
+
+/**
+ *
+ * Default constructor
+ *
+ */
+template<class T> inline multi_iterator_permutation<T>::multi_iterator_permutation(void) : basic_multi_iterator<T>()
+{}
+
+/**
+ *
+ * Construct a multi_iterator with upper limit N and size k .
+ *
+ */
+template<class T> inline multi_iterator_permutation<T>::multi_iterator_permutation(T B, T N, size_t k) : basic_multi_iterator<T>(B,N,k)
+{}
+
+/**
+ *
+ * Construct from a vector.
+ *
+ */
+template<class T> inline multi_iterator_permutation<T>::multi_iterator_permutation(T B, T N, const std::vector<T> & v) : basic_multi_iterator<T>(B,N,v)
+{}
+
+// functions
+
+/**
+ *
+ * Initialize the multi-index to
+ * \f[
+ * (n_1,n_2,n_3,...,n_k) = (B+0,B+1,B+2,...,B+k-1)
+ * \f]
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_permutation<T>::init(void)
+{
+ this->flag_overflow = false;
+ T it = this->B;
+
+ for ( size_t i=0; i < this->v.size(); i++) {
+ this->v[i] = it;
+ it++;
+ }
+ return *this;
+}
+
+/**
+ *
+ * The postfix increment operator allows to
+ * write for a multi-index n++, which will
+ * update n to the next configuration.
+ *
+ * If n is in the last configuration and the
+ * increment operator ++ is applied to n,
+ * the overflow flag will be raised.
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_permutation<T>::operator++ (int)
+{
+ int k = this->size();
+ int j = k - 1;
+
+ while ( j>=0 ) {
+ bool flag_have_already = true;
+ while ( flag_have_already ) {
+ this->v[j]++;
+
+ // update flag_have_already
+ flag_have_already = false;
+ for (int ii=0; ii<j; ii++) {
+ if (this->v[j] == this->v[ii]) {
+ flag_have_already = true;
+ }
+ }
+ }
+
+ if ( this->v[j] == this->N ) {
+ j--;
+ }
+ else {
+ break;
+ }
+ }
+
+ for (int l=j+1; l<k; l++) {
+ this->v[l] = this->B;
+
+ bool flag_have_already;
+ do {
+ flag_have_already = false;
+ for (int ii=0; ii<l; ii++) {
+ if (this->v[l] == this->v[ii]) {
+ flag_have_already = true;
+ }
+ }
+ if (flag_have_already) {
+ this->v[l]++;
+ }
+ }
+ while (flag_have_already);
+ }
+
+ // check for overflow
+ this->flag_overflow = true;
+ T it = this->B;
+ for (int ii=0; ii<k; ii++) {
+ if (this->v[ii] != it) {
+ this->flag_overflow = false;
+ }
+ it++;
+ }
+
+ return *this;
+}
+
+/**
+ *
+ * Returns the sign of the permutation, defined by
+ * \f[
+ * \left(-1\right)^{n_{inv}},
+ * \f]
+ * where \f$ n_{inv} \f$ is the number of inversions, e.g. the
+ * number of pairs \f$ i < j \f$ for which
+ * \f[
+ * n_i > n_j.
+ * \f]
+ *
+ */
+template<class T> inline int multi_iterator_permutation<T>::get_sign() const
+{
+ int sign = 1;
+ int k = this->size();
+
+ for ( int i=0; i<k; i++) {
+ for ( int j=i+1; j<k; j++) {
+ // works only for random-access iterators
+ if ( this->v[i] > this->v[j] ) {
+ sign = -sign;
+ }
+ }
+ }
+
+ return sign;
+}
+
+
+// I/O operators
+
+/**
+ *
+ * Output operator. A multi_iterator_permutation prints out as
+ * multi_iterator_permutation(\f$n_0,n_1,...\f$).
+ *
+ */
+template<class T> inline std::ostream & operator<< (std::ostream & os, const multi_iterator_permutation<T> & v)
+{
+ os << "multi_iterator_permutation(";
+ for ( size_t i=0; i<v.size(); i++) {
+ if (i>0) {
+ os << ",";
+ }
+ os << format_index_value(v.B,v(i));
+ }
+
+ return os << ")";
+}
+
+
+
+// ctors
+
+/**
+ *
+ * Default constructor
+ *
+ */
+template<class T> inline multi_iterator_shuffle<T>::multi_iterator_shuffle(void) : basic_multi_iterator<T>(), N_internal(), v_internal(), v_orig()
+{}
+
+/**
+ *
+ * Construct from a vector.
+ *
+ */
+template<class T> inline multi_iterator_shuffle<T>::multi_iterator_shuffle(const std::vector<T> & a, const std::vector<T> & b) : basic_multi_iterator<T>(), N_internal(), v_internal(), v_orig()
+{
+ this->B = a[0];
+
+ for (size_t i=0; i<a.size(); i++) {
+ this->v.push_back( a[i] );
+ this->v_orig.push_back( a[i] );
+ this->v_internal.push_back( i );
+ }
+ for (size_t i=0; i<b.size(); i++) {
+ this->v.push_back( b[i] );
+ this->v_orig.push_back( b[i] );
+ }
+ this->N_internal = this->v.size();
+}
+
+// functions
+
+/**
+ *
+ * Initialize the multi-index to the first shuffle.
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_shuffle<T>::init(void)
+{
+ this->flag_overflow = false;
+
+ for ( size_t i=0; i < this->v_internal.size(); i++) {
+ this->v_internal[i] = i;
+ }
+ for ( size_t i=0; i < this->v.size(); i++) {
+ this->v[i] = this->v_orig[i];
+ }
+ return *this;
+}
+
+/**
+ *
+ * The postfix increment operator allows to
+ * write for a multi-index n++, which will
+ * update n to the next configuration.
+ *
+ * If n is in the last configuration and the
+ * increment operator ++ is applied to n,
+ * the overflow flag will be raised.
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_shuffle<T>::operator++ (int)
+{
+ int k = this->v_internal.size();
+ int j = k - 1;
+ size_t Upper_limit = this->N_internal;
+
+ while ( j>0 ) {
+ this->v_internal[j]++;
+ if ( this->v_internal[j] == Upper_limit ) {
+ j--;
+ Upper_limit--;
+ }
+ else {
+ break;
+ }
+ }
+
+ if (j==0) {
+ this->v_internal[j]++;
+ if (this->v_internal[j] == Upper_limit) {
+ this->flag_overflow=true;
+ }
+ }
+
+ if ( j>= 0) {
+ for (int jj=j+1;jj<k;jj++) {
+ this->v_internal[jj] = this->v_internal[jj-1];
+ this->v_internal[jj]++;
+ }
+ }
+
+ // update v
+ if ( !(this->flag_overflow) ) {
+ size_t i_a = 0;
+ size_t i_b = 0;
+ size_t i_all = 0;
+ for (size_t j=0; j<k; j++) {
+ for (size_t i=i_all; i < this->v_internal[j]; i++) {
+ this->v[i_all] = this->v_orig[k+i_b];
+ i_b++;
+ i_all++;
+ }
+ this->v[i_all] = this->v_orig[i_a];
+ i_a++;
+ i_all++;
+ }
+ for (size_t i = this->v_internal[k-1]+1; i < this->v.size(); i++) {
+ this->v[i_all] = this->v_orig[k+i_b];
+ i_b++;
+ i_all++;
+ }
+ }
+
+ return *this;
+}
+
+// I/O operators
+
+/**
+ *
+ * Output operator. A multi_iterator_shuffle prints out as
+ * multi_iterator_shuffle(\f$n_0,n_1,...\f$).
+ *
+ */
+template<class T> inline std::ostream & operator<< (std::ostream & os, const multi_iterator_shuffle<T> & v)
+{
+ os << "multi_iterator_shuffle(";
+ for ( size_t i=0; i<v.size(); i++) {
+ if (i>0) {
+ os << ",";
+ }
+ os << format_index_value(v.B,v(i));
+ }
+
+ return os << ")";
+}
+
+
+
+// ctors
+
+/**
+ *
+ * Default constructor
+ *
+ */
+template<class T> inline multi_iterator_shuffle_prime<T>::multi_iterator_shuffle_prime(void) : multi_iterator_shuffle<T>()
+{}
+
+/**
+ *
+ * Construct from a vector.
+ *
+ */
+template<class T> inline multi_iterator_shuffle_prime<T>::multi_iterator_shuffle_prime(const std::vector<T> & a, const std::vector<T> & b) : multi_iterator_shuffle<T>(a,b)
+{}
+
+// functions
+
+/**
+ *
+ * Initialize the multi-index to the first shuffle.
+ *
+ */
+template<class T> inline basic_multi_iterator<T> & multi_iterator_shuffle_prime<T>::init(void)
+{
+ this->flag_overflow = false;
+
+ for ( size_t i=0; i < this->v_internal.size(); i++) {
+ this->v_internal[i] = i;
+ }
+ for ( size_t i=0; i < this->v.size(); i++) {
+ this->v[i] = this->v_orig[i];
+ }
+
+ (*this)++;
+
+ return *this;
+}
+
+// I/O operators
+
+/**
+ *
+ * Output operator. A multi_iterator_shuffle_prime prints out as
+ * multi_iterator_shuffle_prime(\f$n_0,n_1,...\f$).
+ *
+ */
+template<class T> inline std::ostream & operator<< (std::ostream & os, const multi_iterator_shuffle_prime<T> & v)
+{
+ os << "multi_iterator_shuffle_prime(";
+ for ( size_t i=0; i<v.size(); i++) {
+ if (i>0) {
+ os << ",";
+ }
+ os << format_index_value(v.B,v(i));
+ }
+
+ return os << ")";
+}
+
+} // namespace GiNaC
+
+#endif // ndef GINAC_UTILS_MULTI_ITERATOR_H
return e[0].unit(e[1]);
}
+static ex f_basic_log_kernel(const exprseq &e)
+{
+ return basic_log_kernel();
+}
+
+static ex f_multiple_polylog_kernel(const exprseq &e)
+{
+ return multiple_polylog_kernel(e[0]);
+}
+
+static ex f_ELi_kernel(const exprseq &e)
+{
+ return ELi_kernel(e[0],e[1],e[2],e[3]);
+}
+
+static ex f_Ebar_kernel(const exprseq &e)
+{
+ return Ebar_kernel(e[0],e[1],e[2],e[3]);
+}
+
+static ex f_Kronecker_dtau_kernel_4(const exprseq &e)
+{
+ return Kronecker_dtau_kernel(e[0],e[1],e[2],e[3]);
+}
+
+static ex f_Kronecker_dtau_kernel_3(const exprseq &e)
+{
+ return Kronecker_dtau_kernel(e[0],e[1],e[2]);
+}
+
+static ex f_Kronecker_dtau_kernel_2(const exprseq &e)
+{
+ return Kronecker_dtau_kernel(e[0],e[1]);
+}
+
+static ex f_Kronecker_dz_kernel_5(const exprseq &e)
+{
+ return Kronecker_dz_kernel(e[0],e[1],e[2],e[3],e[4]);
+}
+
+static ex f_Kronecker_dz_kernel_4(const exprseq &e)
+{
+ return Kronecker_dz_kernel(e[0],e[1],e[2],e[3]);
+}
+
+static ex f_Kronecker_dz_kernel_3(const exprseq &e)
+{
+ return Kronecker_dz_kernel(e[0],e[1],e[2]);
+}
+
+static ex f_Eisenstein_kernel_6(const exprseq &e)
+{
+ return Eisenstein_kernel(e[0],e[1],e[2],e[3],e[4],e[5]);
+}
+
+static ex f_Eisenstein_kernel_5(const exprseq &e)
+{
+ return Eisenstein_kernel(e[0],e[1],e[2],e[3],e[4]);
+}
+
+static ex f_Eisenstein_h_kernel_5(const exprseq &e)
+{
+ return Eisenstein_h_kernel(e[0],e[1],e[2],e[3],e[4]);
+}
+
+static ex f_Eisenstein_h_kernel_4(const exprseq &e)
+{
+ return Eisenstein_h_kernel(e[0],e[1],e[2],e[3]);
+}
+
+static ex f_modular_form_kernel_3(const exprseq &e)
+{
+ return modular_form_kernel(e[0],e[1],e[2]);
+}
+
+static ex f_modular_form_kernel_2(const exprseq &e)
+{
+ return modular_form_kernel(e[0],e[1]);
+}
+
+static ex f_user_defined_kernel(const exprseq &e)
+{
+ return user_defined_kernel(e[0],e[1]);
+}
+
+static ex f_q_expansion_modular_form(const exprseq &e)
+{
+ if ( is_a<Eisenstein_kernel>(e[0]) ) {
+ return ex_to<Eisenstein_kernel>(e[0]).q_expansion_modular_form(e[1], ex_to<numeric>(e[2]).to_int());
+ }
+ if ( is_a<Eisenstein_h_kernel>(e[0]) ) {
+ return ex_to<Eisenstein_h_kernel>(e[0]).q_expansion_modular_form(e[1], ex_to<numeric>(e[2]).to_int());
+ }
+ if ( is_a<modular_form_kernel>(e[0]) ) {
+ return ex_to<modular_form_kernel>(e[0]).q_expansion_modular_form(e[1], ex_to<numeric>(e[2]).to_int());
+ }
+ throw(std::invalid_argument("first argument must be a modular form"));
+}
+
static ex f_dummy(const exprseq &e)
{
throw(std::logic_error("dummy function called (shouldn't happen)"));
{"transpose", f_transpose, 1},
{"unassign", f_unassign, 1},
{"unit", f_unit, 2},
+ {"basic_log_kernel", f_basic_log_kernel, 0},
+ {"multiple_polylog_kernel", f_multiple_polylog_kernel, 1},
+ {"ELi_kernel", f_ELi_kernel, 4},
+ {"Ebar_kernel", f_Ebar_kernel, 4},
+ {"Kronecker_dtau_kernel", f_Kronecker_dtau_kernel_4, 4},
+ {"Kronecker_dtau_kernel", f_Kronecker_dtau_kernel_3, 3},
+ {"Kronecker_dtau_kernel", f_Kronecker_dtau_kernel_2, 2},
+ {"Kronecker_dz_kernel", f_Kronecker_dz_kernel_5, 5},
+ {"Kronecker_dz_kernel", f_Kronecker_dz_kernel_4, 4},
+ {"Kronecker_dz_kernel", f_Kronecker_dz_kernel_3, 3},
+ {"Eisenstein_kernel", f_Eisenstein_kernel_6, 6},
+ {"Eisenstein_kernel", f_Eisenstein_kernel_5, 5},
+ {"Eisenstein_h_kernel", f_Eisenstein_h_kernel_5, 5},
+ {"Eisenstein_h_kernel", f_Eisenstein_h_kernel_4, 4},
+ {"modular_form_kernel", f_modular_form_kernel_3, 3},
+ {"modular_form_kernel", f_modular_form_kernel_2, 2},
+ {"user_defined_kernel", f_user_defined_kernel, 2},
+ {"q_expansion_modular_form", f_q_expansion_modular_form, 3},
{nullptr, f_dummy, 0} // End marker
};
{"tan", "tangent function"},
{"tanh", "hyperbolic tangent function"},
{"zeta", "zeta function\nzeta(x) is Riemann's zeta function, zetaderiv(n,x) its nth derivative.\nIf x is a GiNaC::lst, it is a multiple zeta value\nzeta(x,s) is an alternating Euler sum"},
+ {"G", "multiple polylogarithm (integral representation)"},
{"Li2", "dilogarithm"},
{"Li3", "trilogarithm"},
{"Li", "(multiple) polylogarithm"},
{"S", "Nielsen's generalized polylogarithm"},
{"H", "harmonic polylogarithm"},
+ {"EllipticK", "complete elliptic integral of the first kind"},
+ {"EllipticE", "complete elliptic integral of the second kind"},
+ {"iterated_integral", "iterated integral"},
{"Order", "order term function (for truncated power series)"},
{"Derivative", "inert differential operator"},
{nullptr, nullptr} // End marker
git://www.ginac.de / ginac.git / commitdiff