ginac/integration_kernel.cpp

   1 /** @file integration_kernel.cpp
   2  *
   3  *  Implementation of GiNaC's integration kernels for iterated integrals. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2020 Johannes Gutenberg University Mainz, Germany
   7  *
   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.
  12  *
  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.
  17  *
  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
  21  */
  22
  23 #include "integration_kernel.h"
  24 #include "add.h"
  25 #include "mul.h"
  26 #include "operators.h"
  27 #include "power.h"
  28 #include "relational.h"
  29 #include "symbol.h"
  30 #include "constant.h"
  31 #include "numeric.h"
  32 #include "function.h"
  33 #include "pseries.h"
  34 #include "utils.h"
  35 #include "inifcns.h"
  36
  37 #include <iostream>
  38 #include <stdexcept>
  39 #include <cln/cln.h>
  40
  41
  42 namespace GiNaC {
  43
  44 // anonymous namespace for helper function
  45 namespace {
  46
  47 /**
  48  *
  49  * Returns the Kronecker symbol in the case where
  50  *  a: integer
  51  *  n: unit or prime number
  52  *
  53  * If n is an odd prime number, the routine returns the Legendre symbol.
  54  *
  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.
  57  *
  58  * Ref.: Toshitsune Miyake, Modular Forms, Chapter 3.1
  59  *
  60  */
  61 numeric kronecker_symbol_prime(const numeric & a, const numeric & n)
  62 {
  63         if ( n == 1 ) {
  64                 return 1;
  65         }
  66         else if ( n == -1 ) {
  67                 if ( a >= 0 ) {
  68                         return 1;
  69                 }
  70                 else {
  71                         return -1;
  72                 }
  73         }
  74         else if ( n == 2 ) {
  75                 if ( GiNaC::smod(a,8) == 1 ) {
  76                         return 1;
  77                 }
  78                 else if ( GiNaC::smod(a,8) == -1 ) {
  79                         return 1;
  80                 }
  81                 else if  ( GiNaC::smod(a,8) == 3 ) {
  82                         return -1;
  83                 }
  84                 else if ( GiNaC::smod(a,8) == -3 ) {
  85                         return -1;
  86                 }
  87                 else {
  88                         return 0;
  89                 }
  90         }
  91
  92         // n is an odd prime number
  93         return GiNaC::smod( pow(a,(n-1)/2), n);
  94 }
  95
  96 /**
  97  *
  98  * n:     positive integer
  99  *
 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
 104  *
 105  * k:     modular weight
 106  *
 107  * This function computes
 108  * \f[
 109  *      \sum\limits_{d | n} \psi(d) \phi(n/d) d^{k-1}
 110  * \f]
 111  *
 112  * Ref.: Eq.(5.3.1), William A. Stein, Modular Forms, A computational Approach;
 113  *
 114  *       Eq.(32), arXiv:1704.08895
 115  *
 116  */
 117 numeric divisor_function(const numeric & n, const numeric & a, const numeric & b, const numeric & k)
 118 {
 119         ex res = 0;
 120
 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);
 125                 }
 126         }
 127
 128         return ex_to<numeric>(res);
 129 }
 130
 131 /**
 132  *
 133  * k:     modular weight
 134  *
 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
 139  *
 140  * This function computes the constant term of the q-expansion of an Eisenstein series.
 141  *
 142  * The coefficient is given by the equation below eq.(5.3.1) in William A. Stein, Modular Forms, A computational Approach;
 143  *
 144  * or by eq.(33), arXiv:1704.08895
 145  *
 146  */
 147 numeric coefficient_a0(const numeric & k, const numeric & a, const numeric & b)
 148 {
 149         ex conductor = abs(a);
 150
 151         numeric a0;
 152         if ( conductor == 1 ) {
 153                 a0 = -numeric(1,2)/k*generalised_Bernoulli_number(k,b);
 154         }
 155         else {
 156                 a0 = 0;
 157         }
 158
 159         return a0;
 160 }
 161
 162 /**
 163  *
 164  * k:     modular weight
 165  * q:     exp(2 Pi I tau/M)
 166  *
 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
 171  *
 172  * N:     truncation order
 173  *
 174  * Returns the q-expansion of an Eisenstein to order N (the q^(N-1)-term is included, q^N is neglected).
 175  *
 176  * Ref.: Eq.(5.3.1), William A. Stein, Modular Forms, A computational Approach;
 177  *
 178  *       Eq.(32), arXiv:1704.08895
 179  *
 180  */
 181 ex eisenstein_series(const numeric & k, const ex & q, const numeric & a, const numeric & b, const numeric & N)
 182 {
 183         ex res = coefficient_a0(k,a,b);
 184
 185         for (numeric i1=1; i1<N; i1++) {
 186                 res += divisor_function(i1,a,b,k) * pow(q,i1);
 187         }
 188
 189         return res;
 190 }
 191
 192 /**
 193  *
 194  * Returns the q_N-expansion of the Eisenstein series E_{k,a,b}(K tau_N)
 195  *
 196  */
 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)
 198 {
 199         int N_order_int = N_order.to_int();
 200
 201         ex res = eisenstein_series(k,pow(q,K),a,b,iquo(N_order,K));
 202
 203         res += Order(pow(q,N_order_int));
 204         res = res.series(q,N_order_int);
 205
 206         return res;
 207 }
 208
 209 /**
 210  *
 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).
 213  *
 214  */
 215 ex B_eisenstein_series(const ex & q, const numeric & N_level, const numeric & K, const numeric & N_order)
 216 {
 217         int N_order_int = N_order.to_int();
 218
 219         ex res = eisenstein_series(2,q,1,1,N_order) - K*eisenstein_series(2,pow(q,K),1,1,iquo(N_order,K));
 220
 221         res += Order(pow(q,N_order_int));
 222         res = res.series(q,N_order_int);
 223
 224         return res;
 225 }
 226
 227 /**
 228  *
 229  * A helper function to expand polynomials in Eisenstein series.
 230  *
 231  */
 232 struct subs_q_expansion : public map_function
 233 {
 234         subs_q_expansion(const ex & arg_qbar, int arg_order) : qbar(arg_qbar), order(arg_order)
 235                 {}
 236
 237         ex operator()(const ex & e)
 238                 {
 239                         if ( is_a<Eisenstein_kernel>(e) || is_a<Eisenstein_h_kernel>(e) ) {
 240                                 return series_to_poly(e.series(qbar,order));
 241                         }
 242                         else {
 243                                 return e.map(*this);
 244                         }
 245                 }
 246
 247         ex qbar;
 248         int order;
 249 };
 250
 251 /**
 252  *
 253  * \f[
 254  *     Li_{-n}(x),  n>=0
 255  * \f]
 256  *
 257  * This is a rational function in x.
 258  *
 259  * To speed things up, we cache it.
 260  *
 261  */
 262 class Li_negative
 263 {
 264         // ctors
 265 public:
 266         Li_negative();
 267
 268         // non-virtual functions
 269 public:
 270         ex get_symbolic_value(int n, const ex & x_val);
 271         ex get_numerical_value(int n, const ex & x_val);
 272
 273         // member variables :
 274 protected:
 275         static std::vector<ex> cache_vec;
 276         static symbol x;
 277 };
 278
 279
 280 Li_negative::Li_negative() {}
 281
 282 ex Li_negative::get_symbolic_value(int n, const ex & x_val)
 283 {
 284         int n_cache = cache_vec.size();
 285
 286         if ( n >= n_cache ) {
 287                 for (int j=n_cache; j<=n; j++) {
 288                         ex f;
 289                         if ( j == 0 ) {
 290                                 f = x/(1-x);
 291                         }
 292                         else {
 293                                 f = normal( x*diff(cache_vec[j-1],x));
 294                         }
 295                         cache_vec.push_back( f );
 296                 }
 297         }
 298
 299         return cache_vec[n].subs(x==x_val);
 300 }
 301
 302 ex Li_negative::get_numerical_value(int n, const ex & x_val)
 303 {
 304         symbol x_symb("x_symb");
 305
 306         ex f = this->get_symbolic_value(n,x_symb);
 307
 308         ex res = f.subs(x_symb==x_val).evalf();
 309
 310         return res;
 311 }
 312
 313 // initialise static data members
 314 std::vector<ex> Li_negative::cache_vec;
 315 symbol Li_negative::x = symbol("x");
 316
 317
 318 } // end of anonymous namespace
 319
 320 /**
 321  *
 322  * Returns the decomposition of the positive integer n into prime numbers
 323  * in the form
 324  *  lst( lst(p1,...,pr), lst(a1,...,ar) )
 325  * such that
 326  *  n = p1^a1 ... pr^ar.
 327  *
 328  */
 329 ex ifactor(const numeric & n)
 330 {
 331         if ( !n.is_pos_integer() ) throw (std::runtime_error("ifactor(): argument not a positive integer"));
 332
 333         lst p_lst, exp_lst;
 334
 335         // implementation for small integers
 336         numeric n_temp = n;
 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 ) {
 341                                 n_temp = n_temp/p;
 342                                 exp_temp++;
 343                         }
 344                         if ( exp_temp>0 ) {
 345                                 p_lst.append(p);
 346                                 exp_lst.append(exp_temp);
 347                         }
 348                 }
 349                 if ( n_temp == 1 ) break;
 350         }
 351
 352         if ( n_temp != 1 ) throw (std::runtime_error("ifactor(): probabilistic primality test failed"));
 353
 354         lst res = {p_lst,exp_lst};
 355
 356         return res;
 357 }
 358
 359 /**
 360  *
 361  * Returns true if the integer n is either one or the discriminant of a quadratic number field.
 362  *
 363  * Returns false otherwise.
 364  *
 365  * Ref.: Toshitsune Miyake, Modular Forms, Chapter 3.1
 366  *
 367  */
 368 bool is_discriminant_of_quadratic_number_field(const numeric & n)
 369 {
 370         if ( n == 0 ) {
 371                 return false;
 372         }
 373
 374         if ( n == 1 ) {
 375                 return true;
 376         }
 377
 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));
 381
 382         size_t n_primes = p_lst.nops();
 383
 384         if ( n_primes > 0 ) {
 385                 // take the last prime
 386                 numeric p = ex_to<numeric>(p_lst.op(n_primes-1));
 387
 388                 if ( p.is_odd() ) {
 389                         if ( e_lst.op(n_primes-1) != 1 ) {
 390                                 return false;
 391                         }
 392
 393                         numeric pstar = p;
 394                         if ( mod(p,4) == 3 ) {
 395                                 pstar = -p;
 396                         }
 397                         return is_discriminant_of_quadratic_number_field(n/pstar);
 398                 }
 399         }
 400         // power of two now
 401         if ( (n==-4) || (n==-8) || (n==8) || (n==-32) || (n==32) || (n==-64) || (n==128) ) {
 402                 return true;
 403         }
 404
 405         return false;
 406 }
 407
 408 /**
 409  *
 410  * Returns the Kronecker symbol
 411  *  a: integer
 412  *  n: integer
 413  *
 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
 418  *
 419  * In particular
 420  *  kronecker_symbol(-1,0) = 1 (in agreement with Sage)
 421  *
 422  * Ref.: Toshitsune Miyake, Modular Forms, Chapter 3.1
 423  *
 424  */
 425 numeric kronecker_symbol(const numeric & a, const numeric & n)
 426 {
 427         // case n=0 first, include kronecker_symbol(0,0)=0
 428         if ( n == 0 ) {
 429                 if ( (a == 1) || (a == -1) ) {
 430                         return 1;
 431                 }
 432                 else {
 433                         return 0;
 434                 }
 435         }
 436
 437         numeric unit = 1;
 438         numeric n_pos = n;
 439         if ( n_pos<0 ) {
 440                 unit = -1;
 441                 n_pos = -n;
 442         }
 443
 444         ex res = kronecker_symbol_prime(a,unit);
 445
 446         numeric n_odd = n_pos;
 447         numeric alpha = 0;
 448         while ( n_odd.is_even() ) {
 449                 alpha++;
 450                 n_odd = n_odd/2;
 451         }
 452         if ( alpha>0 ) {
 453                 res *= pow(kronecker_symbol_prime(a,2),alpha);
 454         }
 455
 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));
 459
 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));
 462         }
 463
 464         return ex_to<numeric>(res);
 465 }
 466
 467 /**
 468  *
 469  * Defines a primitive Dirichlet character through the Kronecker symbol.
 470  *
 471  *  n:  integer
 472  *  a:  discriminant of a quadratic field
 473  * |a|: conductor
 474  *
 475  * The character takes the values -1,0,1.
 476  *
 477  */
 478 numeric primitive_dirichlet_character(const numeric & n, const numeric & a)
 479 {
 480         return kronecker_symbol(a,n);
 481 }
 482
 483 /**
 484  *
 485  * Defines a Dirichlet character through the Kronecker symbol.
 486  *
 487  *  n:  integer
 488  *  a:  discriminant of a quadratic field
 489  * |a|: conductor
 490  *  N:  modulus, needs to be multiple of |a|
 491  *
 492  * The character takes the values -1,0,1.
 493  *
 494  */
 495 numeric dirichlet_character(const numeric & n, const numeric & a, const numeric & N)
 496 {
 497         if ( gcd(n,N) == 1 ) {
 498                 return primitive_dirichlet_character(n,a);
 499         }
 500
 501         return 0;
 502 }
 503
 504 /**
 505  *
 506  * The generalised Bernoulli number.
 507  *
 508  * k:     index / modular weight
 509  *
 510  * b:     discriminant of a quadratic field, defines primitive character psi
 511  * M=|b|: conductor of primitive character psi
 512  *
 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
 515  *
 516  */
 517 numeric generalised_Bernoulli_number(const numeric & k, const numeric & b)
 518 {
 519         int k_int = k.to_int();
 520
 521         symbol x("x");
 522
 523         numeric conductor = abs(b);
 524
 525         ex gen_fct = 0;
 526         for (numeric i1=1; i1<=conductor; i1++) {
 527                 gen_fct += primitive_dirichlet_character(i1,b) * x*exp(i1*x)/(exp(conductor*x)-1);
 528         }
 529
 530         gen_fct = series_to_poly(gen_fct.series(x,k_int+1));
 531
 532         ex B = factorial(k) * gen_fct.coeff(x,k_int);
 533
 534         return ex_to<numeric>(B);
 535 }
 536
 537 /**
 538  *
 539  * The Bernoulli polynomials
 540  *
 541  */
 542 ex Bernoulli_polynomial(const numeric & k, const ex & x)
 543 {
 544         int k_int = k.to_int();
 545
 546         symbol t("t");
 547
 548         ex gen_fct = t*exp(x*t)/(exp(t)-1);
 549
 550         gen_fct = series_to_poly(gen_fct.series(t,k_int+1));
 551
 552         ex B = factorial(k) * gen_fct.coeff(t,k_int);
 553
 554         return B;
 555 }
 556
 557
 558
 559 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integration_kernel, basic,
 560   print_func<print_context>(&integration_kernel::do_print))
 561
 562 integration_kernel::integration_kernel() : inherited(), cache_step_size(100), series_vec()
 563 {
 564 }
 565
 566 int integration_kernel::compare_same_type(const basic &other) const
 567 {
 568         return 0;
 569 }
 570
 571 ex integration_kernel::series(const relational & r, int order, unsigned options) const
 572 {
 573         if ( r.rhs() != 0 ) {
 574                 throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
 575         }
 576
 577         return Laurent_series(r.lhs(),order);
 578 }
 579
 580 /**
 581  *
 582  * This routine returns true, if the integration kernel has a trailing zero.
 583  *
 584  */
 585 bool integration_kernel::has_trailing_zero(void) const
 586 {
 587         if ( cln::zerop( series_coeff(0) ) ) {
 588                 return false;
 589         }
 590
 591         return true;
 592 }
 593
 594 /**
 595  *
 596  * This routine returns true, if the integration kernel can be evaluated numerically.
 597  *
 598  */
 599 bool integration_kernel::is_numeric(void) const
 600 {
 601         return true;
 602 }
 603
 604 /**
 605  *
 606  * Subclasses have either to implement series_coeff_impl
 607  * or the two methods Laurent_series and uses_Laurent_series.
 608  *
 609  * The method series_coeff_impl can be used, if a single coefficient can be computed
 610  * independently of the others.
 611  *
 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.
 614  *
 615  */
 616 cln::cl_N integration_kernel::series_coeff(int i) const
 617 {
 618         int n_vec = series_vec.size();
 619
 620         if ( i >= n_vec ) {
 621                 int N = cache_step_size*(i/cache_step_size+1);
 622
 623                 if ( uses_Laurent_series() ) {
 624                         symbol x("x");
 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() );
 630                         }
 631                 }
 632                 else {
 633                         for (int j=n_vec; j<N; j++) {
 634                                 series_vec.push_back( series_coeff_impl(j) );
 635                         }
 636                 }
 637         }
 638
 639         return series_vec[i];
 640 }
 641
 642 /**
 643  *
 644  * For \f$ \omega = d\lambda \f$ only the coefficient of \f$ \lambda^0 \f$ is non-zero.
 645  *
 646  * The i-th coefficient corresponds to the power \f$ \lambda^{i-1} \f$.
 647  *
 648  */
 649 cln::cl_N integration_kernel::series_coeff_impl(int i) const
 650 {
 651         if ( i == 1 ) {
 652                 return 1;
 653         }
 654
 655         return 0;
 656 }
 657
 658 /**
 659  *
 660  * Returns the Laurent series, starting possibly with the pole term.
 661  * Neglected terms are of order \f$ O(x^order) \f$.
 662  *
 663  */
 664 ex integration_kernel::Laurent_series(const ex & x, int order) const
 665 {
 666         ex res = 0;
 667         for (int n=-1; n<order; n++) {
 668                 res += numeric(series_coeff(n+1)) * pow(x,n);
 669         }
 670         res += Order(pow(x,order));
 671         res = res.series(x,order);
 672
 673         return res;
 674 }
 675
 676 /**
 677  *
 678  * Evaluates the integrand at lambda.
 679  *
 680  */
 681 ex  integration_kernel::get_numerical_value(const ex & lambda, int N_trunc) const
 682 {
 683         return get_numerical_value_impl(lambda, 1, 0, N_trunc);
 684 }
 685
 686 /**
 687  *
 688  * Returns true, if the coefficients are computed from the Laurent series
 689  * (in which case the method Laurent_series needs to be implemented).
 690  *
 691  * Returns false if this is not the case
 692  * (and the class has an implementation of series_coeff_impl).
 693  *
 694  */
 695 bool integration_kernel::uses_Laurent_series() const
 696 {
 697         return false;
 698 }
 699
 700 /**
 701  *
 702  * Returns the current size of the cache.
 703  *
 704  */
 705 size_t integration_kernel::get_cache_size(void) const
 706 {
 707         return series_vec.size();
 708 }
 709
 710 /**
 711  *
 712  * Sets the step size by which the cache is increased.
 713  *
 714  */
 715 void integration_kernel::set_cache_step(int cache_steps) const
 716 {
 717         cache_step_size = cache_steps;
 718 }
 719
 720 /**
 721  *
 722  * Wrapper around series_coeff(i), converts cl_N to numeric.
 723  *
 724  */
 725 ex integration_kernel::get_series_coeff(int i) const
 726 {
 727         return numeric(series_coeff(i));
 728 }
 729
 730 /**
 731  *
 732  * The actual implementation for computing a numerical value for the integrand.
 733  *
 734  */
 735 ex  integration_kernel::get_numerical_value_impl(const ex & lambda, const ex & pre, int shift, int N_trunc) const
 736 {
 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();
 739
 740         cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
 741
 742         cln::cl_N res = 0;
 743         cln::cl_N resbuf;
 744         cln::cl_N subexpr;
 745
 746         if ( N_trunc == 0 ) {
 747                 // sum until precision is reached
 748                 bool flag_accidental_zero = false;
 749
 750                 int N = 0;
 751
 752                 do {
 753                         resbuf = res;
 754
 755                         subexpr = series_coeff(N);
 756
 757                         res += pre_cln * subexpr * cln::expt(lambda_cln,N-1+shift);
 758
 759                         flag_accidental_zero = cln::zerop(subexpr);
 760
 761                         N++;
 762                 } while ( (res != resbuf) || flag_accidental_zero );
 763         }
 764         else {
 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);
 768
 769                         res += pre_cln * subexpr * cln::expt(lambda_cln,N-1+shift);
 770                 }
 771         }
 772
 773         return numeric(res);
 774 }
 775
 776 void integration_kernel::do_print(const print_context & c, unsigned level) const
 777 {
 778         c.s << "integration_kernel()";
 779 }
 780
 781 GINAC_BIND_UNARCHIVER(integration_kernel);
 782
 783
 784 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic_log_kernel, integration_kernel,
 785                                      print_func<print_context>(&basic_log_kernel::do_print))
 786
 787 basic_log_kernel::basic_log_kernel() : inherited()
 788 {
 789 }
 790
 791 int basic_log_kernel::compare_same_type(const basic &other) const
 792 {
 793         return 0;
 794 }
 795
 796 cln::cl_N basic_log_kernel::series_coeff_impl(int i) const
 797 {
 798         if ( i == 0 ) {
 799                 return 1;
 800         }
 801
 802         return 0;
 803 }
 804
 805 void basic_log_kernel::do_print(const print_context & c, unsigned level) const
 806 {
 807         c.s << "basic_log_kernel()";
 808 }
 809
 810 GINAC_BIND_UNARCHIVER(basic_log_kernel);
 811
 812
 813 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(multiple_polylog_kernel, integration_kernel,
 814                                      print_func<print_context>(&multiple_polylog_kernel::do_print))
 815
 816 multiple_polylog_kernel::multiple_polylog_kernel() : inherited(), z(_ex1)
 817 {
 818 }
 819
 820 multiple_polylog_kernel::multiple_polylog_kernel(const ex & arg_z) : inherited(), z(arg_z)
 821 {
 822 }
 823
 824 int multiple_polylog_kernel::compare_same_type(const basic &other) const
 825 {
 826         const multiple_polylog_kernel &o = static_cast<const multiple_polylog_kernel &>(other);
 827
 828         return z.compare(o.z);
 829 }
 830
 831 size_t multiple_polylog_kernel::nops() const
 832 {
 833         return 1;
 834 }
 835
 836 ex multiple_polylog_kernel::op(size_t i) const
 837 {
 838         if ( i != 0 ) {
 839                 throw(std::range_error("multiple_polylog_kernel::op(): out of range"));
 840         }
 841
 842         return z;
 843 }
 844
 845 ex & multiple_polylog_kernel::let_op(size_t i)
 846 {
 847         ensure_if_modifiable();
 848
 849         if ( i != 0 ) {
 850                 throw(std::range_error("multiple_polylog_kernel::let_op(): out of range"));
 851         }
 852
 853         return z;
 854 }
 855
 856 bool multiple_polylog_kernel::is_numeric(void) const
 857 {
 858         return z.evalf().info(info_flags::numeric);
 859 }
 860
 861 cln::cl_N multiple_polylog_kernel::series_coeff_impl(int i) const
 862 {
 863         if ( i == 0 ) {
 864                 return 0;
 865         }
 866
 867         return -cln::expt(ex_to<numeric>(z.evalf()).to_cl_N(),-i);
 868 }
 869
 870 void multiple_polylog_kernel::do_print(const print_context & c, unsigned level) const
 871 {
 872         c.s << "multiple_polylog_kernel(";
 873         z.print(c);
 874         c.s << ")";
 875 }
 876
 877 GINAC_BIND_UNARCHIVER(multiple_polylog_kernel);
 878
 879
 880 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(ELi_kernel, integration_kernel,
 881                                      print_func<print_context>(&ELi_kernel::do_print))
 882
 883 ELi_kernel::ELi_kernel() : inherited(), n(_ex0), m(_ex0), x(_ex0), y(_ex0)
 884 {
 885 }
 886
 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)
 888 {
 889 }
 890
 891 int ELi_kernel::compare_same_type(const basic &other) const
 892 {
 893         const ELi_kernel &o = static_cast<const ELi_kernel &>(other);
 894         int cmpval;
 895
 896         cmpval = n.compare(o.n);
 897         if ( cmpval) {
 898                 return cmpval;
 899         }
 900
 901         cmpval = m.compare(o.m);
 902         if ( cmpval) {
 903                 return cmpval;
 904         }
 905
 906         cmpval = x.compare(o.x);
 907         if ( cmpval) {
 908                 return cmpval;
 909         }
 910
 911         return y.compare(o.y);
 912 }
 913
 914 size_t ELi_kernel::nops() const
 915 {
 916         return 4;
 917 }
 918
 919 ex ELi_kernel::op(size_t i) const
 920 {
 921         switch (i) {
 922         case 0:
 923                 return n;
 924         case 1:
 925                 return m;
 926         case 2:
 927                 return x;
 928         case 3:
 929                 return y;
 930         default:
 931                 throw (std::out_of_range("ELi_kernel::op() out of range"));
 932         }
 933 }
 934
 935 ex & ELi_kernel::let_op(size_t i)
 936 {
 937         ensure_if_modifiable();
 938
 939         switch (i) {
 940         case 0:
 941                 return n;
 942         case 1:
 943                 return m;
 944         case 2:
 945                 return x;
 946         case 3:
 947                 return y;
 948         default:
 949                 throw (std::out_of_range("ELi_kernel::let_op() out of range"));
 950         }
 951 }
 952
 953 bool ELi_kernel::is_numeric(void) const
 954 {
 955         return (n.info(info_flags::nonnegint) && m.info(info_flags::numeric) && x.evalf().info(info_flags::numeric) && y.evalf().info(info_flags::numeric));
 956 }
 957
 958 cln::cl_N ELi_kernel::series_coeff_impl(int i) const
 959 {
 960         if ( i == 0 ) {
 961                 return 0;
 962         }
 963
 964         int n_int = ex_to<numeric>(n).to_int();
 965         int m_int = ex_to<numeric>(m).to_int();
 966
 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();
 969
 970         cln::cl_N res_cln = 0;
 971
 972         for (int j=1; j<=i; j++) {
 973                 if ( (i % j) == 0 ) {
 974                         int k = i/j;
 975
 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);
 977                 }
 978         }
 979
 980         return res_cln;
 981 }
 982
 983 /**
 984  *
 985  * Returns the value of ELi_{n,m}(x,y,qbar)
 986  *
 987  */
 988 ex  ELi_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
 989 {
 990         return get_numerical_value_impl(qbar, 1, 1, N_trunc);
 991 }
 992
 993 void ELi_kernel::do_print(const print_context & c, unsigned level) const
 994 {
 995         c.s << "ELi_kernel(";
 996         n.print(c);
 997         c.s << ",";
 998         m.print(c);
 999         c.s << ",";
1000         x.print(c);
1001         c.s << ",";
1002         y.print(c);
1003         c.s << ")";
1004 }
1005
1006 GINAC_BIND_UNARCHIVER(ELi_kernel);
1007
1008
1009 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Ebar_kernel, integration_kernel,
1010                                      print_func<print_context>(&Ebar_kernel::do_print))
1011
1012 Ebar_kernel::Ebar_kernel() : inherited(), n(_ex0), m(_ex0), x(_ex0), y(_ex0)
1013 {
1014 }
1015
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)
1017 {
1018 }
1019
1020 int Ebar_kernel::compare_same_type(const basic &other) const
1021 {
1022         const Ebar_kernel &o = static_cast<const Ebar_kernel &>(other);
1023         int cmpval;
1024
1025         cmpval = n.compare(o.n);
1026         if ( cmpval) {
1027                 return cmpval;
1028         }
1029
1030         cmpval = m.compare(o.m);
1031         if ( cmpval) {
1032                 return cmpval;
1033         }
1034
1035         cmpval = x.compare(o.x);
1036         if ( cmpval) {
1037                 return cmpval;
1038         }
1039
1040         return y.compare(o.y);
1041 }
1042
1043 size_t Ebar_kernel::nops() const
1044 {
1045         return 4;
1046 }
1047
1048 ex Ebar_kernel::op(size_t i) const
1049 {
1050         switch (i) {
1051         case 0:
1052                 return n;
1053         case 1:
1054                 return m;
1055         case 2:
1056                 return x;
1057         case 3:
1058                 return y;
1059         default:
1060                 throw (std::out_of_range("Ebar_kernel::op() out of range"));
1061         }
1062 }
1063
1064 ex & Ebar_kernel::let_op(size_t i)
1065 {
1066         ensure_if_modifiable();
1067
1068         switch (i) {
1069         case 0:
1070                 return n;
1071         case 1:
1072                 return m;
1073         case 2:
1074                 return x;
1075         case 3:
1076                 return y;
1077         default:
1078                 throw (std::out_of_range("Ebar_kernel::let_op() out of range"));
1079         }
1080 }
1081
1082 bool Ebar_kernel::is_numeric(void) const
1083 {
1084         return (n.info(info_flags::nonnegint) && m.info(info_flags::numeric) && x.evalf().info(info_flags::numeric) && y.evalf().info(info_flags::numeric));
1085 }
1086
1087 cln::cl_N Ebar_kernel::series_coeff_impl(int i) const
1088 {
1089         if ( i == 0 ) {
1090                 return 0;
1091         }
1092
1093         int n_int = ex_to<numeric>(n).to_int();
1094         int m_int = ex_to<numeric>(m).to_int();
1095
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();
1098
1099         cln::cl_N res_cln = 0;
1100
1101         for (int j=1; j<=i; j++) {
1102                 if ( (i % j) == 0 ) {
1103                         int k = i/j;
1104
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);
1106                 }
1107         }
1108
1109         return res_cln;
1110 }
1111
1112 /**
1113  *
1114  * Returns the value of Ebar_{n,m}(x,y,qbar)
1115  *
1116  */
1117 ex  Ebar_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
1118 {
1119         return get_numerical_value_impl(qbar, 1, 1, N_trunc);
1120 }
1121
1122 void Ebar_kernel::do_print(const print_context & c, unsigned level) const
1123 {
1124         c.s << "Ebar_kernel(";
1125         n.print(c);
1126         c.s << ",";
1127         m.print(c);
1128         c.s << ",";
1129         x.print(c);
1130         c.s << ",";
1131         y.print(c);
1132         c.s << ")";
1133 }
1134
1135 GINAC_BIND_UNARCHIVER(Ebar_kernel);
1136
1137
1138 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Kronecker_dtau_kernel, integration_kernel,
1139                                      print_func<print_context>(&Kronecker_dtau_kernel::do_print))
1140
1141 Kronecker_dtau_kernel::Kronecker_dtau_kernel() : inherited(), n(_ex0), z(_ex0), K(_ex1), C_norm(_ex1)
1142 {
1143 }
1144
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)
1146 {
1147 }
1148
1149 int Kronecker_dtau_kernel::compare_same_type(const basic &other) const
1150 {
1151         const Kronecker_dtau_kernel &o = static_cast<const Kronecker_dtau_kernel &>(other);
1152         int cmpval;
1153
1154         cmpval = n.compare(o.n);
1155         if ( cmpval) {
1156                 return cmpval;
1157         }
1158
1159         cmpval = z.compare(o.z);
1160         if ( cmpval) {
1161                 return cmpval;
1162         }
1163
1164         cmpval = K.compare(o.K);
1165         if ( cmpval) {
1166                 return cmpval;
1167         }
1168
1169         return C_norm.compare(o.C_norm);
1170 }
1171
1172 size_t Kronecker_dtau_kernel::nops() const
1173 {
1174         return 4;
1175 }
1176
1177 ex Kronecker_dtau_kernel::op(size_t i) const
1178 {
1179         switch (i) {
1180         case 0:
1181                 return n;
1182         case 1:
1183                 return z;
1184         case 2:
1185                 return K;
1186         case 3:
1187                 return C_norm;
1188         default:
1189                 throw (std::out_of_range("Kronecker_dtau_kernel::op() out of range"));
1190         }
1191 }
1192
1193 ex & Kronecker_dtau_kernel::let_op(size_t i)
1194 {
1195         ensure_if_modifiable();
1196
1197         switch (i) {
1198         case 0:
1199                 return n;
1200         case 1:
1201                 return z;
1202         case 2:
1203                 return K;
1204         case 3:
1205                 return C_norm;
1206         default:
1207                 throw (std::out_of_range("Kronecker_dtau_kernel::let_op() out of range"));
1208         }
1209 }
1210
1211 bool Kronecker_dtau_kernel::is_numeric(void) const
1212 {
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));
1214 }
1215
1216 cln::cl_N Kronecker_dtau_kernel::series_coeff_impl(int i) const
1217 {
1218         numeric n_num = ex_to<numeric>(n);
1219         int n_int = n_num.to_int();
1220
1221         // case n=0
1222         if ( n_num == 0 ) {
1223                 if ( i == 0 ) {
1224                         ex res = -C_norm*K;
1225
1226                         return ex_to<numeric>(res.evalf()).to_cl_N();
1227                 }
1228
1229                 return 0;
1230         }
1231
1232         // case n=1
1233         if ( n_num == 1 ) {
1234                 return 0;
1235         }
1236
1237         // case n>1
1238         if ( i == 0 ) {
1239                 ex res = C_norm*K / factorial(n_num-2) * bernoulli(n_num)/n_num;
1240
1241                 return ex_to<numeric>(res.evalf()).to_cl_N();
1242         }
1243
1244         // n>1, i>0
1245
1246         // if K>1 the variable i needs to be a multiple of K
1247         int K_int = ex_to<numeric>(K).to_int();
1248
1249         if ( (i % K_int) != 0 ) {
1250                 return 0;
1251         }
1252         int i_local = i/K_int;
1253
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);
1260                 }
1261         }
1262         ex pre = -C_norm*K/factorial(n_num-2);
1263
1264         return ex_to<numeric>(pre.evalf()).to_cl_N() * res_cln;
1265 }
1266
1267 /**
1268  *
1269  * Returns the value of the g^(n)(z,K*tau), where tau is given by qbar.
1270  *
1271  */
1272 ex  Kronecker_dtau_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
1273 {
1274         numeric n_num = ex_to<numeric>(n);
1275
1276         if ( n_num == 0 ) {
1277                 return 1;
1278         }
1279
1280         // use the direct formula here
1281         if ( n_num == 1 ) {
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));
1284
1285                 return ex_to<numeric>(res.evalf());
1286         }
1287
1288         ex pre = pow(2*Pi*I,n_num)/C_norm/K/(n_num-1);
1289
1290         return get_numerical_value_impl(qbar, pre, 1, N_trunc);
1291 }
1292
1293 void Kronecker_dtau_kernel::do_print(const print_context & c, unsigned level) const
1294 {
1295         c.s << "Kronecker_dtau_kernel(";
1296         n.print(c);
1297         c.s << ",";
1298         z.print(c);
1299         c.s << ",";
1300         K.print(c);
1301         c.s << ",";
1302         C_norm.print(c);
1303         c.s << ")";
1304 }
1305
1306 GINAC_BIND_UNARCHIVER(Kronecker_dtau_kernel);
1307
1308
1309 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Kronecker_dz_kernel, integration_kernel,
1310                                      print_func<print_context>(&Kronecker_dz_kernel::do_print))
1311
1312 Kronecker_dz_kernel::Kronecker_dz_kernel() : inherited(), n(_ex0), z_j(_ex0), tau(_ex0), K(_ex1), C_norm(_ex1)
1313 {
1314 }
1315
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)
1317 {
1318 }
1319
1320 int Kronecker_dz_kernel::compare_same_type(const basic &other) const
1321 {
1322         const Kronecker_dz_kernel &o = static_cast<const Kronecker_dz_kernel &>(other);
1323         int cmpval;
1324
1325         cmpval = n.compare(o.n);
1326         if ( cmpval) {
1327                 return cmpval;
1328         }
1329
1330         cmpval = z_j.compare(o.z_j);
1331         if ( cmpval) {
1332                 return cmpval;
1333         }
1334
1335         cmpval = tau.compare(o.tau);
1336         if ( cmpval) {
1337                 return cmpval;
1338         }
1339
1340         cmpval = K.compare(o.K);
1341         if ( cmpval) {
1342                 return cmpval;
1343         }
1344
1345         return C_norm.compare(o.C_norm);
1346 }
1347
1348 size_t Kronecker_dz_kernel::nops() const
1349 {
1350         return 5;
1351 }
1352
1353 ex Kronecker_dz_kernel::op(size_t i) const
1354 {
1355         switch (i) {
1356         case 0:
1357                 return n;
1358         case 1:
1359                 return z_j;
1360         case 2:
1361                 return tau;
1362         case 3:
1363                 return K;
1364         case 4:
1365                 return C_norm;
1366         default:
1367                 throw (std::out_of_range("Kronecker_dz_kernel::op() out of range"));
1368         }
1369 }
1370
1371 ex & Kronecker_dz_kernel::let_op(size_t i)
1372 {
1373         ensure_if_modifiable();
1374
1375         switch (i) {
1376         case 0:
1377                 return n;
1378         case 1:
1379                 return z_j;
1380         case 2:
1381                 return tau;
1382         case 3:
1383                 return K;
1384         case 4:
1385                 return C_norm;
1386         default:
1387                 throw (std::out_of_range("Kronecker_dz_kernel::let_op() out of range"));
1388         }
1389 }
1390
1391 bool Kronecker_dz_kernel::is_numeric(void) const
1392 {
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));
1394 }
1395
1396 cln::cl_N Kronecker_dz_kernel::series_coeff_impl(int i) const
1397 {
1398         numeric n_num = ex_to<numeric>(n);
1399         int n_int = n_num.to_int();
1400
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();
1403
1404         ex qbar = exp(ex_to<numeric>((2*Pi*I*K*tau).evalf()));
1405
1406         // case n=0
1407         if ( n_num == 0 ) {
1408                 return 0;
1409         }
1410
1411         // case n=1
1412         if ( n_num == 1 ) {
1413                 if ( i == 1 ) {
1414                         return ex_to<numeric>((C_norm * 2*Pi*I).evalf()).to_cl_N();
1415                 }
1416
1417                 return 0;
1418         }
1419
1420         // case n=2
1421         if ( n_num == 2 ) {
1422                 if ( ex_to<numeric>(z_j.evalf()).is_zero() ) {
1423                         if ( i == 0 ) {
1424                                 return ex_to<numeric>((C_norm).evalf()).to_cl_N();
1425                         }
1426                         else if ( i == 1 ) {
1427                                 return 0;
1428                         }
1429                         else {
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);
1434                                 }
1435
1436                                 res *= -pow(2*Pi*I,i)*C_norm/factorial(i-1);
1437
1438                                 return ex_to<numeric>(res.evalf()).to_cl_N();
1439                         }
1440                 }
1441                 else {
1442                         // z_j is not zero
1443                         if ( i == 0 ) {
1444                                 return 0;
1445                         }
1446                         else {
1447                                 Li_negative my_Li_negative;
1448
1449                                 ex res = 0;
1450                                 if ( i == 1 ) {
1451                                         res = numeric(1,2);
1452                                 }
1453
1454                                 Ebar_kernel Ebar = Ebar_kernel( 1-i, 0, w_j_inv, numeric(1) );
1455
1456                                 res += my_Li_negative.get_numerical_value(i-1,w_j_inv) + Ebar.get_numerical_value(qbar);
1457
1458                                 res *= -pow(2*Pi*I,i)*C_norm/factorial(i-1);
1459
1460                                 return ex_to<numeric>(res.evalf()).to_cl_N();
1461                         }
1462                 }
1463         }
1464
1465         // case n>2
1466         ex res = 0;
1467         if ( i == 1 ) {
1468                 res += - bernoulli(n_num-1)/(n_num-1);
1469         }
1470         if ( i > 0 ) {
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) );
1474
1475                                 res += pow(2*Pi*I,i-1)/factorial(i-1) * Ebar.get_numerical_value(qbar);
1476                         }
1477                 }
1478                 else {
1479                         // z_j is not zero
1480                         Ebar_kernel Ebar = Ebar_kernel( 1-i, 2-n_num, w_j_inv, numeric(1) );
1481
1482                         res += pow(2*Pi*I,i-1)/factorial(i-1) * Ebar.get_numerical_value(qbar);
1483                 }
1484         }
1485
1486         res *= - C_norm * 2*Pi*I/factorial(n_num-2);
1487
1488         return ex_to<numeric>(res.evalf()).to_cl_N();
1489 }
1490
1491 /**
1492  *
1493  * Returns the value of the g^(n-1)(z-z_j,K*tau).
1494  *
1495  */
1496 ex  Kronecker_dz_kernel::get_numerical_value(const ex & z, int N_trunc) const
1497 {
1498         numeric n_num = ex_to<numeric>(n);
1499
1500         if ( n_num == 1 ) {
1501                 return 1;
1502         }
1503
1504         ex pre = pow(2*Pi*I,n-2)/C_norm;
1505
1506         return get_numerical_value_impl(z, pre, 0, N_trunc);
1507 }
1508
1509 void Kronecker_dz_kernel::do_print(const print_context & c, unsigned level) const
1510 {
1511         c.s << "Kronecker_dz_kernel(";
1512         n.print(c);
1513         c.s << ",";
1514         z_j.print(c);
1515         c.s << ",";
1516         tau.print(c);
1517         c.s << ",";
1518         K.print(c);
1519         c.s << ",";
1520         C_norm.print(c);
1521         c.s << ")";
1522 }
1523
1524 GINAC_BIND_UNARCHIVER(Kronecker_dz_kernel);
1525
1526
1527 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eisenstein_kernel, integration_kernel,
1528                                      print_func<print_context>(&Eisenstein_kernel::do_print))
1529
1530 Eisenstein_kernel::Eisenstein_kernel() : inherited(), k(_ex0), N(_ex0), a(_ex0), b(_ex0), K(_ex0), C_norm(_ex1)
1531 {
1532 }
1533
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)
1535 {
1536 }
1537
1538 int Eisenstein_kernel::compare_same_type(const basic &other) const
1539 {
1540         const Eisenstein_kernel &o = static_cast<const Eisenstein_kernel &>(other);
1541         int cmpval;
1542
1543         cmpval = k.compare(o.k);
1544         if ( cmpval) {
1545                 return cmpval;
1546         }
1547
1548         cmpval = N.compare(o.N);
1549         if ( cmpval) {
1550                 return cmpval;
1551         }
1552
1553         cmpval = a.compare(o.a);
1554         if ( cmpval) {
1555                 return cmpval;
1556         }
1557
1558         cmpval = b.compare(o.b);
1559         if ( cmpval) {
1560                 return cmpval;
1561         }
1562
1563         cmpval = K.compare(o.K);
1564         if ( cmpval) {
1565                 return cmpval;
1566         }
1567
1568         return C_norm.compare(o.C_norm);
1569 }
1570
1571 /**
1572  *
1573  * The series method for this class returns the qbar-expansion of the modular form,
1574  * without an additional factor of C_norm/qbar.
1575  *
1576  * This allows for easy use in the class modular_form_kernel.
1577  *
1578  */
1579 ex Eisenstein_kernel::series(const relational & r, int order, unsigned options) const
1580 {
1581         if ( r.rhs() != 0 ) {
1582                 throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
1583         }
1584
1585         ex qbar = r.lhs();
1586         ex res = q_expansion_modular_form(qbar, order);
1587         res = res.series(qbar,order);
1588
1589         return res;
1590 }
1591
1592 size_t Eisenstein_kernel::nops() const
1593 {
1594         return 6;
1595 }
1596
1597 ex Eisenstein_kernel::op(size_t i) const
1598 {
1599         switch (i) {
1600         case 0:
1601                 return k;
1602         case 1:
1603                 return N;
1604         case 2:
1605                 return a;
1606         case 3:
1607                 return b;
1608         case 4:
1609                 return K;
1610         case 5:
1611                 return C_norm;
1612         default:
1613                 throw (std::out_of_range("Eisenstein_kernel::op() out of range"));
1614         }
1615 }
1616
1617 ex & Eisenstein_kernel::let_op(size_t i)
1618 {
1619         ensure_if_modifiable();
1620
1621         switch (i) {
1622         case 0:
1623                 return k;
1624         case 1:
1625                 return N;
1626         case 2:
1627                 return a;
1628         case 3:
1629                 return b;
1630         case 4:
1631                 return K;
1632         case 5:
1633                 return C_norm;
1634         default:
1635                 throw (std::out_of_range("Eisenstein_kernel::let_op() out of range"));
1636         }
1637 }
1638
1639 bool Eisenstein_kernel::is_numeric(void) const
1640 {
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));
1642 }
1643
1644 ex Eisenstein_kernel::Laurent_series(const ex & x, int order) const
1645 {
1646         ex res = C_norm * q_expansion_modular_form(x, order+1)/x;
1647         res = res.series(x,order);
1648
1649         return res;
1650 }
1651
1652 /**
1653  *
1654  * Returns the value of the modular form.
1655  *
1656  */
1657 ex  Eisenstein_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
1658 {
1659         ex pre = numeric(1)/C_norm;
1660
1661         return get_numerical_value_impl(qbar, pre, 1, N_trunc);
1662 }
1663
1664 bool Eisenstein_kernel::uses_Laurent_series() const
1665 {
1666         return true;
1667 }
1668
1669 ex Eisenstein_kernel::q_expansion_modular_form(const ex & q, int order) const
1670 {
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);
1676
1677         if ( (k==2) && (a==1) && (b==1) ) {
1678                 return B_eisenstein_series(q, N_num, K_num, order);
1679         }
1680
1681         return E_eisenstein_series(q, k_num, N_num, a_num, b_num, K_num, order);
1682 }
1683
1684 void Eisenstein_kernel::do_print(const print_context & c, unsigned level) const
1685 {
1686         c.s << "Eisenstein_kernel(";
1687         k.print(c);
1688         c.s << ",";
1689         N.print(c);
1690         c.s << ",";
1691         a.print(c);
1692         c.s << ",";
1693         b.print(c);
1694         c.s << ",";
1695         K.print(c);
1696         c.s << ",";
1697         C_norm.print(c);
1698         c.s << ")";
1699 }
1700
1701 GINAC_BIND_UNARCHIVER(Eisenstein_kernel);
1702
1703
1704 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eisenstein_h_kernel, integration_kernel,
1705                                      print_func<print_context>(&Eisenstein_h_kernel::do_print))
1706
1707 Eisenstein_h_kernel::Eisenstein_h_kernel() : inherited(), k(_ex0), N(_ex0), r(_ex0), s(_ex0), C_norm(_ex1)
1708 {
1709 }
1710
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)
1712 {
1713 }
1714
1715 int Eisenstein_h_kernel::compare_same_type(const basic &other) const
1716 {
1717         const Eisenstein_h_kernel &o = static_cast<const Eisenstein_h_kernel &>(other);
1718         int cmpval;
1719
1720         cmpval = k.compare(o.k);
1721         if ( cmpval) {
1722                 return cmpval;
1723         }
1724
1725         cmpval = N.compare(o.N);
1726         if ( cmpval) {
1727                 return cmpval;
1728         }
1729
1730         cmpval = r.compare(o.r);
1731         if ( cmpval) {
1732                 return cmpval;
1733         }
1734
1735         cmpval = s.compare(o.s);
1736         if ( cmpval) {
1737                 return cmpval;
1738         }
1739
1740         return C_norm.compare(o.C_norm);
1741 }
1742
1743 /**
1744  *
1745  * The series method for this class returns the qbar-expansion of the modular form,
1746  * without an additional factor of C_norm/qbar.
1747  *
1748  * This allows for easy use in the class modular_form_kernel.
1749  *
1750  */
1751 ex Eisenstein_h_kernel::series(const relational & r, int order, unsigned options) const
1752 {
1753         if ( r.rhs() != 0 ) {
1754                 throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
1755         }
1756
1757         ex qbar = r.lhs();
1758         ex res = q_expansion_modular_form(qbar, order);
1759         res = res.series(qbar,order);
1760
1761         return res;
1762 }
1763
1764 size_t Eisenstein_h_kernel::nops() const
1765 {
1766         return 5;
1767 }
1768
1769 ex Eisenstein_h_kernel::op(size_t i) const
1770 {
1771         switch (i) {
1772         case 0:
1773                 return k;
1774         case 1:
1775                 return N;
1776         case 2:
1777                 return r;
1778         case 3:
1779                 return s;
1780         case 4:
1781                 return C_norm;
1782         default:
1783                 throw (std::out_of_range("Eisenstein_h_kernel::op() out of range"));
1784         }
1785 }
1786
1787 ex & Eisenstein_h_kernel::let_op(size_t i)
1788 {
1789         ensure_if_modifiable();
1790
1791         switch (i) {
1792         case 0:
1793                 return k;
1794         case 1:
1795                 return N;
1796         case 2:
1797                 return r;
1798         case 3:
1799                 return s;
1800         case 4:
1801                 return C_norm;
1802         default:
1803                 throw (std::out_of_range("Eisenstein_h_kernel::let_op() out of range"));
1804         }
1805 }
1806
1807 bool Eisenstein_h_kernel::is_numeric(void) const
1808 {
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));
1810 }
1811
1812 ex Eisenstein_h_kernel::Laurent_series(const ex & x, int order) const
1813 {
1814         ex res = C_norm * q_expansion_modular_form(x, order+1)/x;
1815         res = res.series(x,order);
1816
1817         return res;
1818 }
1819
1820 /**
1821  *
1822  * Returns the value of the modular form.
1823  *
1824  */
1825 ex  Eisenstein_h_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
1826 {
1827         ex pre = numeric(1)/C_norm;
1828
1829         return get_numerical_value_impl(qbar, pre, 1, N_trunc);
1830 }
1831
1832 bool Eisenstein_h_kernel::uses_Laurent_series() const
1833 {
1834         return true;
1835 }
1836
1837 /**
1838  *
1839  * The constant coefficient in the Fourier expansion.
1840  *
1841  */
1842 ex Eisenstein_h_kernel::coefficient_a0(const numeric & k, const numeric & r, const numeric & s, const numeric & N) const
1843 {
1844         if ( k == 1 ) {
1845                 if ( irem(s,N) != 0 ) {
1846                         return numeric(1,4) - mod(s,N)/numeric(2)/N;
1847                 }
1848                 else if ( (irem(r,N)==0) && (irem(s,N)==0) ) {
1849                         return 0;
1850                 }
1851                 else {
1852                         return I*numeric(1,4)*cos(Pi*mod(r,N)/N)/sin(Pi*mod(r,N)/N);
1853                 }
1854         }
1855
1856         // case k > 1
1857         return -Bernoulli_polynomial(k,mod(s,N)/N)/numeric(2)/k;
1858 }
1859
1860 /**
1861  *
1862  * The higher coefficients in the Fourier expansion.
1863  *
1864  */
1865 ex Eisenstein_h_kernel::coefficient_an(const numeric & n, const numeric & k, const numeric & r, const numeric & s, const numeric & N) const
1866 {
1867         ex res = 0;
1868
1869         for (numeric m=1; m<=n; m++) {
1870                 if ( irem(n,m) == 0 ) {
1871                         for (numeric c1=0; c1<N; c1++) {
1872                                 numeric c2 = n/m;
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));
1874                         }
1875                 }
1876         }
1877
1878         return res/numeric(2)/pow(N,k);
1879 }
1880
1881 ex Eisenstein_h_kernel::q_expansion_modular_form(const ex & q, int N_order) const
1882 {
1883         numeric N_order_num = numeric(N_order);
1884
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);
1889
1890         ex res = coefficient_a0(k_num,r_num,s_num,N_num);
1891
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);
1894         }
1895
1896         res += Order(pow(q,N_order));
1897         res = res.series(q,N_order);
1898
1899         return res;
1900 }
1901
1902 void Eisenstein_h_kernel::do_print(const print_context & c, unsigned level) const
1903 {
1904         c.s << "Eisenstein_h_kernel(";
1905         k.print(c);
1906         c.s << ",";
1907         N.print(c);
1908         c.s << ",";
1909         r.print(c);
1910         c.s << ",";
1911         s.print(c);
1912         c.s << ",";
1913         C_norm.print(c);
1914         c.s << ")";
1915 }
1916
1917 GINAC_BIND_UNARCHIVER(Eisenstein_h_kernel);
1918
1919
1920 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(modular_form_kernel, integration_kernel,
1921                                      print_func<print_context>(&modular_form_kernel::do_print))
1922
1923 modular_form_kernel::modular_form_kernel() : inherited(), k(_ex0), P(_ex0), C_norm(_ex1)
1924 {
1925 }
1926
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)
1928 {
1929 }
1930
1931 int modular_form_kernel::compare_same_type(const basic &other) const
1932 {
1933         const modular_form_kernel &o = static_cast<const modular_form_kernel &>(other);
1934         int cmpval;
1935
1936         cmpval = k.compare(o.k);
1937         if ( cmpval) {
1938                 return cmpval;
1939         }
1940
1941         cmpval = P.compare(o.P);
1942         if ( cmpval) {
1943                 return cmpval;
1944         }
1945
1946         return C_norm.compare(o.C_norm);
1947 }
1948
1949 /**
1950  *
1951  * The series method for this class returns the qbar-expansion of the modular form,
1952  * without an additional factor of C_norm/qbar.
1953  *
1954  */
1955 ex modular_form_kernel::series(const relational & r, int order, unsigned options) const
1956 {
1957         if ( r.rhs() != 0 ) {
1958                 throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
1959         }
1960
1961         ex qbar = r.lhs();
1962
1963         subs_q_expansion do_subs_q_expansion(qbar, order);
1964
1965         ex res = do_subs_q_expansion(P).series(qbar,order);
1966
1967         return res;
1968 }
1969
1970 size_t modular_form_kernel::nops() const
1971 {
1972         return 3;
1973 }
1974
1975 ex modular_form_kernel::op(size_t i) const
1976 {
1977         switch (i) {
1978         case 0:
1979                 return k;
1980         case 1:
1981                 return P;
1982         case 2:
1983                 return C_norm;
1984         default:
1985                 throw (std::out_of_range("modular_form_kernel::op() out of range"));
1986         }
1987 }
1988
1989 ex & modular_form_kernel::let_op(size_t i)
1990 {
1991         ensure_if_modifiable();
1992
1993         switch (i) {
1994         case 0:
1995                 return k;
1996         case 1:
1997                 return P;
1998         case 2:
1999                 return C_norm;
2000         default:
2001                 throw (std::out_of_range("modular_form_kernel::let_op() out of range"));
2002         }
2003 }
2004
2005 bool modular_form_kernel::is_numeric(void) const
2006 {
2007         bool flag = (k.info(info_flags::nonnegint) && C_norm.evalf().info(info_flags::numeric));
2008         if ( !flag ) {
2009                 return false;
2010         }
2011
2012         symbol qbar("qbar");
2013
2014         // test with a random number and random expansion
2015         return series_to_poly(P.series(qbar,18)).subs(qbar==numeric(1,937)).evalf().info(info_flags::numeric);
2016 }
2017
2018 ex modular_form_kernel::Laurent_series(const ex & qbar, int order) const
2019 {
2020         ex res = series_to_poly(P.series(qbar,order+1));
2021         res = C_norm * res/qbar;
2022         res = res.series(qbar,order);
2023         return res;
2024 }
2025
2026 /**
2027  *
2028  * Returns the value of the modular form.
2029  *
2030  */
2031 ex  modular_form_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
2032 {
2033         ex pre = numeric(1)/C_norm;
2034
2035         return get_numerical_value_impl(qbar, pre, 1, N_trunc);
2036 }
2037
2038 bool modular_form_kernel::uses_Laurent_series() const
2039 {
2040         return true;
2041 }
2042
2043 ex modular_form_kernel::q_expansion_modular_form(const ex & q, int N_order) const
2044 {
2045         return this->series(q==0,N_order);
2046 }
2047
2048 void modular_form_kernel::do_print(const print_context & c, unsigned level) const
2049 {
2050         c.s << "modular_form_kernel(";
2051         k.print(c);
2052         c.s << ",";
2053         P.print(c);
2054         c.s << ",";
2055         C_norm.print(c);
2056         c.s << ")";
2057 }
2058
2059 GINAC_BIND_UNARCHIVER(modular_form_kernel);
2060
2061
2062 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(user_defined_kernel, integration_kernel,
2063                                      print_func<print_context>(&user_defined_kernel::do_print))
2064
2065 user_defined_kernel::user_defined_kernel() : inherited(), f(_ex0), x(_ex0)
2066 {
2067 }
2068
2069 user_defined_kernel::user_defined_kernel(const ex & arg_f, const ex & arg_x) : inherited(), f(arg_f), x(arg_x)
2070 {
2071 }
2072
2073 int user_defined_kernel::compare_same_type(const basic &other) const
2074 {
2075         const user_defined_kernel &o = static_cast<const user_defined_kernel &>(other);
2076         int cmpval;
2077
2078         cmpval = f.compare(o.f);
2079         if ( cmpval) {
2080                 return cmpval;
2081         }
2082
2083         return x.compare(o.x);
2084 }
2085
2086 size_t user_defined_kernel::nops() const
2087 {
2088         return 2;
2089 }
2090
2091 ex user_defined_kernel::op(size_t i) const
2092 {
2093         switch (i) {
2094         case 0:
2095                 return f;
2096         case 1:
2097                 return x;
2098         default:
2099                 throw (std::out_of_range("user_defined_kernel::op() out of range"));
2100         }
2101 }
2102
2103 ex & user_defined_kernel::let_op(size_t i)
2104 {
2105         ensure_if_modifiable();
2106
2107         switch (i) {
2108         case 0:
2109                 return f;
2110         case 1:
2111                 return x;
2112         default:
2113                 throw (std::out_of_range("user_defined_kernel::let_op() out of range"));
2114         }
2115 }
2116
2117 bool user_defined_kernel::is_numeric(void) const
2118 {
2119         // test with a random number
2120         return f.subs(x==numeric(1,937)).evalf().info(info_flags::numeric);
2121 }
2122
2123 ex user_defined_kernel::Laurent_series(const ex & x_up, int order) const
2124 {
2125         ex res = f.series(x,order).subs(x==x_up);
2126
2127         return res;
2128 }
2129
2130 bool user_defined_kernel::uses_Laurent_series() const
2131 {
2132         return true;
2133 }
2134
2135 void user_defined_kernel::do_print(const print_context & c, unsigned level) const
2136 {
2137         c.s << "user_defined_kernel(";
2138         f.print(c);
2139         c.s << ",";
2140         x.print(c);
2141         c.s << ")";
2142 }
2143
2144 GINAC_BIND_UNARCHIVER(user_defined_kernel);
2145
2146 } // namespace GiNaC