1 /** @file integration_kernel.cpp
3 * Implementation of GiNaC's integration kernels for iterated integrals. */
6 * GiNaC Copyright (C) 1999-2021 Johannes Gutenberg University Mainz, Germany
8 * This program is free software; you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation; either version 2 of the License, or
11 * (at your option) any later version.
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
18 * You should have received a copy of the GNU General Public License
19 * along with this program; if not, write to the Free Software
20 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 #include "integration_kernel.h"
26 #include "operators.h"
28 #include "relational.h"
44 // anonymous namespace for helper function
49 * Returns the Kronecker symbol in the case where
51 * n: unit or prime number
53 * If n is an odd prime number, the routine returns the Legendre symbol.
55 * 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)
56 * the routine returns the special values defined below.
58 * Ref.: Toshitsune Miyake, Modular Forms, Chapter 3.1
61 numeric kronecker_symbol_prime(const numeric & a, const numeric & n)
75 if ( GiNaC::smod(a,8) == 1 ) {
78 else if ( GiNaC::smod(a,8) == -1 ) {
81 else if ( GiNaC::smod(a,8) == 3 ) {
84 else if ( GiNaC::smod(a,8) == -3 ) {
92 // n is an odd prime number
93 return GiNaC::smod( pow(a,(n-1)/2), n);
100 * a: discriminant of a quadratic field, defines primitive character phi
101 * b: discriminant of a quadratic field, defines primitive character psi
102 * L=|a|: conductor of primitive character phi
103 * M=|b|: conductor of primitive character psi
107 * This function computes
109 * \sum\limits_{d | n} \psi(d) \phi(n/d) d^{k-1}
112 * Ref.: Eq.(5.3.1), William A. Stein, Modular Forms, A computational Approach;
114 * Eq.(32), arXiv:1704.08895
117 numeric divisor_function(const numeric & n, const numeric & a, const numeric & b, const numeric & k)
121 for (numeric i1=1; i1<=n; i1++) {
122 if ( irem(n,i1) == 0 ) {
123 numeric ratio = n/i1;
124 res += primitive_dirichlet_character(ratio,a) * primitive_dirichlet_character(i1,b) * pow(i1,k-1);
128 return ex_to<numeric>(res);
135 * a: discriminant of a quadratic field, defines primitive character phi
136 * b: discriminant of a quadratic field, defines primitive character psi
137 * L=|a|: conductor of primitive character phi
138 * M=|b|: conductor of primitive character psi
140 * This function computes the constant term of the q-expansion of an Eisenstein series.
142 * The coefficient is given by the equation below eq.(5.3.1) in William A. Stein, Modular Forms, A computational Approach;
144 * or by eq.(33), arXiv:1704.08895
147 numeric coefficient_a0(const numeric & k, const numeric & a, const numeric & b)
149 ex conductor = abs(a);
152 if ( conductor == 1 ) {
153 a0 = -numeric(1,2)/k*generalised_Bernoulli_number(k,b);
165 * q: exp(2 Pi I tau/M)
167 * a: discriminant of a quadratic field, defines primitive character phi
168 * b: discriminant of a quadratic field, defines primitive character psi
169 * L=|a|: conductor of primitive character phi
170 * M=|b|: conductor of primitive character psi
172 * N: truncation order
174 * Returns the q-expansion of an Eisenstein to order N (the q^(N-1)-term is included, q^N is neglected).
176 * Ref.: Eq.(5.3.1), William A. Stein, Modular Forms, A computational Approach;
178 * Eq.(32), arXiv:1704.08895
181 ex eisenstein_series(const numeric & k, const ex & q, const numeric & a, const numeric & b, const numeric & N)
183 ex res = coefficient_a0(k,a,b);
185 for (numeric i1=1; i1<N; i1++) {
186 res += divisor_function(i1,a,b,k) * pow(q,i1);
194 * Returns the q_N-expansion of the Eisenstein series E_{k,a,b}(K tau_N)
197 ex E_eisenstein_series(const ex & q, const numeric & k, const numeric & N_level, const numeric & a, const numeric & b, const numeric & K, const numeric & N_order)
199 int N_order_int = N_order.to_int();
201 ex res = eisenstein_series(k,pow(q,K),a,b,iquo(N_order,K));
203 res += Order(pow(q,N_order_int));
204 res = res.series(q,N_order_int);
211 * In weight 2 we have a special case for trivial Dirichlet characters:
212 * Returns the q_N-expansion of the Eisenstein series E_{2,1,1}(tau_N) - K E_{2,1,1}(K tau_N).
215 ex B_eisenstein_series(const ex & q, const numeric & N_level, const numeric & K, const numeric & N_order)
217 int N_order_int = N_order.to_int();
219 ex res = eisenstein_series(2,q,1,1,N_order) - K*eisenstein_series(2,pow(q,K),1,1,iquo(N_order,K));
221 res += Order(pow(q,N_order_int));
222 res = res.series(q,N_order_int);
229 * A helper function to expand polynomials in Eisenstein series.
232 struct subs_q_expansion : public map_function
234 subs_q_expansion(const ex & arg_qbar, int arg_order) : qbar(arg_qbar), order(arg_order)
237 ex operator()(const ex & e)
239 if ( is_a<Eisenstein_kernel>(e) || is_a<Eisenstein_h_kernel>(e) ) {
240 return series_to_poly(e.series(qbar,order));
257 * This is a rational function in x.
259 * To speed things up, we cache it.
268 // non-virtual functions
270 ex get_symbolic_value(int n, const ex & x_val);
271 ex get_numerical_value(int n, const ex & x_val);
273 // member variables :
275 static std::vector<ex> cache_vec;
280 Li_negative::Li_negative() {}
282 ex Li_negative::get_symbolic_value(int n, const ex & x_val)
284 int n_cache = cache_vec.size();
286 if ( n >= n_cache ) {
287 for (int j=n_cache; j<=n; j++) {
293 f = normal( x*diff(cache_vec[j-1],x));
295 cache_vec.push_back( f );
299 return cache_vec[n].subs(x==x_val);
302 ex Li_negative::get_numerical_value(int n, const ex & x_val)
304 symbol x_symb("x_symb");
306 ex f = this->get_symbolic_value(n,x_symb);
308 ex res = f.subs(x_symb==x_val).evalf();
313 // initialise static data members
314 std::vector<ex> Li_negative::cache_vec;
315 symbol Li_negative::x = symbol("x");
318 } // end of anonymous namespace
322 * Returns the decomposition of the positive integer n into prime numbers
324 * lst( lst(p1,...,pr), lst(a1,...,ar) )
326 * n = p1^a1 ... pr^ar.
329 ex ifactor(const numeric & n)
331 if ( !n.is_pos_integer() ) throw (std::runtime_error("ifactor(): argument not a positive integer"));
335 // implementation for small integers
337 for (numeric p=2; p<=n; p++) {
338 if ( p.info(info_flags::prime) ) {
339 numeric exp_temp = 0;
340 while ( irem(n_temp, p) == 0 ) {
346 exp_lst.append(exp_temp);
349 if ( n_temp == 1 ) break;
352 if ( n_temp != 1 ) throw (std::runtime_error("ifactor(): probabilistic primality test failed"));
354 lst res = {p_lst,exp_lst};
361 * Returns true if the integer n is either one or the discriminant of a quadratic number field.
363 * Returns false otherwise.
365 * Ref.: Toshitsune Miyake, Modular Forms, Chapter 3.1
368 bool is_discriminant_of_quadratic_number_field(const numeric & n)
378 lst prime_factorisation = ex_to<lst>(ifactor(abs(n)));
379 lst p_lst = ex_to<lst>(prime_factorisation.op(0));
380 lst e_lst = ex_to<lst>(prime_factorisation.op(1));
382 size_t n_primes = p_lst.nops();
384 if ( n_primes > 0 ) {
385 // take the last prime
386 numeric p = ex_to<numeric>(p_lst.op(n_primes-1));
389 if ( e_lst.op(n_primes-1) != 1 ) {
394 if ( mod(p,4) == 3 ) {
397 return is_discriminant_of_quadratic_number_field(n/pstar);
401 if ( (n==-4) || (n==-8) || (n==8) || (n==-32) || (n==32) || (n==-64) || (n==128) ) {
410 * Returns the Kronecker symbol
414 * This routine defines
415 * kronecker_symbol(1,0) = 1
416 * kronecker_symbol(-1,0) = 1
417 * kronecker_symbol(a,0) = 0, a != 1,-1
420 * kronecker_symbol(-1,0) = 1 (in agreement with Sage)
422 * Ref.: Toshitsune Miyake, Modular Forms, Chapter 3.1
425 numeric kronecker_symbol(const numeric & a, const numeric & n)
427 // case n=0 first, include kronecker_symbol(0,0)=0
429 if ( (a == 1) || (a == -1) ) {
444 ex res = kronecker_symbol_prime(a,unit);
446 numeric n_odd = n_pos;
448 while ( n_odd.is_even() ) {
453 res *= pow(kronecker_symbol_prime(a,2),alpha);
456 lst temp_lst = ex_to<lst>(ifactor(n_odd));
457 lst prime_lst = ex_to<lst>(temp_lst.op(0));
458 lst expo_lst = ex_to<lst>(temp_lst.op(1));
460 for (auto it_p = prime_lst.begin(), it_e = expo_lst.begin(); it_p != prime_lst.end(); it_p++, it_e++) {
461 res *= pow(kronecker_symbol_prime(a,ex_to<numeric>(*it_p)),ex_to<numeric>(*it_e));
464 return ex_to<numeric>(res);
469 * Defines a primitive Dirichlet character through the Kronecker symbol.
472 * a: discriminant of a quadratic field
475 * The character takes the values -1,0,1.
478 numeric primitive_dirichlet_character(const numeric & n, const numeric & a)
480 return kronecker_symbol(a,n);
485 * Defines a Dirichlet character through the Kronecker symbol.
488 * a: discriminant of a quadratic field
490 * N: modulus, needs to be multiple of |a|
492 * The character takes the values -1,0,1.
495 numeric dirichlet_character(const numeric & n, const numeric & a, const numeric & N)
497 if ( gcd(n,N) == 1 ) {
498 return primitive_dirichlet_character(n,a);
506 * The generalised Bernoulli number.
508 * k: index / modular weight
510 * b: discriminant of a quadratic field, defines primitive character psi
511 * M=|b|: conductor of primitive character psi
513 * The generalised Bernoulli number is computed from the series expansion of the generating function.
514 * The generating function is given in eq.(34), arXiv:1704.08895
517 numeric generalised_Bernoulli_number(const numeric & k, const numeric & b)
519 int k_int = k.to_int();
523 numeric conductor = abs(b);
526 for (numeric i1=1; i1<=conductor; i1++) {
527 gen_fct += primitive_dirichlet_character(i1,b) * x*exp(i1*x)/(exp(conductor*x)-1);
530 gen_fct = series_to_poly(gen_fct.series(x,k_int+1));
532 ex B = factorial(k) * gen_fct.coeff(x,k_int);
534 return ex_to<numeric>(B);
539 * The Bernoulli polynomials
542 ex Bernoulli_polynomial(const numeric & k, const ex & x)
544 int k_int = k.to_int();
548 ex gen_fct = t*exp(x*t)/(exp(t)-1);
550 gen_fct = series_to_poly(gen_fct.series(t,k_int+1));
552 ex B = factorial(k) * gen_fct.coeff(t,k_int);
559 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integration_kernel, basic,
560 print_func<print_context>(&integration_kernel::do_print))
562 integration_kernel::integration_kernel() : inherited(), cache_step_size(100), series_vec()
566 int integration_kernel::compare_same_type(const basic &other) const
571 ex integration_kernel::series(const relational & r, int order, unsigned options) const
573 if ( r.rhs() != 0 ) {
574 throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
577 return Laurent_series(r.lhs(),order);
582 * This routine returns true, if the integration kernel has a trailing zero.
585 bool integration_kernel::has_trailing_zero(void) const
587 if ( cln::zerop( series_coeff(0) ) ) {
596 * This routine returns true, if the integration kernel can be evaluated numerically.
599 bool integration_kernel::is_numeric(void) const
606 * Subclasses have either to implement series_coeff_impl
607 * or the two methods Laurent_series and uses_Laurent_series.
609 * The method series_coeff_impl can be used, if a single coefficient can be computed
610 * independently of the others.
612 * The method Laurent_series can be used, if it is more efficient to compute a Laurent series
613 * in one shot and to determine a range of coefficients from this Laurent series.
616 cln::cl_N integration_kernel::series_coeff(int i) const
618 int n_vec = series_vec.size();
621 int N = cache_step_size*(i/cache_step_size+1);
623 if ( uses_Laurent_series() ) {
625 // series_vec[0] gives coefficient of 1/z, series_vec[N-1] coefficient of z^(N-2),
626 // thus expansion up to order (N-1) is required
627 ex temp = Laurent_series(x, N-1);
628 for (int j=n_vec; j<N; j++) {
629 series_vec.push_back( ex_to<numeric>(temp.coeff(x,j-1).evalf()).to_cl_N() );
633 for (int j=n_vec; j<N; j++) {
634 series_vec.push_back( series_coeff_impl(j) );
639 return series_vec[i];
644 * For \f$ \omega = d\lambda \f$ only the coefficient of \f$ \lambda^0 \f$ is non-zero.
646 * The i-th coefficient corresponds to the power \f$ \lambda^{i-1} \f$.
649 cln::cl_N integration_kernel::series_coeff_impl(int i) const
660 * Returns the Laurent series, starting possibly with the pole term.
661 * Neglected terms are of order \f$ O(x^order) \f$.
664 ex integration_kernel::Laurent_series(const ex & x, int order) const
667 for (int n=-1; n<order; n++) {
668 res += numeric(series_coeff(n+1)) * pow(x,n);
670 res += Order(pow(x,order));
671 res = res.series(x,order);
678 * Evaluates the integrand at lambda.
681 ex integration_kernel::get_numerical_value(const ex & lambda, int N_trunc) const
683 return get_numerical_value_impl(lambda, 1, 0, N_trunc);
688 * Returns true, if the coefficients are computed from the Laurent series
689 * (in which case the method Laurent_series needs to be implemented).
691 * Returns false if this is not the case
692 * (and the class has an implementation of series_coeff_impl).
695 bool integration_kernel::uses_Laurent_series() const
702 * Returns the current size of the cache.
705 size_t integration_kernel::get_cache_size(void) const
707 return series_vec.size();
712 * Sets the step size by which the cache is increased.
715 void integration_kernel::set_cache_step(int cache_steps) const
717 cache_step_size = cache_steps;
722 * Wrapper around series_coeff(i), converts cl_N to numeric.
725 ex integration_kernel::get_series_coeff(int i) const
727 return numeric(series_coeff(i));
732 * The actual implementation for computing a numerical value for the integrand.
735 ex integration_kernel::get_numerical_value_impl(const ex & lambda, const ex & pre, int shift, int N_trunc) const
737 cln::cl_N lambda_cln = ex_to<numeric>(lambda.evalf()).to_cl_N();
738 cln::cl_N pre_cln = ex_to<numeric>(pre.evalf()).to_cl_N();
740 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
746 if ( N_trunc == 0 ) {
747 // sum until precision is reached
748 bool flag_accidental_zero = false;
755 subexpr = series_coeff(N);
757 res += pre_cln * subexpr * cln::expt(lambda_cln,N-1+shift);
759 flag_accidental_zero = cln::zerop(subexpr);
762 } while ( (res != resbuf) || flag_accidental_zero );
765 // N_trunc > 0, sum up the first N_trunc terms
766 for (int N=0; N<N_trunc; N++) {
767 subexpr = series_coeff(N);
769 res += pre_cln * subexpr * cln::expt(lambda_cln,N-1+shift);
776 void integration_kernel::do_print(const print_context & c, unsigned level) const
778 c.s << "integration_kernel()";
781 GINAC_BIND_UNARCHIVER(integration_kernel);
784 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic_log_kernel, integration_kernel,
785 print_func<print_context>(&basic_log_kernel::do_print))
787 basic_log_kernel::basic_log_kernel() : inherited()
791 int basic_log_kernel::compare_same_type(const basic &other) const
796 cln::cl_N basic_log_kernel::series_coeff_impl(int i) const
805 void basic_log_kernel::do_print(const print_context & c, unsigned level) const
807 c.s << "basic_log_kernel()";
810 GINAC_BIND_UNARCHIVER(basic_log_kernel);
813 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(multiple_polylog_kernel, integration_kernel,
814 print_func<print_context>(&multiple_polylog_kernel::do_print))
816 multiple_polylog_kernel::multiple_polylog_kernel() : inherited(), z(_ex1)
820 multiple_polylog_kernel::multiple_polylog_kernel(const ex & arg_z) : inherited(), z(arg_z)
824 int multiple_polylog_kernel::compare_same_type(const basic &other) const
826 const multiple_polylog_kernel &o = static_cast<const multiple_polylog_kernel &>(other);
828 return z.compare(o.z);
831 size_t multiple_polylog_kernel::nops() const
836 ex multiple_polylog_kernel::op(size_t i) const
839 throw(std::range_error("multiple_polylog_kernel::op(): out of range"));
845 ex & multiple_polylog_kernel::let_op(size_t i)
847 ensure_if_modifiable();
850 throw(std::range_error("multiple_polylog_kernel::let_op(): out of range"));
856 bool multiple_polylog_kernel::is_numeric(void) const
858 return z.evalf().info(info_flags::numeric);
861 cln::cl_N multiple_polylog_kernel::series_coeff_impl(int i) const
867 return -cln::expt(ex_to<numeric>(z.evalf()).to_cl_N(),-i);
870 void multiple_polylog_kernel::do_print(const print_context & c, unsigned level) const
872 c.s << "multiple_polylog_kernel(";
877 GINAC_BIND_UNARCHIVER(multiple_polylog_kernel);
880 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(ELi_kernel, integration_kernel,
881 print_func<print_context>(&ELi_kernel::do_print))
883 ELi_kernel::ELi_kernel() : inherited(), n(_ex0), m(_ex0), x(_ex0), y(_ex0)
887 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)
891 int ELi_kernel::compare_same_type(const basic &other) const
893 const ELi_kernel &o = static_cast<const ELi_kernel &>(other);
896 cmpval = n.compare(o.n);
901 cmpval = m.compare(o.m);
906 cmpval = x.compare(o.x);
911 return y.compare(o.y);
914 size_t ELi_kernel::nops() const
919 ex ELi_kernel::op(size_t i) const
931 throw (std::out_of_range("ELi_kernel::op() out of range"));
935 ex & ELi_kernel::let_op(size_t i)
937 ensure_if_modifiable();
949 throw (std::out_of_range("ELi_kernel::let_op() out of range"));
953 bool ELi_kernel::is_numeric(void) const
955 return (n.info(info_flags::nonnegint) && m.info(info_flags::numeric) && x.evalf().info(info_flags::numeric) && y.evalf().info(info_flags::numeric));
958 cln::cl_N ELi_kernel::series_coeff_impl(int i) const
964 int n_int = ex_to<numeric>(n).to_int();
965 int m_int = ex_to<numeric>(m).to_int();
967 cln::cl_N x_cln = ex_to<numeric>(x.evalf()).to_cl_N();
968 cln::cl_N y_cln = ex_to<numeric>(y.evalf()).to_cl_N();
970 cln::cl_N res_cln = 0;
972 for (int j=1; j<=i; j++) {
973 if ( (i % j) == 0 ) {
976 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);
985 * Returns the value of ELi_{n,m}(x,y,qbar)
988 ex ELi_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
990 return get_numerical_value_impl(qbar, 1, 1, N_trunc);
993 void ELi_kernel::do_print(const print_context & c, unsigned level) const
995 c.s << "ELi_kernel(";
1006 GINAC_BIND_UNARCHIVER(ELi_kernel);
1009 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Ebar_kernel, integration_kernel,
1010 print_func<print_context>(&Ebar_kernel::do_print))
1012 Ebar_kernel::Ebar_kernel() : inherited(), n(_ex0), m(_ex0), x(_ex0), y(_ex0)
1016 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)
1020 int Ebar_kernel::compare_same_type(const basic &other) const
1022 const Ebar_kernel &o = static_cast<const Ebar_kernel &>(other);
1025 cmpval = n.compare(o.n);
1030 cmpval = m.compare(o.m);
1035 cmpval = x.compare(o.x);
1040 return y.compare(o.y);
1043 size_t Ebar_kernel::nops() const
1048 ex Ebar_kernel::op(size_t i) const
1060 throw (std::out_of_range("Ebar_kernel::op() out of range"));
1064 ex & Ebar_kernel::let_op(size_t i)
1066 ensure_if_modifiable();
1078 throw (std::out_of_range("Ebar_kernel::let_op() out of range"));
1082 bool Ebar_kernel::is_numeric(void) const
1084 return (n.info(info_flags::nonnegint) && m.info(info_flags::numeric) && x.evalf().info(info_flags::numeric) && y.evalf().info(info_flags::numeric));
1087 cln::cl_N Ebar_kernel::series_coeff_impl(int i) const
1093 int n_int = ex_to<numeric>(n).to_int();
1094 int m_int = ex_to<numeric>(m).to_int();
1096 cln::cl_N x_cln = ex_to<numeric>(x.evalf()).to_cl_N();
1097 cln::cl_N y_cln = ex_to<numeric>(y.evalf()).to_cl_N();
1099 cln::cl_N res_cln = 0;
1101 for (int j=1; j<=i; j++) {
1102 if ( (i % j) == 0 ) {
1105 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);
1114 * Returns the value of Ebar_{n,m}(x,y,qbar)
1117 ex Ebar_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
1119 return get_numerical_value_impl(qbar, 1, 1, N_trunc);
1122 void Ebar_kernel::do_print(const print_context & c, unsigned level) const
1124 c.s << "Ebar_kernel(";
1135 GINAC_BIND_UNARCHIVER(Ebar_kernel);
1138 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Kronecker_dtau_kernel, integration_kernel,
1139 print_func<print_context>(&Kronecker_dtau_kernel::do_print))
1141 Kronecker_dtau_kernel::Kronecker_dtau_kernel() : inherited(), n(_ex0), z(_ex0), K(_ex1), C_norm(_ex1)
1145 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)
1149 int Kronecker_dtau_kernel::compare_same_type(const basic &other) const
1151 const Kronecker_dtau_kernel &o = static_cast<const Kronecker_dtau_kernel &>(other);
1154 cmpval = n.compare(o.n);
1159 cmpval = z.compare(o.z);
1164 cmpval = K.compare(o.K);
1169 return C_norm.compare(o.C_norm);
1172 size_t Kronecker_dtau_kernel::nops() const
1177 ex Kronecker_dtau_kernel::op(size_t i) const
1189 throw (std::out_of_range("Kronecker_dtau_kernel::op() out of range"));
1193 ex & Kronecker_dtau_kernel::let_op(size_t i)
1195 ensure_if_modifiable();
1207 throw (std::out_of_range("Kronecker_dtau_kernel::let_op() out of range"));
1211 bool Kronecker_dtau_kernel::is_numeric(void) const
1213 return (n.info(info_flags::nonnegint) && z.evalf().info(info_flags::numeric) && K.info(info_flags::posint) && C_norm.evalf().info(info_flags::numeric));
1216 cln::cl_N Kronecker_dtau_kernel::series_coeff_impl(int i) const
1218 numeric n_num = ex_to<numeric>(n);
1219 int n_int = n_num.to_int();
1226 return ex_to<numeric>(res.evalf()).to_cl_N();
1239 ex res = C_norm*K / factorial(n_num-2) * bernoulli(n_num)/n_num;
1241 return ex_to<numeric>(res.evalf()).to_cl_N();
1246 // if K>1 the variable i needs to be a multiple of K
1247 int K_int = ex_to<numeric>(K).to_int();
1249 if ( (i % K_int) != 0 ) {
1252 int i_local = i/K_int;
1254 ex w = exp(ex_to<numeric>((2*Pi*I*z).evalf()));
1255 cln::cl_N w_cln = ex_to<numeric>(w).to_cl_N();
1256 cln::cl_N res_cln = 0;
1257 for (int j=1; j<=i_local; j++) {
1258 if ( (i_local % j) == 0 ) {
1259 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);
1262 ex pre = -C_norm*K/factorial(n_num-2);
1264 return ex_to<numeric>(pre.evalf()).to_cl_N() * res_cln;
1269 * Returns the value of the g^(n)(z,K*tau), where tau is given by qbar.
1272 ex Kronecker_dtau_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
1274 numeric n_num = ex_to<numeric>(n);
1280 // use the direct formula here
1282 ex wbar = exp(ex_to<numeric>((2*Pi*I*z).evalf()));
1283 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));
1285 return ex_to<numeric>(res.evalf());
1288 ex pre = pow(2*Pi*I,n_num)/C_norm/K/(n_num-1);
1290 return get_numerical_value_impl(qbar, pre, 1, N_trunc);
1293 void Kronecker_dtau_kernel::do_print(const print_context & c, unsigned level) const
1295 c.s << "Kronecker_dtau_kernel(";
1306 GINAC_BIND_UNARCHIVER(Kronecker_dtau_kernel);
1309 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Kronecker_dz_kernel, integration_kernel,
1310 print_func<print_context>(&Kronecker_dz_kernel::do_print))
1312 Kronecker_dz_kernel::Kronecker_dz_kernel() : inherited(), n(_ex0), z_j(_ex0), tau(_ex0), K(_ex1), C_norm(_ex1)
1316 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)
1320 int Kronecker_dz_kernel::compare_same_type(const basic &other) const
1322 const Kronecker_dz_kernel &o = static_cast<const Kronecker_dz_kernel &>(other);
1325 cmpval = n.compare(o.n);
1330 cmpval = z_j.compare(o.z_j);
1335 cmpval = tau.compare(o.tau);
1340 cmpval = K.compare(o.K);
1345 return C_norm.compare(o.C_norm);
1348 size_t Kronecker_dz_kernel::nops() const
1353 ex Kronecker_dz_kernel::op(size_t i) const
1367 throw (std::out_of_range("Kronecker_dz_kernel::op() out of range"));
1371 ex & Kronecker_dz_kernel::let_op(size_t i)
1373 ensure_if_modifiable();
1387 throw (std::out_of_range("Kronecker_dz_kernel::let_op() out of range"));
1391 bool Kronecker_dz_kernel::is_numeric(void) const
1393 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));
1396 cln::cl_N Kronecker_dz_kernel::series_coeff_impl(int i) const
1398 numeric n_num = ex_to<numeric>(n);
1399 int n_int = n_num.to_int();
1401 ex w_j_inv = exp(ex_to<numeric>((-2*Pi*I*z_j).evalf()));
1402 cln::cl_N w_j_inv_cln = ex_to<numeric>(w_j_inv).to_cl_N();
1404 ex qbar = exp(ex_to<numeric>((2*Pi*I*K*tau).evalf()));
1414 return ex_to<numeric>((C_norm * 2*Pi*I).evalf()).to_cl_N();
1422 if ( ex_to<numeric>(z_j.evalf()).is_zero() ) {
1424 return ex_to<numeric>((C_norm).evalf()).to_cl_N();
1426 else if ( i == 1 ) {
1430 ex res = -bernoulli(i)/numeric(i);
1431 if ( numeric(i).is_even() ) {
1432 Ebar_kernel Ebar = Ebar_kernel( 1-i, 0, numeric(1), numeric(1) );
1433 res += Ebar.get_numerical_value(qbar);
1436 res *= -pow(2*Pi*I,i)*C_norm/factorial(i-1);
1438 return ex_to<numeric>(res.evalf()).to_cl_N();
1447 Li_negative my_Li_negative;
1454 Ebar_kernel Ebar = Ebar_kernel( 1-i, 0, w_j_inv, numeric(1) );
1456 res += my_Li_negative.get_numerical_value(i-1,w_j_inv) + Ebar.get_numerical_value(qbar);
1458 res *= -pow(2*Pi*I,i)*C_norm/factorial(i-1);
1460 return ex_to<numeric>(res.evalf()).to_cl_N();
1468 res += - bernoulli(n_num-1)/(n_num-1);
1471 if ( ex_to<numeric>(z_j.evalf()).is_zero() ) {
1472 if ( (numeric(i)+n_num).is_even() ) {
1473 Ebar_kernel Ebar = Ebar_kernel( 1-i, 2-n_num, numeric(1), numeric(1) );
1475 res += pow(2*Pi*I,i-1)/factorial(i-1) * Ebar.get_numerical_value(qbar);
1480 Ebar_kernel Ebar = Ebar_kernel( 1-i, 2-n_num, w_j_inv, numeric(1) );
1482 res += pow(2*Pi*I,i-1)/factorial(i-1) * Ebar.get_numerical_value(qbar);
1486 res *= - C_norm * 2*Pi*I/factorial(n_num-2);
1488 return ex_to<numeric>(res.evalf()).to_cl_N();
1493 * Returns the value of the g^(n-1)(z-z_j,K*tau).
1496 ex Kronecker_dz_kernel::get_numerical_value(const ex & z, int N_trunc) const
1498 numeric n_num = ex_to<numeric>(n);
1504 ex pre = pow(2*Pi*I,n-2)/C_norm;
1506 return get_numerical_value_impl(z, pre, 0, N_trunc);
1509 void Kronecker_dz_kernel::do_print(const print_context & c, unsigned level) const
1511 c.s << "Kronecker_dz_kernel(";
1524 GINAC_BIND_UNARCHIVER(Kronecker_dz_kernel);
1527 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eisenstein_kernel, integration_kernel,
1528 print_func<print_context>(&Eisenstein_kernel::do_print))
1530 Eisenstein_kernel::Eisenstein_kernel() : inherited(), k(_ex0), N(_ex0), a(_ex0), b(_ex0), K(_ex0), C_norm(_ex1)
1534 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)
1538 int Eisenstein_kernel::compare_same_type(const basic &other) const
1540 const Eisenstein_kernel &o = static_cast<const Eisenstein_kernel &>(other);
1543 cmpval = k.compare(o.k);
1548 cmpval = N.compare(o.N);
1553 cmpval = a.compare(o.a);
1558 cmpval = b.compare(o.b);
1563 cmpval = K.compare(o.K);
1568 return C_norm.compare(o.C_norm);
1573 * The series method for this class returns the qbar-expansion of the modular form,
1574 * without an additional factor of C_norm/qbar.
1576 * This allows for easy use in the class modular_form_kernel.
1579 ex Eisenstein_kernel::series(const relational & r, int order, unsigned options) const
1581 if ( r.rhs() != 0 ) {
1582 throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
1586 ex res = q_expansion_modular_form(qbar, order);
1587 res = res.series(qbar,order);
1592 size_t Eisenstein_kernel::nops() const
1597 ex Eisenstein_kernel::op(size_t i) const
1613 throw (std::out_of_range("Eisenstein_kernel::op() out of range"));
1617 ex & Eisenstein_kernel::let_op(size_t i)
1619 ensure_if_modifiable();
1635 throw (std::out_of_range("Eisenstein_kernel::let_op() out of range"));
1639 bool Eisenstein_kernel::is_numeric(void) const
1641 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));
1644 ex Eisenstein_kernel::Laurent_series(const ex & x, int order) const
1646 ex res = C_norm * q_expansion_modular_form(x, order+1)/x;
1647 res = res.series(x,order);
1654 * Returns the value of the modular form.
1657 ex Eisenstein_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
1659 ex pre = numeric(1)/C_norm;
1661 return get_numerical_value_impl(qbar, pre, 1, N_trunc);
1664 bool Eisenstein_kernel::uses_Laurent_series() const
1669 ex Eisenstein_kernel::q_expansion_modular_form(const ex & q, int order) const
1671 numeric k_num = ex_to<numeric>(k);
1672 numeric N_num = ex_to<numeric>(N);
1673 numeric a_num = ex_to<numeric>(a);
1674 numeric b_num = ex_to<numeric>(b);
1675 numeric K_num = ex_to<numeric>(K);
1677 if ( (k==2) && (a==1) && (b==1) ) {
1678 return B_eisenstein_series(q, N_num, K_num, order);
1681 return E_eisenstein_series(q, k_num, N_num, a_num, b_num, K_num, order);
1684 void Eisenstein_kernel::do_print(const print_context & c, unsigned level) const
1686 c.s << "Eisenstein_kernel(";
1701 GINAC_BIND_UNARCHIVER(Eisenstein_kernel);
1704 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eisenstein_h_kernel, integration_kernel,
1705 print_func<print_context>(&Eisenstein_h_kernel::do_print))
1707 Eisenstein_h_kernel::Eisenstein_h_kernel() : inherited(), k(_ex0), N(_ex0), r(_ex0), s(_ex0), C_norm(_ex1)
1711 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)
1715 int Eisenstein_h_kernel::compare_same_type(const basic &other) const
1717 const Eisenstein_h_kernel &o = static_cast<const Eisenstein_h_kernel &>(other);
1720 cmpval = k.compare(o.k);
1725 cmpval = N.compare(o.N);
1730 cmpval = r.compare(o.r);
1735 cmpval = s.compare(o.s);
1740 return C_norm.compare(o.C_norm);
1745 * The series method for this class returns the qbar-expansion of the modular form,
1746 * without an additional factor of C_norm/qbar.
1748 * This allows for easy use in the class modular_form_kernel.
1751 ex Eisenstein_h_kernel::series(const relational & r, int order, unsigned options) const
1753 if ( r.rhs() != 0 ) {
1754 throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
1758 ex res = q_expansion_modular_form(qbar, order);
1759 res = res.series(qbar,order);
1764 size_t Eisenstein_h_kernel::nops() const
1769 ex Eisenstein_h_kernel::op(size_t i) const
1783 throw (std::out_of_range("Eisenstein_h_kernel::op() out of range"));
1787 ex & Eisenstein_h_kernel::let_op(size_t i)
1789 ensure_if_modifiable();
1803 throw (std::out_of_range("Eisenstein_h_kernel::let_op() out of range"));
1807 bool Eisenstein_h_kernel::is_numeric(void) const
1809 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));
1812 ex Eisenstein_h_kernel::Laurent_series(const ex & x, int order) const
1814 ex res = C_norm * q_expansion_modular_form(x, order+1)/x;
1815 res = res.series(x,order);
1822 * Returns the value of the modular form.
1825 ex Eisenstein_h_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
1827 ex pre = numeric(1)/C_norm;
1829 return get_numerical_value_impl(qbar, pre, 1, N_trunc);
1832 bool Eisenstein_h_kernel::uses_Laurent_series() const
1839 * The constant coefficient in the Fourier expansion.
1842 ex Eisenstein_h_kernel::coefficient_a0(const numeric & k, const numeric & r, const numeric & s, const numeric & N) const
1845 if ( irem(s,N) != 0 ) {
1846 return numeric(1,4) - mod(s,N)/numeric(2)/N;
1848 else if ( (irem(r,N)==0) && (irem(s,N)==0) ) {
1852 return I*numeric(1,4)*cos(Pi*mod(r,N)/N)/sin(Pi*mod(r,N)/N);
1857 return -Bernoulli_polynomial(k,mod(s,N)/N)/numeric(2)/k;
1862 * The higher coefficients in the Fourier expansion.
1865 ex Eisenstein_h_kernel::coefficient_an(const numeric & n, const numeric & k, const numeric & r, const numeric & s, const numeric & N) const
1869 for (numeric m=1; m<=n; m++) {
1870 if ( irem(n,m) == 0 ) {
1871 for (numeric c1=0; c1<N; c1++) {
1873 res += pow(m,k-1)*exp(2*Pi*I/N*mod(r*c2-(s-m)*c1,N)) - pow(-m,k-1)*exp(2*Pi*I/N*mod(-r*c2+(s+m)*c1,N));
1878 return res/numeric(2)/pow(N,k);
1881 ex Eisenstein_h_kernel::q_expansion_modular_form(const ex & q, int N_order) const
1883 numeric N_order_num = numeric(N_order);
1885 numeric k_num = ex_to<numeric>(k);
1886 numeric r_num = ex_to<numeric>(r);
1887 numeric s_num = ex_to<numeric>(s);
1888 numeric N_num = ex_to<numeric>(N);
1890 ex res = coefficient_a0(k_num,r_num,s_num,N_num);
1892 for (numeric i1=1; i1<N_order_num; i1++) {
1893 res += coefficient_an(i1,k_num,r_num,s_num,N_num) * pow(q,i1);
1896 res += Order(pow(q,N_order));
1897 res = res.series(q,N_order);
1902 void Eisenstein_h_kernel::do_print(const print_context & c, unsigned level) const
1904 c.s << "Eisenstein_h_kernel(";
1917 GINAC_BIND_UNARCHIVER(Eisenstein_h_kernel);
1920 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(modular_form_kernel, integration_kernel,
1921 print_func<print_context>(&modular_form_kernel::do_print))
1923 modular_form_kernel::modular_form_kernel() : inherited(), k(_ex0), P(_ex0), C_norm(_ex1)
1927 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)
1931 int modular_form_kernel::compare_same_type(const basic &other) const
1933 const modular_form_kernel &o = static_cast<const modular_form_kernel &>(other);
1936 cmpval = k.compare(o.k);
1941 cmpval = P.compare(o.P);
1946 return C_norm.compare(o.C_norm);
1951 * The series method for this class returns the qbar-expansion of the modular form,
1952 * without an additional factor of C_norm/qbar.
1955 ex modular_form_kernel::series(const relational & r, int order, unsigned options) const
1957 if ( r.rhs() != 0 ) {
1958 throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
1963 subs_q_expansion do_subs_q_expansion(qbar, order);
1965 ex res = do_subs_q_expansion(P).series(qbar,order);
1966 res += Order(pow(qbar,order));
1967 res = res.series(qbar,order);
1972 size_t modular_form_kernel::nops() const
1977 ex modular_form_kernel::op(size_t i) const
1987 throw (std::out_of_range("modular_form_kernel::op() out of range"));
1991 ex & modular_form_kernel::let_op(size_t i)
1993 ensure_if_modifiable();
2003 throw (std::out_of_range("modular_form_kernel::let_op() out of range"));
2007 bool modular_form_kernel::is_numeric(void) const
2009 bool flag = (k.info(info_flags::nonnegint) && C_norm.evalf().info(info_flags::numeric));
2014 symbol qbar("qbar");
2016 // test with a random number and random expansion
2017 return series_to_poly(q_expansion_modular_form(qbar,18)).subs(qbar==numeric(1,937)).evalf().info(info_flags::numeric);
2020 ex modular_form_kernel::Laurent_series(const ex & qbar, int order) const
2022 ex res = series_to_poly(q_expansion_modular_form(qbar,order+1));
2023 res = C_norm * res/qbar;
2024 res = res.series(qbar,order);
2030 * Returns the value of the modular form.
2033 ex modular_form_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
2035 ex pre = numeric(1)/C_norm;
2037 return get_numerical_value_impl(qbar, pre, 1, N_trunc);
2040 bool modular_form_kernel::uses_Laurent_series() const
2045 ex modular_form_kernel::q_expansion_modular_form(const ex & q, int N_order) const
2047 return this->series(q==0,N_order);
2050 void modular_form_kernel::do_print(const print_context & c, unsigned level) const
2052 c.s << "modular_form_kernel(";
2061 GINAC_BIND_UNARCHIVER(modular_form_kernel);
2064 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(user_defined_kernel, integration_kernel,
2065 print_func<print_context>(&user_defined_kernel::do_print))
2067 user_defined_kernel::user_defined_kernel() : inherited(), f(_ex0), x(_ex0)
2071 user_defined_kernel::user_defined_kernel(const ex & arg_f, const ex & arg_x) : inherited(), f(arg_f), x(arg_x)
2075 int user_defined_kernel::compare_same_type(const basic &other) const
2077 const user_defined_kernel &o = static_cast<const user_defined_kernel &>(other);
2080 cmpval = f.compare(o.f);
2085 return x.compare(o.x);
2088 size_t user_defined_kernel::nops() const
2093 ex user_defined_kernel::op(size_t i) const
2101 throw (std::out_of_range("user_defined_kernel::op() out of range"));
2105 ex & user_defined_kernel::let_op(size_t i)
2107 ensure_if_modifiable();
2115 throw (std::out_of_range("user_defined_kernel::let_op() out of range"));
2119 bool user_defined_kernel::is_numeric(void) const
2121 // test with a random number
2122 return f.subs(x==numeric(1,937)).evalf().info(info_flags::numeric);
2125 ex user_defined_kernel::Laurent_series(const ex & x_up, int order) const
2127 ex res = f.series(x,order).subs(x==x_up);
2132 bool user_defined_kernel::uses_Laurent_series() const
2137 void user_defined_kernel::do_print(const print_context & c, unsigned level) const
2139 c.s << "user_defined_kernel(";
2146 GINAC_BIND_UNARCHIVER(user_defined_kernel);
2148 } // namespace GiNaC