1 /** @file inifcns_elliptic.cpp
3 * Implementation of some special functions related to elliptic curves
6 * complete elliptic integral of the first kind EllipticK(k)
7 * complete elliptic integral of the second kind EllipticE(k)
8 * iterated integral iterated_integral(a,y) or iterated_integral(a,y,N_trunc)
12 * - All formulae used can be looked up in the following publication:
13 * [WW] Numerical evaluation of iterated integrals related to elliptic Feynman integrals, M.Walden, S.Weinzierl, arXiv:2010.05271
15 * - When these routines and methods are used for scientific work that leads to publication in a scientific journal,
16 * please refer to this program as :
17 * M.Walden, S.Weinzierl, "Numerical evaluation of iterated integrals related to elliptic Feynman integrals", arXiv:2010.05271
19 * - As these routines build on the core part of GiNaC, it is also polite to acknowledge
20 * C. Bauer, A. Frink, R. Kreckel, "Introduction to the GiNaC Framework for Symbolic Computation within the C++ Programming Language",
21 * J. Symbolic Computations 33, 1 (2002), cs.sc/0004015
26 * GiNaC Copyright (C) 1999-2021 Johannes Gutenberg University Mainz, Germany
28 * This program is free software; you can redistribute it and/or modify
29 * it under the terms of the GNU General Public License as published by
30 * the Free Software Foundation; either version 2 of the License, or
31 * (at your option) any later version.
33 * This program is distributed in the hope that it will be useful,
34 * but WITHOUT ANY WARRANTY; without even the implied warranty of
35 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
36 * GNU General Public License for more details.
38 * You should have received a copy of the GNU General Public License
39 * along with this program; if not, write to the Free Software
40 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
50 #include "operators.h"
53 #include "relational.h"
58 #include "integration_kernel.h"
59 #include "utils_multi_iterator.h"
70 //////////////////////////////////////////////////////////////////////
72 // Complete elliptic integrals
76 //////////////////////////////////////////////////////////////////////
79 // anonymous namespace for helper function
82 // Computes the arithmetic geometric of two numbers a_0 and b_0
83 cln::cl_N arithmetic_geometric_mean(const cln::cl_N & a_0, const cln::cl_N & b_0)
85 cln::cl_N a_old = a_0 * cln::cl_float(1, cln::float_format(Digits));
86 cln::cl_N b_old = b_0 * cln::cl_float(1, cln::float_format(Digits));
89 cln::cl_N res = a_old;
94 a_new = (a_old+b_old)/2;
95 b_new = sqrt(a_old*b_old);
97 if ( ( abs(a_new-b_new) > abs(a_new+b_new) )
99 ( (abs(a_new-b_new) == abs(a_new+b_new)) && (imagpart(b_new/a_new) <= 0) ) ) {
107 } while (res != resbuf);
112 // a0^2 - sum_{n=0}^infinity 2^{n-1}*c_n^2
114 // c_{n+1} = c_n^2/4/a_{n+1}
116 // Needed for the complete elliptic integral of the second kind.
118 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)
120 cln::cl_N a_old = a_0 * cln::cl_float(1, cln::float_format(Digits));
121 cln::cl_N b_old = b_0 * cln::cl_float(1, cln::float_format(Digits));
122 cln::cl_N c_old = c_0 * cln::cl_float(1, cln::float_format(Digits));
126 cln::cl_N res = square(a_old)-square(c_old)/2;
128 cln::cl_N pre = cln::cl_float(1, cln::float_format(Digits));
132 a_new = (a_old+b_old)/2;
133 b_new = sqrt(a_old*b_old);
135 if ( ( abs(a_new-b_new) > abs(a_new+b_new) )
137 ( (abs(a_new-b_new) == abs(a_new+b_new)) && (imagpart(b_new/a_new) <= 0) ) ) {
141 c_new = square(c_old)/4/a_new;
143 res -= pre*square(c_new);
150 } while (res != resbuf);
155 } // end of anonymous namespace
158 //////////////////////////////////////////////////////////////////////
160 // Complete elliptic integrals
164 //////////////////////////////////////////////////////////////////////
166 static ex EllipticK_evalf(const ex& k)
168 if ( !k.info(info_flags::numeric) ) {
169 return EllipticK(k).hold();
172 cln::cl_N kbar = sqrt(1-square(ex_to<numeric>(k).to_cl_N()));
174 ex result = Pi/2/numeric(arithmetic_geometric_mean(1,kbar));
176 return result.evalf();
180 static ex EllipticK_eval(const ex& k)
186 if ( k.info(info_flags::numeric) && !k.info(info_flags::crational) ) {
187 return EllipticK(k).evalf();
190 return EllipticK(k).hold();
194 static ex EllipticK_deriv(const ex& k, unsigned deriv_param)
196 return -EllipticK(k)/k + EllipticE(k)/k/(1-k*k);
200 static ex EllipticK_series(const ex& k, const relational& rel, int order, unsigned options)
202 const ex k_pt = k.subs(rel, subs_options::no_pattern);
207 // manually construct the primitive expansion
208 for (int i=0; i<(order+1)/2; ++i)
210 ser += Pi/2 * numeric(cln::square(cln::binomial(2*i,i))) * pow(s/4,2*i);
212 // substitute the argument's series expansion
213 ser = ser.subs(s==k.series(rel, order), subs_options::no_pattern);
214 // maybe that was terminating, so add a proper order term
215 epvector nseq { expair(Order(_ex1), order) };
216 ser += pseries(rel, std::move(nseq));
217 // reexpanding it will collapse the series again
218 return ser.series(rel, order);
221 if ( (k_pt == _ex1) || (k_pt == _ex_1) ) {
222 throw std::runtime_error("EllipticK_series: don't know how to do the series expansion at this point!");
229 static void EllipticK_print_latex(const ex& k, const print_context& c)
231 c.s << "\\mathrm{K}(";
237 REGISTER_FUNCTION(EllipticK,
238 evalf_func(EllipticK_evalf).
239 eval_func(EllipticK_eval).
240 derivative_func(EllipticK_deriv).
241 series_func(EllipticK_series).
242 print_func<print_latex>(EllipticK_print_latex).
243 do_not_evalf_params());
246 static ex EllipticE_evalf(const ex& k)
248 if ( !k.info(info_flags::numeric) ) {
249 return EllipticE(k).hold();
252 cln::cl_N kbar = sqrt(1-square(ex_to<numeric>(k).to_cl_N()));
254 ex result = Pi/2 * numeric( agm_helper_second_kind(1,kbar,ex_to<numeric>(k).to_cl_N()) / arithmetic_geometric_mean(1,kbar) );
256 return result.evalf();
260 static ex EllipticE_eval(const ex& k)
266 if ( (k == _ex1) || (k == _ex_1) ) {
270 if ( k.info(info_flags::numeric) && !k.info(info_flags::crational) ) {
271 return EllipticE(k).evalf();
274 return EllipticE(k).hold();
278 static ex EllipticE_deriv(const ex& k, unsigned deriv_param)
280 return -EllipticK(k)/k + EllipticE(k)/k;
284 static ex EllipticE_series(const ex& k, const relational& rel, int order, unsigned options)
286 const ex k_pt = k.subs(rel, subs_options::no_pattern);
291 // manually construct the primitive expansion
292 for (int i=0; i<(order+1)/2; ++i)
294 ser -= Pi/2 * numeric(cln::square(cln::binomial(2*i,i)))/(2*i-1) * pow(s/4,2*i);
296 // substitute the argument's series expansion
297 ser = ser.subs(s==k.series(rel, order), subs_options::no_pattern);
298 // maybe that was terminating, so add a proper order term
299 epvector nseq { expair(Order(_ex1), order) };
300 ser += pseries(rel, std::move(nseq));
301 // reexpanding it will collapse the series again
302 return ser.series(rel, order);
305 if ( (k_pt == _ex1) || (k_pt == _ex_1) ) {
306 throw std::runtime_error("EllipticE_series: don't know how to do the series expansion at this point!");
313 static void EllipticE_print_latex(const ex& k, const print_context& c)
315 c.s << "\\mathrm{K}(";
321 REGISTER_FUNCTION(EllipticE,
322 evalf_func(EllipticE_evalf).
323 eval_func(EllipticE_eval).
324 derivative_func(EllipticE_deriv).
325 series_func(EllipticE_series).
326 print_func<print_latex>(EllipticE_print_latex).
327 do_not_evalf_params());
330 //////////////////////////////////////////////////////////////////////
332 // Iterated integrals
336 //////////////////////////////////////////////////////////////////////
338 // anonymous namespace for helper function
341 // performs the actual series summation for an iterated integral
342 cln::cl_N iterated_integral_do_sum(const std::vector<int> & m, const std::vector<const integration_kernel *> & kernel, const cln::cl_N & lambda, int N_trunc)
344 if ( cln::zerop(lambda) ) {
348 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
350 const int depth = m.size();
356 if ( N_trunc == 0 ) {
357 // sum until precision is reached
358 bool flag_accidental_zero = false;
367 multi_iterator_ordered_eq<int> i_multi(1,N+1,depth-1);
368 for( i_multi.init(); !i_multi.overflow(); i_multi++) {
370 for (int j=1; j<depth; j++) {
372 tt = tt * kernel[0]->series_coeff(N-i_multi[depth-2]) / cln::expt(cln::cl_I(i_multi[depth-2]),m[1]);
375 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]);
378 tt = tt * kernel[depth-1]->series_coeff(i_multi[0]);
384 subexpr = kernel[0]->series_coeff(N) * one;
386 flag_accidental_zero = cln::zerop(subexpr);
387 res += cln::expt(lambda, N) / cln::expt(cln::cl_I(N),m[0]) * subexpr;
390 } while ( (res != resbuf) || flag_accidental_zero );
393 // N_trunc > 0, sum up the first N_trunc terms
394 for (int N=1; N<=N_trunc; N++) {
397 multi_iterator_ordered_eq<int> i_multi(1,N+1,depth-1);
398 for( i_multi.init(); !i_multi.overflow(); i_multi++) {
400 for (int j=1; j<depth; j++) {
402 tt = tt * kernel[0]->series_coeff(N-i_multi[depth-2]) / cln::expt(cln::cl_I(i_multi[depth-2]),m[1]);
405 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]);
408 tt = tt * kernel[depth-1]->series_coeff(i_multi[0]);
414 subexpr = kernel[0]->series_coeff(N) * one;
416 res += cln::expt(lambda, N) / cln::expt(cln::cl_I(N),m[0]) * subexpr;
423 // figure out the number of basic_log_kernels before a non-basic_log_kernel
424 cln::cl_N iterated_integral_prepare_m_lst(const std::vector<const integration_kernel *> & kernel_in, const cln::cl_N & lambda, int N_trunc)
426 size_t depth = kernel_in.size();
431 std::vector<const integration_kernel *> kernel;
432 kernel.reserve(depth);
436 for (const auto & it : kernel_in) {
437 if ( is_a<basic_log_kernel>(*it) ) {
442 kernel.push_back( &ex_to<integration_kernel>(*it) );
447 cln::cl_N result = iterated_integral_do_sum(m, kernel, lambda, N_trunc);
452 // shuffle to remove trailing zeros,
453 // integration kernels, which are not basic_log_kernels, are treated as regularised kernels
454 cln::cl_N iterated_integral_shuffle(const std::vector<const integration_kernel *> & kernel, const cln::cl_N & lambda, int N_trunc)
456 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
458 const size_t depth = kernel.size();
460 size_t i_trailing = 0;
461 for (size_t i=0; i<depth; i++) {
462 if ( !(is_a<basic_log_kernel>(*(kernel[i]))) ) {
467 if ( i_trailing == 0 ) {
468 return cln::expt(cln::log(lambda), depth) / cln::factorial(depth) * one;
471 if ( i_trailing == depth ) {
472 return iterated_integral_prepare_m_lst(kernel, lambda, N_trunc);
476 std::vector<const integration_kernel *> a,b;
477 for (size_t i=0; i<i_trailing; i++) {
478 a.push_back(kernel[i]);
480 for (size_t i=i_trailing; i<depth; i++) {
481 b.push_back(kernel[i]);
484 cln::cl_N result = iterated_integral_prepare_m_lst(a, lambda, N_trunc) * cln::expt(cln::log(lambda), depth-i_trailing) / cln::factorial(depth-i_trailing);
485 multi_iterator_shuffle_prime<const integration_kernel *> i_shuffle(a,b);
486 for( i_shuffle.init(); !i_shuffle.overflow(); i_shuffle++) {
487 std::vector<const integration_kernel *> new_kernel;
488 new_kernel.reserve(depth);
489 for (size_t i=0; i<depth; i++) {
490 new_kernel.push_back(i_shuffle[i]);
493 result -= iterated_integral_shuffle(new_kernel, lambda, N_trunc);
499 } // end of anonymous namespace
501 //////////////////////////////////////////////////////////////////////
503 // Iterated integrals
507 //////////////////////////////////////////////////////////////////////
509 static ex iterated_integral_evalf_impl(const ex& kernel_lst, const ex& lambda, const ex& N_trunc)
512 if ((!kernel_lst.info(info_flags::list)) || (!lambda.evalf().info(info_flags::numeric)) || (!N_trunc.info(info_flags::nonnegint))) {
513 return iterated_integral(kernel_lst,lambda,N_trunc).hold();
516 lst k_lst = ex_to<lst>(kernel_lst);
518 bool flag_not_numeric = false;
519 for (const auto & it : k_lst) {
520 if ( !is_a<integration_kernel>(it) ) {
521 flag_not_numeric = true;
524 if ( flag_not_numeric) {
525 return iterated_integral(kernel_lst,lambda,N_trunc).hold();
528 for (const auto & it : k_lst) {
529 if ( !(ex_to<integration_kernel>(it).is_numeric()) ) {
530 flag_not_numeric = true;
533 if ( flag_not_numeric) {
534 return iterated_integral(kernel_lst,lambda,N_trunc).hold();
537 // now we know that iterated_integral gives a number
539 int N_trunc_int = ex_to<numeric>(N_trunc).to_int();
541 // check trailing zeros
542 const size_t depth = k_lst.nops();
544 std::vector<const integration_kernel *> kernel_vec;
545 kernel_vec.reserve(depth);
547 for (const auto & it : k_lst) {
548 kernel_vec.push_back( &ex_to<integration_kernel>(it) );
551 size_t i_trailing = 0;
552 for (size_t i=0; i<depth; i++) {
553 if ( !(kernel_vec[i]->has_trailing_zero()) ) {
558 // split integral into regularised integrals and trailing basic_log_kernels
559 // non-basic_log_kernels are treated as regularised kernels in call to iterated_integral_shuffle
560 cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
561 cln::cl_N lambda_cln = ex_to<numeric>(lambda.evalf()).to_cl_N();
562 basic_log_kernel L0 = basic_log_kernel();
565 if ( is_a<basic_log_kernel>(*(kernel_vec[depth-1])) ) {
569 result = iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
572 cln::cl_N coeff = one;
573 for (size_t i_plus=depth; i_plus>i_trailing; i_plus--) {
574 coeff = coeff * kernel_vec[i_plus-1]->series_coeff(0);
575 kernel_vec[i_plus-1] = &L0;
576 if ( i_plus==i_trailing+1 ) {
577 result += coeff * iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
580 if ( !(is_a<basic_log_kernel>(*(kernel_vec[i_plus-2]))) ) {
581 result += coeff * iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
586 return numeric(result);
589 static ex iterated_integral2_evalf(const ex& kernel_lst, const ex& lambda)
591 return iterated_integral_evalf_impl(kernel_lst,lambda,0);
594 static ex iterated_integral3_evalf(const ex& kernel_lst, const ex& lambda, const ex& N_trunc)
596 return iterated_integral_evalf_impl(kernel_lst,lambda,N_trunc);
599 static ex iterated_integral2_eval(const ex& kernel_lst, const ex& lambda)
601 if ( lambda.info(info_flags::numeric) && !lambda.info(info_flags::crational) ) {
602 return iterated_integral(kernel_lst,lambda).evalf();
605 return iterated_integral(kernel_lst,lambda).hold();
608 static ex iterated_integral3_eval(const ex& kernel_lst, const ex& lambda, const ex& N_trunc)
610 if ( lambda.info(info_flags::numeric) && !lambda.info(info_flags::crational) ) {
611 return iterated_integral(kernel_lst,lambda,N_trunc).evalf();
614 return iterated_integral(kernel_lst,lambda,N_trunc).hold();
617 unsigned iterated_integral2_SERIAL::serial =
618 function::register_new(function_options("iterated_integral", 2).
619 eval_func(iterated_integral2_eval).
620 evalf_func(iterated_integral2_evalf).
621 do_not_evalf_params().
624 unsigned iterated_integral3_SERIAL::serial =
625 function::register_new(function_options("iterated_integral", 3).
626 eval_func(iterated_integral3_eval).
627 evalf_func(iterated_integral3_evalf).
628 do_not_evalf_params().