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-2021 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