Minor bug fix for the class modular_form_kernel: Ensure that the series-method includ...
[ginac.git] / 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-2021 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         res += Order(pow(qbar,order));
1967         res = res.series(qbar,order);
1968
1969         return res;
1970 }
1971
1972 size_t modular_form_kernel::nops() const
1973 {
1974         return 3;
1975 }
1976
1977 ex modular_form_kernel::op(size_t i) const
1978 {
1979         switch (i) {
1980         case 0:
1981                 return k;
1982         case 1:
1983                 return P;
1984         case 2:
1985                 return C_norm;
1986         default:
1987                 throw (std::out_of_range("modular_form_kernel::op() out of range"));
1988         }
1989 }
1990
1991 ex & modular_form_kernel::let_op(size_t i)
1992 {
1993         ensure_if_modifiable();
1994
1995         switch (i) {
1996         case 0:
1997                 return k;
1998         case 1:
1999                 return P;
2000         case 2:
2001                 return C_norm;
2002         default:
2003                 throw (std::out_of_range("modular_form_kernel::let_op() out of range"));
2004         }
2005 }
2006
2007 bool modular_form_kernel::is_numeric(void) const
2008 {
2009         bool flag = (k.info(info_flags::nonnegint) && C_norm.evalf().info(info_flags::numeric));
2010         if ( !flag ) {
2011                 return false;
2012         }
2013
2014         symbol qbar("qbar");
2015
2016         // test with a random number and random expansion
2017         return series_to_poly(q_expansion_modular_form(qbar,18)).subs(qbar==numeric(1,937)).evalf().info(info_flags::numeric);
2018 }
2019
2020 ex modular_form_kernel::Laurent_series(const ex & qbar, int order) const
2021 {
2022         ex res = series_to_poly(q_expansion_modular_form(qbar,order+1));
2023         res = C_norm * res/qbar;
2024         res = res.series(qbar,order);
2025         return res;
2026 }
2027
2028 /**
2029  *
2030  * Returns the value of the modular form.
2031  *
2032  */
2033 ex  modular_form_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
2034 {
2035         ex pre = numeric(1)/C_norm;
2036
2037         return get_numerical_value_impl(qbar, pre, 1, N_trunc);
2038 }
2039
2040 bool modular_form_kernel::uses_Laurent_series() const
2041 {
2042         return true;
2043 }
2044
2045 ex modular_form_kernel::q_expansion_modular_form(const ex & q, int N_order) const
2046 {
2047         return this->series(q==0,N_order);
2048 }
2049
2050 void modular_form_kernel::do_print(const print_context & c, unsigned level) const
2051 {
2052         c.s << "modular_form_kernel(";
2053         k.print(c);
2054         c.s << ",";
2055         P.print(c);
2056         c.s << ",";
2057         C_norm.print(c);
2058         c.s << ")";
2059 }
2060
2061 GINAC_BIND_UNARCHIVER(modular_form_kernel);
2062
2063
2064 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(user_defined_kernel, integration_kernel,
2065                                      print_func<print_context>(&user_defined_kernel::do_print))
2066
2067 user_defined_kernel::user_defined_kernel() : inherited(), f(_ex0), x(_ex0)
2068
2069 }
2070
2071 user_defined_kernel::user_defined_kernel(const ex & arg_f, const ex & arg_x) : inherited(), f(arg_f), x(arg_x)
2072
2073 }
2074
2075 int user_defined_kernel::compare_same_type(const basic &other) const
2076 {
2077         const user_defined_kernel &o = static_cast<const user_defined_kernel &>(other);
2078         int cmpval;
2079
2080         cmpval = f.compare(o.f);
2081         if ( cmpval) {
2082                 return cmpval;
2083         }
2084
2085         return x.compare(o.x);
2086 }
2087
2088 size_t user_defined_kernel::nops() const
2089 {
2090         return 2;
2091 }
2092
2093 ex user_defined_kernel::op(size_t i) const
2094 {
2095         switch (i) {
2096         case 0:
2097                 return f;
2098         case 1:
2099                 return x;
2100         default:
2101                 throw (std::out_of_range("user_defined_kernel::op() out of range"));
2102         }
2103 }
2104
2105 ex & user_defined_kernel::let_op(size_t i)
2106 {
2107         ensure_if_modifiable();
2108
2109         switch (i) {
2110         case 0:
2111                 return f;
2112         case 1:
2113                 return x;
2114         default:
2115                 throw (std::out_of_range("user_defined_kernel::let_op() out of range"));
2116         }
2117 }
2118
2119 bool user_defined_kernel::is_numeric(void) const
2120 {
2121         // test with a random number
2122         return f.subs(x==numeric(1,937)).evalf().info(info_flags::numeric);
2123 }
2124
2125 ex user_defined_kernel::Laurent_series(const ex & x_up, int order) const
2126 {
2127         ex res = f.series(x,order).subs(x==x_up);
2128
2129         return res;
2130 }
2131
2132 bool user_defined_kernel::uses_Laurent_series() const
2133 {
2134         return true;
2135 }
2136
2137 void user_defined_kernel::do_print(const print_context & c, unsigned level) const
2138 {
2139         c.s << "user_defined_kernel(";
2140         f.print(c);
2141         c.s << ",";
2142         x.print(c);
2143         c.s << ")";
2144 }
2145
2146 GINAC_BIND_UNARCHIVER(user_defined_kernel);
2147
2148 } // namespace GiNaC