Added new routines for the numerical evaluation of iterated integrals like
authorStefan Weinzierl <weinzierl@uni-mainz.de>
Sat, 10 Oct 2020 15:11:09 +0000 (17:11 +0200)
committerStefan Weinzierl <weinzierl@uni-mainz.de>
Sat, 10 Oct 2020 15:11:09 +0000 (17:11 +0200)
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

13 files changed:
check/CMakeLists.txt
check/Makefile.am
check/exam_inifcns_elliptic.cpp [new file with mode: 0644]
doc/tutorial/ginac.texi
ginac/CMakeLists.txt
ginac/Makefile.am
ginac/ginac.h
ginac/inifcns.h
ginac/inifcns_elliptic.cpp [new file with mode: 0644]
ginac/integration_kernel.cpp [new file with mode: 0644]
ginac/integration_kernel.h [new file with mode: 0644]
ginac/utils_multi_iterator.h [new file with mode: 0644]
ginsh/ginsh_parser.ypp

index caabbfc3240beb27d45e3b8f784091f673537bc6..6b8bd340ac48747fd6ca71f846f52537da7c7b32 100644 (file)
@@ -14,6 +14,7 @@ set(ginac_tests
        exam_powerlaws
        exam_inifcns
        exam_inifcns_nstdsums
+       exam_inifcns_elliptic
        exam_differentiation
        exam_polygcd
        exam_collect_common_factors
index 1fffd9f77ef1df9710bdc15c76992e650ddddb3b..d1881520f2702efb73bad022c7ab61ab43fe98ab 100644 (file)
@@ -14,6 +14,7 @@ EXAMS = exam_paranoia \
        exam_powerlaws  \
        exam_inifcns \
        exam_inifcns_nstdsums \
+       exam_inifcns_elliptic \
        exam_differentiation  \
        exam_polygcd  \
        exam_collect_common_factors \
@@ -105,6 +106,9 @@ exam_inifcns_nstdsums_SOURCES = exam_inifcns_nstdsums.cpp \
                                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
 
diff --git a/check/exam_inifcns_elliptic.cpp b/check/exam_inifcns_elliptic.cpp
new file mode 100644 (file)
index 0000000..7539a92
--- /dev/null
@@ -0,0 +1,419 @@
+/** @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();
+}
index 1ba42151e86a39898c9c2dbef89cc4268c1d5a3f..520a13dfb1f47c303c5014fd32de518dbecbc7a5 100644 (file)
@@ -3826,6 +3826,7 @@ avoided.
 * 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.
@@ -5873,6 +5874,12 @@ GiNaC contains the following predefined mathematical functions:
 @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()}
@@ -5888,6 +5895,12 @@ GiNaC contains the following predefined mathematical functions:
 @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()}
@@ -5968,7 +5981,7 @@ GiNaC uses the opposite order: firstly expands the function and then its
 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
 
@@ -6151,7 +6164,71 @@ J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001)
 @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
index 086665511c974f92d12a1151b0cf745b57ee4f04..5b0937c4089a93b3831201ae073a7838d60fe73a 100644 (file)
@@ -20,7 +20,9 @@ set(ginaclib_sources
     inifcns.cpp
     inifcns_gamma.cpp
     inifcns_nstdsums.cpp
+    inifcns_elliptic.cpp
     inifcns_trans.cpp
+    integration_kernel.cpp
     integral.cpp
     lst.cpp
     matrix.cpp
@@ -88,6 +90,7 @@ set(ginaclib_public_headers
     indexed.h 
     inifcns.h
     integral.h
+    integration_kernel.h
     lst.h
     matrix.h
     mul.h
@@ -116,6 +119,7 @@ set(ginaclib_private_headers
     utils.h
     crc32.h
     hash_seed.h
+    utils_multi_iterator.h
     parser/lexer.h
     parser/debug.h
     polynomial/gcd_euclid.h
index cb76517aa89e017c198232323b4dd4ccf2e1d84f..0dc82c1503d6e253340701330eaffd4c58190a8a 100644 (file)
@@ -4,12 +4,13 @@ lib_LTLIBRARIES = libginac.la
 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 \
@@ -64,7 +65,7 @@ ginacincludedir = $(includedir)/ginac
 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 \
index 84a26f478184fdeb207b6f081ce50d76247f391a..f3905dfefee97a53e42a8605d1dfea45334d9872 100644 (file)
@@ -68,6 +68,8 @@
 
 #include "factor.h"
 
+#include "integration_kernel.h"
+
 #include "excompiler.h"
 
 #ifndef IN_GINAC
index 15503454b4260adf3b7c333cbae12b6ed98e937f..41e222fb4d64e4f17ec10cd63e899613db78275b 100644 (file)
@@ -177,6 +177,32 @@ template<> inline bool is_the_function<psi_SERIAL>(const ex & x)
        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)
 
diff --git a/ginac/inifcns_elliptic.cpp b/ginac/inifcns_elliptic.cpp
new file mode 100644 (file)
index 0000000..07b4061
--- /dev/null
@@ -0,0 +1,632 @@
+/** @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
+
diff --git a/ginac/integration_kernel.cpp b/ginac/integration_kernel.cpp
new file mode 100644 (file)
index 0000000..4372f5c
--- /dev/null
@@ -0,0 +1,2146 @@
+/** @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
diff --git a/ginac/integration_kernel.h b/ginac/integration_kernel.h
new file mode 100644 (file)
index 0000000..73e6c45
--- /dev/null
@@ -0,0 +1,668 @@
+/** @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
diff --git a/ginac/utils_multi_iterator.h b/ginac/utils_multi_iterator.h
new file mode 100644 (file)
index 0000000..8e21960
--- /dev/null
@@ -0,0 +1,1488 @@
+/** @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
index e224fda93cafb5a240b298673799678f04ba6d15..0a056c4d5cac11f546ca05e21100289322e9b7c9 100644 (file)
@@ -574,6 +574,105 @@ static ex f_unit(const exprseq &e)
        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)"));
@@ -650,6 +749,24 @@ static const fcn_init builtin_fcns[] = {
        {"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
 };
 
@@ -681,11 +798,15 @@ static const fcn_help_init builtin_help[] = {
        {"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