From 8bcb361ec7767e33ea9b92141a078827e3f3bb0d Mon Sep 17 00:00:00 2001 From: Stefan Weinzierl Date: Sat, 10 Oct 2020 17:11:09 +0200 Subject: [PATCH] Added new routines for the numerical evaluation of iterated integrals like 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 --- check/CMakeLists.txt | 1 + check/Makefile.am | 4 + check/exam_inifcns_elliptic.cpp | 419 ++++++ doc/tutorial/ginac.texi | 81 +- ginac/CMakeLists.txt | 4 + ginac/Makefile.am | 5 +- ginac/ginac.h | 2 + ginac/inifcns.h | 26 + ginac/inifcns_elliptic.cpp | 632 +++++++++ ginac/integration_kernel.cpp | 2146 +++++++++++++++++++++++++++++++ ginac/integration_kernel.h | 668 ++++++++++ ginac/utils_multi_iterator.h | 1488 +++++++++++++++++++++ ginsh/ginsh_parser.ypp | 121 ++ 13 files changed, 5593 insertions(+), 4 deletions(-) create mode 100644 check/exam_inifcns_elliptic.cpp create mode 100644 ginac/inifcns_elliptic.cpp create mode 100644 ginac/integration_kernel.cpp create mode 100644 ginac/integration_kernel.h create mode 100644 ginac/utils_multi_iterator.h diff --git a/check/CMakeLists.txt b/check/CMakeLists.txt index caabbfc3..6b8bd340 100644 --- a/check/CMakeLists.txt +++ b/check/CMakeLists.txt @@ -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 diff --git a/check/Makefile.am b/check/Makefile.am index 1fffd9f7..d1881520 100644 --- a/check/Makefile.am +++ b/check/Makefile.am @@ -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 index 00000000..7539a921 --- /dev/null +++ b/check/exam_inifcns_elliptic.cpp @@ -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 +#include +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(eta).get_numerical_value(qbar_prime,N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_2,N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_2,N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_4,N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_6,3*N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_12,6*N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_2,N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_6,2*N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_12,4*N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_3,N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(eta).get_numerical_value(qbar_prime_6,2*N_trunc); + res2 = pow(c*tau+d,k)*ex_to(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(); +} diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index 1ba42151..520a13df 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -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 diff --git a/ginac/CMakeLists.txt b/ginac/CMakeLists.txt index 08666551..5b0937c4 100644 --- a/ginac/CMakeLists.txt +++ b/ginac/CMakeLists.txt @@ -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 diff --git a/ginac/Makefile.am b/ginac/Makefile.am index cb76517a..0dc82c15 100644 --- a/ginac/Makefile.am +++ b/ginac/Makefile.am @@ -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 \ diff --git a/ginac/ginac.h b/ginac/ginac.h index 84a26f47..f3905dfe 100644 --- a/ginac/ginac.h +++ b/ginac/ginac.h @@ -68,6 +68,8 @@ #include "factor.h" +#include "integration_kernel.h" + #include "excompiler.h" #ifndef IN_GINAC diff --git a/ginac/inifcns.h b/ginac/inifcns.h index 15503454..41e222fb 100644 --- a/ginac/inifcns.h +++ b/ginac/inifcns.h @@ -177,6 +177,32 @@ template<> inline bool is_the_function(const ex & x) return is_the_function(x) || is_the_function(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 +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 +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(const ex& x) +{ + return is_the_function(x) || is_the_function(x); +} + + /** Factorial function. */ DECLARE_FUNCTION_1P(factorial) diff --git a/ginac/inifcns_elliptic.cpp b/ginac/inifcns_elliptic.cpp new file mode 100644 index 00000000..07b40616 --- /dev/null +++ b/ginac/inifcns_elliptic.cpp @@ -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 +#include +#include +#include +#include + +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(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(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(k).to_cl_N())); + + ex result = Pi/2 * numeric( agm_helper_second_kind(1,kbar,ex_to(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(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 & m, const std::vector & 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 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; jseries_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 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; jseries_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 & kernel_in, const cln::cl_N & lambda, int N_trunc) +{ + size_t depth = kernel_in.size(); + + std::vector m; + m.reserve(depth); + + std::vector kernel; + kernel.reserve(depth); + + int n = 1; + + for (const auto & it : kernel_in) { + if ( is_a(*it) ) { + n++; + } + else { + m.push_back(n); + kernel.push_back( &ex_to(*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 & 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(*(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 a,b; + for (size_t i=0; i i_shuffle(a,b); + for( i_shuffle.init(); !i_shuffle.overflow(); i_shuffle++) { + std::vector new_kernel; + new_kernel.reserve(depth); + for (size_t i=0; i(kernel_lst); + + bool flag_not_numeric = false; + for (const auto & it : k_lst) { + if ( !is_a(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(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(N_trunc).to_int(); + + // check trailing zeros + const size_t depth = k_lst.nops(); + + std::vector kernel_vec; + kernel_vec.reserve(depth); + + for (const auto & it : k_lst) { + kernel_vec.push_back( &ex_to(it) ); + } + + size_t i_trailing = 0; + for (size_t i=0; ihas_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(lambda.evalf()).to_cl_N(); + basic_log_kernel L0 = basic_log_kernel(); + + cln::cl_N result; + if ( is_a(*(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(*(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 index 00000000..4372f5c3 --- /dev/null +++ b/ginac/integration_kernel.cpp @@ -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 +#include +#include + + +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(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(e) || is_a(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 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 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(ifactor(abs(n))); + lst p_lst = ex_to(prime_factorisation.op(0)); + lst e_lst = ex_to(prime_factorisation.op(1)); + + size_t n_primes = p_lst.nops(); + + if ( n_primes > 0 ) { + // take the last prime + numeric p = ex_to(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(ifactor(n_odd)); + lst prime_lst = ex_to(temp_lst.op(0)); + lst expo_lst = ex_to(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(*it_p)),ex_to(*it_e)); + } + + return ex_to(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(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(&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(temp.coeff(x,j-1).evalf()).to_cl_N() ); + } + } + else { + for (int j=n_vec; j(lambda.evalf()).to_cl_N(); + cln::cl_N pre_cln = ex_to(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(&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(&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(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(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(&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(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(n).to_int(); + int m_int = ex_to(m).to_int(); + + cln::cl_N x_cln = ex_to(x.evalf()).to_cl_N(); + cln::cl_N y_cln = ex_to(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(&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(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(n).to_int(); + int m_int = ex_to(m).to_int(); + + cln::cl_N x_cln = ex_to(x.evalf()).to_cl_N(); + cln::cl_N y_cln = ex_to(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(&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(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(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(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(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(K).to_int(); + + if ( (i % K_int) != 0 ) { + return 0; + } + int i_local = i/K_int; + + ex w = exp(ex_to((2*Pi*I*z).evalf())); + cln::cl_N w_cln = ex_to(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(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(n); + + if ( n_num == 0 ) { + return 1; + } + + // use the direct formula here + if ( n_num == 1 ) { + ex wbar = exp(ex_to((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(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(&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(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(n); + int n_int = n_num.to_int(); + + ex w_j_inv = exp(ex_to((-2*Pi*I*z_j).evalf())); + cln::cl_N w_j_inv_cln = ex_to(w_j_inv).to_cl_N(); + + ex qbar = exp(ex_to((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((C_norm * 2*Pi*I).evalf()).to_cl_N(); + } + + return 0; + } + + // case n=2 + if ( n_num == 2 ) { + if ( ex_to(z_j.evalf()).is_zero() ) { + if ( i == 0 ) { + return ex_to((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(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(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(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(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(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(&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(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(k); + numeric N_num = ex_to(N); + numeric a_num = ex_to(a); + numeric b_num = ex_to(b); + numeric K_num = ex_to(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(&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(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(k); + numeric r_num = ex_to(r); + numeric s_num = ex_to(s); + numeric N_num = ex_to(N); + + ex res = coefficient_a0(k_num,r_num,s_num,N_num); + + for (numeric i1=1; i1(&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(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(&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(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 index 00000000..73e6c457 --- /dev/null +++ b/ginac/integration_kernel.h @@ -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 +#include + +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 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 index 00000000..8e219601 --- /dev/null +++ b/ginac/utils_multi_iterator.h @@ -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 +#include +#include +#include + +namespace GiNaC { + +/** + * + * SFINAE test for distance + * + */ +template class has_distance { +private: + typedef char yes_type[1]; + typedef char no_type[2]; + + template static yes_type & test( decltype(std::distance) ) ; + template static no_type & test(...); + +public: + enum { value = sizeof(test(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 std::enable_if::value, typename std::iterator_traits::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 std::enable_if::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 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 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 & vv); + + // dtor + virtual ~basic_multi_iterator(); + + // functions +public : + size_t size(void) const; + bool overflow(void) const; + const std::vector & 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 & init(void); + // postfix increment + virtual basic_multi_iterator & operator++ (int); + + // I/O operators + template friend std::ostream & operator<< (std::ostream & os, const basic_multi_iterator & v); + + // member variables : +protected : + T N; + T B; + std::vector 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 multi_iterator_ordered : public basic_multi_iterator { + + // 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 & vv); + + // overriding virtual functions from base class +public : + // initialization + basic_multi_iterator & init(void); + // postfix increment + basic_multi_iterator & operator++ (int); + + // I/O operators + template friend std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered & 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 multi_iterator_ordered_eq : public basic_multi_iterator { + + // 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 & vv); + + // overriding virtual functions from base class +public : + // initialization + basic_multi_iterator & init(void); + // postfix increment + basic_multi_iterator & operator++ (int); + + // I/O operators + template friend std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered_eq & 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 multi_iterator_ordered_eq_indv : public basic_multi_iterator { + + // ctors +public : + multi_iterator_ordered_eq_indv(void); + explicit multi_iterator_ordered_eq_indv(T B, const std::vector & Nv, size_t k); + explicit multi_iterator_ordered_eq_indv(T B, const std::vector & Nv, const std::vector & vv); + + // overriding virtual functions from base class +public : + // initialization + basic_multi_iterator & init(void); + // postfix increment + basic_multi_iterator & operator++ (int); + + // I/O operators + template friend std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered_eq_indv & v); + + // member variables : +protected : + std::vector 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 multi_iterator_counter : public basic_multi_iterator { + + // 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 & vv); + + // overriding virtual functions from base class +public : + // initialization + basic_multi_iterator & init(void); + // postfix increment + basic_multi_iterator & operator++ (int); + + // I/O operators + template friend std::ostream & operator<< (std::ostream & os, const multi_iterator_counter & 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 multi_iterator_counter_indv : public basic_multi_iterator { + + // ctors +public : + multi_iterator_counter_indv(void); + explicit multi_iterator_counter_indv(T B, const std::vector & Nv, size_t k); + explicit multi_iterator_counter_indv(T B, const std::vector & Nv, const std::vector & vv); + + // overriding virtual functions from base class +public : + // initialization + basic_multi_iterator & init(void); + // postfix increment + basic_multi_iterator & operator++ (int); + + // I/O operators + template friend std::ostream & operator<< (std::ostream & os, const multi_iterator_counter_indv & v); + + // member variables : +protected : + std::vector 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 multi_iterator_permutation : public basic_multi_iterator { + + // 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 & vv); + + // overriding virtual functions from base class +public : + // initialization + basic_multi_iterator & init(void); + // postfix increment + basic_multi_iterator & operator++ (int); + + // new functions in this class + int get_sign(void) const; + + // I/O operators + template friend std::ostream & operator<< (std::ostream & os, const multi_iterator_permutation & v); + +}; + +/** + * + * The class multi_iterator_shuffle defines a multi_iterator, + * which runs over all shuffles of a and b. + * + */ +template class multi_iterator_shuffle : public basic_multi_iterator { + + // ctors +public : + multi_iterator_shuffle(void); + explicit multi_iterator_shuffle(const std::vector & a, const std::vector & b); + + // overriding virtual functions from base class +public : + // initialization + basic_multi_iterator & init(void); + // postfix increment + basic_multi_iterator & operator++ (int); + + // I/O operators + template friend std::ostream & operator<< (std::ostream & os, const multi_iterator_shuffle & v); + + // member variables : +protected : + size_t N_internal; + std::vector v_internal; + std::vector 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 multi_iterator_shuffle_prime : public multi_iterator_shuffle { + + // ctors +public : + multi_iterator_shuffle_prime(void); + explicit multi_iterator_shuffle_prime(const std::vector & a, const std::vector & b); + + // overriding virtual functions from base class +public : + // initialization + basic_multi_iterator & init(void); + + // I/O operators + template friend std::ostream & operator<< (std::ostream & os, const multi_iterator_shuffle_prime & v); +}; + +// ---------------------------------------------------------------------------------------------------------------- + +// ctors + +/** + * + * Default constructor + * + */ +template inline basic_multi_iterator::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 inline basic_multi_iterator::basic_multi_iterator(T BB, T NN, size_t k) : N(NN), B(BB), v(k), flag_overflow(false) +{} + +/** + * + * Construct from a vector. + * + */ +template inline basic_multi_iterator::basic_multi_iterator(T BB, T NN, const std::vector & vv) : N(NN), B(BB), v(vv), flag_overflow(false) +{} + +/** + * + * Destructor + * + */ +template inline basic_multi_iterator::~basic_multi_iterator() +{} + +// functions + +/** + * + * Returns the size of a multi_iterator. + * + */ +template inline size_t basic_multi_iterator::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 inline basic_multi_iterator & basic_multi_iterator::init(void) +{ + flag_overflow = false; + + for ( size_t i=0; i inline bool basic_multi_iterator::overflow(void) const +{ + return flag_overflow; +} + +/** + * + * Returns a reference to the vector v. + * + */ +template inline const std::vector & basic_multi_iterator::get_vector(void) const +{ + return v; +} + +// subscripting + +/** + * + * Subscription via [] + * + */ +template inline T basic_multi_iterator::operator[](size_t i) const +{ + return v[i]; +} + +/** + * + * Subscription via [] + * + */ +template inline T & basic_multi_iterator::operator[](size_t i) +{ + return v[i]; +} + +/** + * + * Subscription via () + * + */ +template inline T basic_multi_iterator::operator()(size_t i) const +{ + return v[i]; +} + +/** + * + * Subscription via () + * + */ +template inline T & basic_multi_iterator::operator()(size_t i) +{ + return v[i]; +} + + +/** + * + * No effect for basic_multi_iterator + * + */ +template inline basic_multi_iterator & basic_multi_iterator::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 inline std::ostream & operator<< (std::ostream & os, const basic_multi_iterator & v) +{ + os << "basic_multi_iterator("; + for ( size_t i=0; i0) { + os << ","; + } + os << format_index_value(v.B,v(i)); + } + + return os << ")"; +} + + + +// ctors + +/** + * + * Default constructor + * + */ +template inline multi_iterator_ordered::multi_iterator_ordered(void) : basic_multi_iterator() +{} + +/** + * + * Construct a multi_iterator with upper limit N and size k . + * + */ +template inline multi_iterator_ordered::multi_iterator_ordered(T B, T N, size_t k) : basic_multi_iterator(B,N,k) +{} + +/** + * + * Construct from a vector. + * + */ +template inline multi_iterator_ordered::multi_iterator_ordered(T B, T N, const std::vector & v) : basic_multi_iterator(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 inline basic_multi_iterator & multi_iterator_ordered::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 inline basic_multi_iterator & multi_iterator_ordered::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;jjv[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 inline std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered & v) +{ + os << "multi_iterator_ordered("; + for ( size_t i=0; i0) { + os << ","; + } + os << format_index_value(v.B,v(i)); + } + + return os << ")"; +} + + + +// ctors + +/** + * + * Default constructor + * + */ +template inline multi_iterator_ordered_eq::multi_iterator_ordered_eq(void) : basic_multi_iterator() +{} + +/** + * + * Construct a multi_iterator with upper limit N and size k . + * + */ +template inline multi_iterator_ordered_eq::multi_iterator_ordered_eq(T B, T N, size_t k) : basic_multi_iterator(B,N,k) +{} + +/** + * + * Construct from a vector. + * + */ +template inline multi_iterator_ordered_eq::multi_iterator_ordered_eq(T B, T N, const std::vector & v) : basic_multi_iterator(B,N,v) +{} + +// functions + +/** + * + * Initialize the multi-index to + * \f[ + * (n_1,n_2,...,n_k) = (B,B,...,B) + * \f] + * + */ +template inline basic_multi_iterator & multi_iterator_ordered_eq::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 inline basic_multi_iterator & multi_iterator_ordered_eq::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;jjv[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 inline std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered_eq & v) +{ + os << "multi_iterator_ordered_eq("; + for ( size_t i=0; i0) { + os << ","; + } + os << format_index_value(v.B,v(i)); + } + + return os << ")"; +} + + + + +// ctors + +/** + * + * Default constructor + * + */ +template inline multi_iterator_ordered_eq_indv::multi_iterator_ordered_eq_indv(void) : basic_multi_iterator(), Nv() +{} + +/** + * + * Construct a multi_iterator with upper limit N and size k . + * + */ +template inline multi_iterator_ordered_eq_indv::multi_iterator_ordered_eq_indv(T B, const std::vector & Nvv, size_t k) : basic_multi_iterator(B,B,k), Nv(Nvv) +{} + +/** + * + * Construct from a vector. + * + */ +template inline multi_iterator_ordered_eq_indv::multi_iterator_ordered_eq_indv(T B, const std::vector & Nvv, const std::vector & v) : basic_multi_iterator(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 inline basic_multi_iterator & multi_iterator_ordered_eq_indv::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 inline basic_multi_iterator & multi_iterator_ordered_eq_indv::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;jjv[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 inline std::ostream & operator<< (std::ostream & os, const multi_iterator_ordered_eq_indv & v) +{ + os << "multi_iterator_ordered_eq_indv("; + for ( size_t i=0; i0) { + os << ","; + } + os << format_index_value(v.B,v(i)); + } + + return os << ")"; +} + + + + +// ctors + +/** + * + * Default constructor + * + */ +template inline multi_iterator_counter::multi_iterator_counter(void) : basic_multi_iterator() +{} + +/** + * + * Construct a multi_iterator with upper limit N and size k . + * + */ +template inline multi_iterator_counter::multi_iterator_counter(T B, T N, size_t k) : basic_multi_iterator(B,N,k) +{} + +/** + * + * Construct from a vector. + * + */ +template inline multi_iterator_counter::multi_iterator_counter(T B, T N, const std::vector & v) : basic_multi_iterator(B,N,v) +{} + +// functions + +/** + * + * Initialize the multi-index to + * \f[ + * (n_1,n_2,n_3,...,n_k) = (B,B,...,B) + * \f] + * + */ +template inline basic_multi_iterator & multi_iterator_counter::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 inline basic_multi_iterator & multi_iterator_counter::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 inline std::ostream & operator<< (std::ostream & os, const multi_iterator_counter & v) +{ + os << "multi_iterator_counter("; + for ( size_t i=0; i0) { + os << ","; + } + os << format_index_value(v.B,v(i)); + } + + return os << ")"; +} + + + + +// ctors + +/** + * + * Default constructor + * + */ +template inline multi_iterator_counter_indv::multi_iterator_counter_indv(void) : basic_multi_iterator(), Nv() +{} + +/** + * + * Construct a multi_iterator with upper limit N and size k . + * + */ +template inline multi_iterator_counter_indv::multi_iterator_counter_indv(T B, const std::vector & Nvv, size_t k) : basic_multi_iterator(B,B,k), Nv(Nvv) +{} + +/** + * + * Construct from a vector. + * + */ +template inline multi_iterator_counter_indv::multi_iterator_counter_indv(T B, const std::vector & Nvv, const std::vector & v) : basic_multi_iterator(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 inline basic_multi_iterator & multi_iterator_counter_indv::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 inline basic_multi_iterator & multi_iterator_counter_indv::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 inline std::ostream & operator<< (std::ostream & os, const multi_iterator_counter_indv & v) +{ + os << "multi_iterator_counter_indv("; + for ( size_t i=0; i0) { + os << ","; + } + os << format_index_value(v.B,v(i)); + } + + return os << ")"; +} + + + + +// ctors + +/** + * + * Default constructor + * + */ +template inline multi_iterator_permutation::multi_iterator_permutation(void) : basic_multi_iterator() +{} + +/** + * + * Construct a multi_iterator with upper limit N and size k . + * + */ +template inline multi_iterator_permutation::multi_iterator_permutation(T B, T N, size_t k) : basic_multi_iterator(B,N,k) +{} + +/** + * + * Construct from a vector. + * + */ +template inline multi_iterator_permutation::multi_iterator_permutation(T B, T N, const std::vector & v) : basic_multi_iterator(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 inline basic_multi_iterator & multi_iterator_permutation::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 inline basic_multi_iterator & multi_iterator_permutation::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; iiv[j] == this->v[ii]) { + flag_have_already = true; + } + } + } + + if ( this->v[j] == this->N ) { + j--; + } + else { + break; + } + } + + for (int l=j+1; lv[l] = this->B; + + bool flag_have_already; + do { + flag_have_already = false; + for (int ii=0; iiv[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; iiv[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 inline int multi_iterator_permutation::get_sign() const +{ + int sign = 1; + int k = this->size(); + + for ( int i=0; iv[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 inline std::ostream & operator<< (std::ostream & os, const multi_iterator_permutation & v) +{ + os << "multi_iterator_permutation("; + for ( size_t i=0; i0) { + os << ","; + } + os << format_index_value(v.B,v(i)); + } + + return os << ")"; +} + + + +// ctors + +/** + * + * Default constructor + * + */ +template inline multi_iterator_shuffle::multi_iterator_shuffle(void) : basic_multi_iterator(), N_internal(), v_internal(), v_orig() +{} + +/** + * + * Construct from a vector. + * + */ +template inline multi_iterator_shuffle::multi_iterator_shuffle(const std::vector & a, const std::vector & b) : basic_multi_iterator(), N_internal(), v_internal(), v_orig() +{ + this->B = a[0]; + + for (size_t i=0; iv.push_back( a[i] ); + this->v_orig.push_back( a[i] ); + this->v_internal.push_back( i ); + } + for (size_t i=0; iv.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 inline basic_multi_iterator & multi_iterator_shuffle::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 inline basic_multi_iterator & multi_iterator_shuffle::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;jjv_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; jv_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 inline std::ostream & operator<< (std::ostream & os, const multi_iterator_shuffle & v) +{ + os << "multi_iterator_shuffle("; + for ( size_t i=0; i0) { + os << ","; + } + os << format_index_value(v.B,v(i)); + } + + return os << ")"; +} + + + +// ctors + +/** + * + * Default constructor + * + */ +template inline multi_iterator_shuffle_prime::multi_iterator_shuffle_prime(void) : multi_iterator_shuffle() +{} + +/** + * + * Construct from a vector. + * + */ +template inline multi_iterator_shuffle_prime::multi_iterator_shuffle_prime(const std::vector & a, const std::vector & b) : multi_iterator_shuffle(a,b) +{} + +// functions + +/** + * + * Initialize the multi-index to the first shuffle. + * + */ +template inline basic_multi_iterator & multi_iterator_shuffle_prime::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 inline std::ostream & operator<< (std::ostream & os, const multi_iterator_shuffle_prime & v) +{ + os << "multi_iterator_shuffle_prime("; + for ( size_t i=0; i0) { + os << ","; + } + os << format_index_value(v.B,v(i)); + } + + return os << ")"; +} + +} // namespace GiNaC + +#endif // ndef GINAC_UTILS_MULTI_ITERATOR_H diff --git a/ginsh/ginsh_parser.ypp b/ginsh/ginsh_parser.ypp index e224fda9..0a056c4d 100644 --- a/ginsh/ginsh_parser.ypp +++ b/ginsh/ginsh_parser.ypp @@ -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(e[0]) ) { + return ex_to(e[0]).q_expansion_modular_form(e[1], ex_to(e[2]).to_int()); + } + if ( is_a(e[0]) ) { + return ex_to(e[0]).q_expansion_modular_form(e[1], ex_to(e[2]).to_int()); + } + if ( is_a(e[0]) ) { + return ex_to(e[0]).q_expansion_modular_form(e[1], ex_to(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 -- 2.45.0