Replaced dummy reference by correct reference.
[ginac.git] / ginac / inifcns_elliptic.cpp
1 /** @file inifcns_elliptic.cpp
2  *
3  *  Implementation of some special functions related to elliptic curves
4  *
5  *  The functions are:
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)
9  *
10  *  Some remarks:
11  *
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
14  *
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
18  *
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
22  *
23  */
24
25 /*
26  *  GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
27  *
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.
32  *
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.
37  *
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
41  */
42
43 #include "inifcns.h"
44
45 #include "add.h"
46 #include "constant.h"
47 #include "lst.h"
48 #include "mul.h"
49 #include "numeric.h"
50 #include "operators.h"
51 #include "power.h"
52 #include "pseries.h"
53 #include "relational.h"
54 #include "symbol.h"
55 #include "utils.h"
56 #include "wildcard.h"
57
58 #include "integration_kernel.h"
59 #include "utils_multi_iterator.h"
60
61 #include <cln/cln.h>
62 #include <sstream>
63 #include <stdexcept>
64 #include <vector>
65 #include <cmath>
66
67 namespace GiNaC {
68
69
70 //////////////////////////////////////////////////////////////////////
71 //
72 // Complete elliptic integrals
73 //
74 // helper functions
75 //
76 //////////////////////////////////////////////////////////////////////
77
78
79 // anonymous namespace for helper function
80 namespace {
81
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)
84 {
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));
87         cln::cl_N a_new;
88         cln::cl_N b_new;
89         cln::cl_N res = a_old;
90         cln::cl_N resbuf;
91         do {
92                 resbuf = res;
93
94                 a_new = (a_old+b_old)/2;
95                 b_new = sqrt(a_old*b_old);
96
97                 if ( ( abs(a_new-b_new) > abs(a_new+b_new) ) 
98                      || 
99                      ( (abs(a_new-b_new) == abs(a_new+b_new)) && (imagpart(b_new/a_new) <= 0) ) ) {
100                         b_new *= -1;
101                 }
102
103                 res = a_new;
104                 a_old = a_new;
105                 b_old = b_new;
106     
107         } while (res != resbuf);
108         return res;
109 }
110
111 // Computes
112 //  a0^2 - sum_{n=0}^infinity 2^{n-1}*c_n^2
113 // with
114 //  c_{n+1} = c_n^2/4/a_{n+1}
115 //
116 // Needed for the complete elliptic integral of the second kind.
117 //
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)
119 {
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));
123         cln::cl_N a_new;
124         cln::cl_N b_new;
125         cln::cl_N c_new;
126         cln::cl_N res = square(a_old)-square(c_old)/2;
127         cln::cl_N resbuf;
128         cln::cl_N pre = cln::cl_float(1, cln::float_format(Digits));
129         do {
130                 resbuf = res;
131
132                 a_new = (a_old+b_old)/2;
133                 b_new = sqrt(a_old*b_old);
134
135                 if ( ( abs(a_new-b_new) > abs(a_new+b_new) ) 
136                      || 
137                      ( (abs(a_new-b_new) == abs(a_new+b_new)) && (imagpart(b_new/a_new) <= 0) ) ) {
138                         b_new *= -1;
139                 }
140
141                 c_new = square(c_old)/4/a_new;
142
143                 res -= pre*square(c_new);
144
145                 a_old = a_new;
146                 b_old = b_new;
147                 c_old = c_new;
148                 pre *= 2;
149     
150         } while (res != resbuf);
151         return res;
152 }
153
154
155 } // end of anonymous namespace
156
157
158 //////////////////////////////////////////////////////////////////////
159 //
160 // Complete elliptic integrals
161 //
162 // GiNaC function
163 //
164 //////////////////////////////////////////////////////////////////////
165
166 static ex EllipticK_evalf(const ex& k)
167 {
168         if ( !k.info(info_flags::numeric) ) {
169                 return EllipticK(k).hold();
170         }
171      
172         cln::cl_N kbar = sqrt(1-square(ex_to<numeric>(k).to_cl_N()));
173
174         ex result = Pi/2/numeric(arithmetic_geometric_mean(1,kbar));
175
176         return result.evalf();
177 }
178
179
180 static ex EllipticK_eval(const ex& k)
181 {
182         if (k == _ex0) {
183                 return Pi/2;
184         }
185
186         if ( k.info(info_flags::numeric) && !k.info(info_flags::crational) ) {
187                 return EllipticK(k).evalf();
188         }
189
190         return EllipticK(k).hold();
191 }
192
193
194 static ex EllipticK_deriv(const ex& k, unsigned deriv_param)
195 {
196         return -EllipticK(k)/k + EllipticE(k)/k/(1-k*k);
197 }
198
199
200 static ex EllipticK_series(const ex& k, const relational& rel, int order, unsigned options)
201 {       
202         const ex k_pt = k.subs(rel, subs_options::no_pattern);
203
204         if (k_pt == _ex0) {
205                 const symbol s;
206                 ex ser;
207                 // manually construct the primitive expansion
208                 for (int i=0; i<(order+1)/2; ++i)
209                 {
210                         ser += Pi/2 * numeric(cln::square(cln::binomial(2*i,i))) * pow(s/4,2*i);
211                 }
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);
219         }
220
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!");
223         }
224
225         // all other cases
226         throw do_taylor();
227 }
228
229 static void EllipticK_print_latex(const ex& k, const print_context& c)
230 {
231         c.s << "\\mathrm{K}(";
232         k.print(c);
233         c.s << ")";
234 }
235
236
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());
244
245
246 static ex EllipticE_evalf(const ex& k)
247 {
248         if ( !k.info(info_flags::numeric) ) {
249                 return EllipticE(k).hold();
250         }
251
252         cln::cl_N kbar = sqrt(1-square(ex_to<numeric>(k).to_cl_N()));
253
254         ex result = Pi/2 * numeric( agm_helper_second_kind(1,kbar,ex_to<numeric>(k).to_cl_N()) / arithmetic_geometric_mean(1,kbar) );
255
256         return result.evalf();
257 }
258
259
260 static ex EllipticE_eval(const ex& k)
261 {
262         if (k == _ex0) {
263                 return Pi/2;
264         }
265
266         if ( (k == _ex1) || (k == _ex_1) ) {
267                 return 1;
268         }
269
270         if ( k.info(info_flags::numeric) && !k.info(info_flags::crational) ) {
271                 return EllipticE(k).evalf();
272         }
273
274         return EllipticE(k).hold();
275 }
276
277
278 static ex EllipticE_deriv(const ex& k, unsigned deriv_param)
279 {
280         return -EllipticK(k)/k + EllipticE(k)/k;
281 }
282
283
284 static ex EllipticE_series(const ex& k, const relational& rel, int order, unsigned options)
285 {       
286         const ex k_pt = k.subs(rel, subs_options::no_pattern);
287
288         if (k_pt == _ex0) {
289                 const symbol s;
290                 ex ser;
291                 // manually construct the primitive expansion
292                 for (int i=0; i<(order+1)/2; ++i)
293                 {
294                         ser -= Pi/2 * numeric(cln::square(cln::binomial(2*i,i)))/(2*i-1) * pow(s/4,2*i);
295                 }
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);
303         }
304
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!");
307         }
308
309         // all other cases
310         throw do_taylor();
311 }
312
313 static void EllipticE_print_latex(const ex& k, const print_context& c)
314 {
315         c.s << "\\mathrm{K}(";
316         k.print(c);
317         c.s << ")";
318 }
319
320
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());
328
329
330 //////////////////////////////////////////////////////////////////////
331 //
332 // Iterated integrals
333 //
334 // helper functions
335 //
336 //////////////////////////////////////////////////////////////////////
337
338 // anonymous namespace for helper function
339 namespace {
340
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)
343 {
344         if ( cln::zerop(lambda) ) {
345                 return 0;
346         }
347
348         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
349
350         const int depth = m.size();
351
352         cln::cl_N res = 0;
353         cln::cl_N resbuf;
354         cln::cl_N subexpr;
355
356         if ( N_trunc == 0 ) {
357                 // sum until precision is reached
358                 bool flag_accidental_zero = false;
359
360                 int N = 1;
361
362                 do {
363                         resbuf = res;
364
365                         if ( depth > 1 ) {
366                                 subexpr = 0;
367                                 multi_iterator_ordered_eq<int> i_multi(1,N+1,depth-1);
368                                 for( i_multi.init(); !i_multi.overflow(); i_multi++) {
369                                         cln::cl_N tt = one;
370                                         for (int j=1; j<depth; j++) {
371                                                 if ( j==1 ) {
372                                                         tt = tt * kernel[0]->series_coeff(N-i_multi[depth-2]) / cln::expt(cln::cl_I(i_multi[depth-2]),m[1]);
373                                                 }
374                                                 else {
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]);
376                                                 }
377                                         }
378                                         tt = tt * kernel[depth-1]->series_coeff(i_multi[0]);
379                                         subexpr += tt;
380                                 }
381                         }
382                         else {
383                                 // depth == 1
384                                 subexpr = kernel[0]->series_coeff(N) * one;
385                         }
386                         flag_accidental_zero = cln::zerop(subexpr);
387                         res += cln::expt(lambda, N) / cln::expt(cln::cl_I(N),m[0]) * subexpr;
388                         N++;
389
390                 } while ( (res != resbuf) || flag_accidental_zero );
391         }
392         else {
393                 // N_trunc > 0, sum up the first N_trunc terms
394                 for (int N=1; N<=N_trunc; N++) {
395                         if ( depth > 1 ) {
396                                 subexpr = 0;
397                                 multi_iterator_ordered_eq<int> i_multi(1,N+1,depth-1);
398                                 for( i_multi.init(); !i_multi.overflow(); i_multi++) {
399                                         cln::cl_N tt = one;
400                                         for (int j=1; j<depth; j++) {
401                                                 if ( j==1 ) {
402                                                         tt = tt * kernel[0]->series_coeff(N-i_multi[depth-2]) / cln::expt(cln::cl_I(i_multi[depth-2]),m[1]);
403                                                 }
404                                                 else {
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]);
406                                                 }
407                                         }
408                                         tt = tt * kernel[depth-1]->series_coeff(i_multi[0]);
409                                         subexpr += tt;
410                                 }
411                         }
412                         else {
413                                 // depth == 1
414                                 subexpr = kernel[0]->series_coeff(N) * one;
415                         }
416                         res += cln::expt(lambda, N) / cln::expt(cln::cl_I(N),m[0]) * subexpr;
417                 }
418         }
419
420         return res;
421 }
422
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)
425 {
426         size_t depth = kernel_in.size();
427
428         std::vector<int> m;
429         m.reserve(depth);
430
431         std::vector<const integration_kernel *> kernel;
432         kernel.reserve(depth);
433
434         int n = 1;
435
436         for (const auto & it : kernel_in) {
437                 if ( is_a<basic_log_kernel>(*it) ) {
438                         n++;
439                 }
440                 else {
441                         m.push_back(n);
442                         kernel.push_back( &ex_to<integration_kernel>(*it) );
443                         n = 1;
444                 }
445         }
446
447         cln::cl_N result = iterated_integral_do_sum(m, kernel, lambda, N_trunc);
448
449         return result;
450 }
451
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)
455 {
456         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
457
458         const size_t depth = kernel.size();
459
460         size_t i_trailing = 0;
461         for (size_t i=0; i<depth; i++) {
462                 if ( !(is_a<basic_log_kernel>(*(kernel[i]))) ) {
463                         i_trailing = i+1;
464                 }
465         }
466
467         if ( i_trailing == 0 ) {
468                 return cln::expt(cln::log(lambda), depth) / cln::factorial(depth) * one;
469         }
470
471         if ( i_trailing == depth ) {
472                 return iterated_integral_prepare_m_lst(kernel, lambda, N_trunc);
473         }
474
475         // shuffle
476         std::vector<const integration_kernel *> a,b;
477         for (size_t i=0; i<i_trailing; i++) {
478                 a.push_back(kernel[i]);
479         }
480         for (size_t i=i_trailing; i<depth; i++) {
481                 b.push_back(kernel[i]);
482         }
483
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]);
491                 }
492
493                 result -= iterated_integral_shuffle(new_kernel, lambda, N_trunc);
494         }
495
496         return result;
497 }
498
499 } // end of anonymous namespace
500
501 //////////////////////////////////////////////////////////////////////
502 //
503 // Iterated integrals
504 //
505 // GiNaC function
506 //
507 //////////////////////////////////////////////////////////////////////
508
509 static ex iterated_integral_evalf_impl(const ex& kernel_lst, const ex& lambda, const ex& N_trunc)
510 {
511         // sanity check
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();
514         }
515
516         lst k_lst = ex_to<lst>(kernel_lst);
517
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;
522                 }
523         }
524         if ( flag_not_numeric) {
525                 return iterated_integral(kernel_lst,lambda,N_trunc).hold();
526         }
527
528         for (const auto & it : k_lst) {
529                 if ( !(ex_to<integration_kernel>(it).is_numeric()) ) {
530                         flag_not_numeric = true;
531                 }
532         }
533         if ( flag_not_numeric) {
534                 return iterated_integral(kernel_lst,lambda,N_trunc).hold();
535         }
536
537         // now we know that iterated_integral gives a number
538
539         int N_trunc_int = ex_to<numeric>(N_trunc).to_int();
540
541         // check trailing zeros
542         const size_t depth = k_lst.nops();
543
544         std::vector<const integration_kernel *> kernel_vec;
545         kernel_vec.reserve(depth);
546
547         for (const auto & it : k_lst) {
548                 kernel_vec.push_back( &ex_to<integration_kernel>(it) );
549         }
550
551         size_t i_trailing = 0;
552         for (size_t i=0; i<depth; i++) {
553                 if ( !(kernel_vec[i]->has_trailing_zero()) ) {
554                         i_trailing = i+1;
555                 }
556         }
557
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();
563
564         cln::cl_N result;
565         if ( is_a<basic_log_kernel>(*(kernel_vec[depth-1])) ) {
566                 result = 0;
567         }
568         else {
569                 result = iterated_integral_shuffle(kernel_vec, lambda_cln, N_trunc_int);
570         }
571
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);
578                 }
579                 else {
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);
582                         }
583                 }
584         }
585
586         return numeric(result);
587 }
588
589 static ex iterated_integral2_evalf(const ex& kernel_lst, const ex& lambda)
590 {
591         return iterated_integral_evalf_impl(kernel_lst,lambda,0);
592 }
593
594 static ex iterated_integral3_evalf(const ex& kernel_lst, const ex& lambda, const ex& N_trunc)
595 {
596         return iterated_integral_evalf_impl(kernel_lst,lambda,N_trunc);
597 }
598
599 static ex iterated_integral2_eval(const ex& kernel_lst, const ex& lambda)
600 {
601         if ( lambda.info(info_flags::numeric) && !lambda.info(info_flags::crational) ) {
602                 return iterated_integral(kernel_lst,lambda).evalf();
603         }
604
605         return iterated_integral(kernel_lst,lambda).hold();
606 }
607
608 static ex iterated_integral3_eval(const ex& kernel_lst, const ex& lambda, const ex& N_trunc)
609 {
610         if ( lambda.info(info_flags::numeric) && !lambda.info(info_flags::crational) ) {
611                 return iterated_integral(kernel_lst,lambda,N_trunc).evalf();
612         }
613
614         return iterated_integral(kernel_lst,lambda,N_trunc).hold();
615 }
616
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().
622                                overloaded(2));
623
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().
629                                overloaded(2));
630
631 } // namespace GiNaC
632