GiNaC  1.8.3
integration_kernel.cpp
Go to the documentation of this file.
1 
5 /*
6  * GiNaC Copyright (C) 1999-2022 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 
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 
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 
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 
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 
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 
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 
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 
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 
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 
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  }
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 
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 
479 {
480  return kronecker_symbol(a,n);
481 }
482 
495 numeric dirichlet_character(const numeric & n, const numeric & a, const numeric & N)
496 {
497  if ( gcd(n,N) == 1 ) {
499  }
500 
501  return 0;
502 }
503 
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 
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 
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 
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 
599 bool integration_kernel::is_numeric(void) const
600 {
601  return true;
602 }
603 
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 
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 
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 
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 
695 bool integration_kernel::uses_Laurent_series() const
696 {
697  return false;
698 }
699 
705 size_t integration_kernel::get_cache_size(void) const
706 {
707  return series_vec.size();
708 }
709 
715 void integration_kernel::set_cache_step(int cache_steps) const
716 {
717  cache_step_size = cache_steps;
718 }
719 
725 ex integration_kernel::get_series_coeff(int i) const
726 {
727  return numeric(series_coeff(i));
728 }
729 
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 
782 
783 
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 
811 
812 
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 
825 {
826  const multiple_polylog_kernel &o = static_cast<const multiple_polylog_kernel &>(other);
827 
828  return z.compare(o.z);
829 }
830 
832 {
833  return 1;
834 }
835 
837 {
838  if ( i != 0 ) {
839  throw(std::range_error("multiple_polylog_kernel::op(): out of range"));
840  }
841 
842  return z;
843 }
844 
846 {
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 
857 {
858  return z.evalf().info(info_flags::numeric);
859 }
860 
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 
878 
879 
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 
936 {
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 {
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 
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 
1007 
1008 
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 
1065 {
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 {
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 
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 
1136 
1137 
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 
1173 {
1174  return 4;
1175 }
1176 
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 
1194 {
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 
1212 {
1214 }
1215 
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 
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 
1307 
1308 
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 
1349 {
1350  return 5;
1351 }
1352 
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 
1372 {
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 
1392 {
1394 }
1395 
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 
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 
1525 
1526 
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 
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();
1587  res = res.series(qbar,order);
1588 
1589  return res;
1590 }
1591 
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 
1618 {
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 
1640 {
1642 }
1643 
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 
1658 {
1659  ex pre = numeric(1)/C_norm;
1660 
1661  return get_numerical_value_impl(qbar, pre, 1, N_trunc);
1662 }
1663 
1665 {
1666  return true;
1667 }
1668 
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 
1702 
1703 
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 
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();
1759  res = res.series(qbar,order);
1760 
1761  return res;
1762 }
1763 
1765 {
1766  return 5;
1767 }
1768 
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 
1788 {
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 
1808 {
1810 }
1811 
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 
1826 {
1827  ex pre = numeric(1)/C_norm;
1828 
1829  return get_numerical_value_impl(qbar, pre, 1, N_trunc);
1830 }
1831 
1833 {
1834  return true;
1835 }
1836 
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 
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 
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 
1918 
1919 
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 
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 
1973 {
1974  return 3;
1975 }
1976 
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 
1992 {
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 
2008 {
2010  if ( !flag ) {
2011  return false;
2012  }
2013 
2014  symbol qbar("qbar");
2015 
2016  // test with a random number and random expansion
2018 }
2019 
2021 {
2023  res = C_norm * res/qbar;
2024  res = res.series(qbar,order);
2025  return res;
2026 }
2027 
2034 {
2035  ex pre = numeric(1)/C_norm;
2036 
2037  return get_numerical_value_impl(qbar, pre, 1, N_trunc);
2038 }
2039 
2041 {
2042  return true;
2043 }
2044 
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 
2062 
2063 
2066 
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 
2089 {
2090  return 2;
2091 }
2092 
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 
2106 {
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 
2120 {
2121  // test with a random number
2122  return f.subs(x==numeric(1,937)).evalf().info(info_flags::numeric);
2123 }
2124 
2126 {
2127  ex res = f.series(x,order).subs(x==x_up);
2128 
2129  return res;
2130 }
2131 
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 
2147 
2148 } // namespace GiNaC
Interface to GiNaC's sums of expressions.
#define GINAC_BIND_UNARCHIVER(classname)
Definition: archive.h:230
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
void do_print(const print_context &c, unsigned level) const
bool is_numeric(void) const override
This routine returns true, if the integration kernel can be evaluated numerically.
ELi_kernel(const ex &n, const ex &m, const ex &x, const ex &y)
size_t nops() const override
Number of operands/members.
ex op(size_t i) const override
Return operand/member at position i.
cln::cl_N series_coeff_impl(int i) const override
For only the coefficient of is non-zero.
ex get_numerical_value(const ex &qbar, int N_trunc=0) const override
Returns the value of ELi_{n,m}(x,y,qbar)
The Ebar-kernel.
size_t nops() const override
Number of operands/members.
void do_print(const print_context &c, unsigned level) const
bool is_numeric(void) const override
This routine returns true, if the integration kernel can be evaluated numerically.
ex op(size_t i) const override
Return operand/member at position i.
Ebar_kernel(const ex &n, const ex &m, const ex &x, const ex &y)
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
cln::cl_N series_coeff_impl(int i) const override
For only the coefficient of is non-zero.
ex get_numerical_value(const ex &qbar, int N_trunc=0) const override
Returns the value of Ebar_{n,m}(x,y,qbar)
The kernel corresponding to the Eisenstein series .
ex get_numerical_value(const ex &qbar, int N_trunc=0) const override
Returns the value of the modular form.
bool is_numeric(void) const override
This routine returns true, if the integration kernel can be evaluated numerically.
Eisenstein_h_kernel(const ex &k, const ex &N, const ex &r, const ex &s, const ex &C_norm=numeric(1))
size_t nops() const override
Number of operands/members.
ex op(size_t i) const override
Return operand/member at position i.
void do_print(const print_context &c, unsigned level) const
ex q_expansion_modular_form(const ex &q, int order) const
ex coefficient_an(const numeric &n, const numeric &k, const numeric &r, const numeric &s, const numeric &N) const
The higher coefficients in the Fourier expansion.
ex coefficient_a0(const numeric &k, const numeric &r, const numeric &s, const numeric &N) const
The constant coefficient in the Fourier expansion.
ex Laurent_series(const ex &x, int order) const override
Returns the Laurent series, starting possibly with the pole term.
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
bool uses_Laurent_series() const override
Returns true, if the coefficients are computed from the Laurent series (in which case the method Laur...
ex series(const relational &r, int order, unsigned options=0) const override
The series method for this class returns the qbar-expansion of the modular form, without an additiona...
The kernel corresponding to the Eisenstein series .
ex get_numerical_value(const ex &qbar, int N_trunc=0) const override
Returns the value of the modular form.
ex Laurent_series(const ex &x, int order) const override
Returns the Laurent series, starting possibly with the pole term.
void do_print(const print_context &c, unsigned level) const
ex op(size_t i) const override
Return operand/member at position i.
ex q_expansion_modular_form(const ex &q, int order) const
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
Eisenstein_kernel(const ex &k, const ex &N, const ex &a, const ex &b, const ex &K, const ex &C_norm=numeric(1))
size_t nops() const override
Number of operands/members.
bool uses_Laurent_series() const override
Returns true, if the coefficients are computed from the Laurent series (in which case the method Laur...
ex series(const relational &r, int order, unsigned options=0) const override
The series method for this class returns the qbar-expansion of the modular form, without an additiona...
bool is_numeric(void) const override
This routine returns true, if the integration kernel can be evaluated numerically.
The kernel corresponding to integrating the Kronecker coefficient function in (or equivalently in )...
size_t nops() const override
Number of operands/members.
Kronecker_dtau_kernel(const ex &n, const ex &z, const ex &K=numeric(1), const ex &C_norm=numeric(1))
ex op(size_t i) const override
Return operand/member at position i.
ex get_numerical_value(const ex &qbar, int N_trunc=0) const override
Returns the value of the g^(n)(z,K*tau), where tau is given by qbar.
void do_print(const print_context &c, unsigned level) const
cln::cl_N series_coeff_impl(int i) const override
For only the coefficient of is non-zero.
bool is_numeric(void) const override
This routine returns true, if the integration kernel can be evaluated numerically.
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
The kernel corresponding to integrating the Kronecker coefficient function in .
Kronecker_dz_kernel(const ex &n, const ex &z_j, const ex &tau, const ex &K=numeric(1), const ex &C_norm=numeric(1))
bool is_numeric(void) const override
This routine returns true, if the integration kernel can be evaluated numerically.
size_t nops() const override
Number of operands/members.
void do_print(const print_context &c, unsigned level) const
cln::cl_N series_coeff_impl(int i) const override
For only the coefficient of is non-zero.
ex op(size_t i) const override
Return operand/member at position i.
ex get_numerical_value(const ex &z, int N_trunc=0) const override
Returns the value of the g^(n-1)(z-z_j,K*tau).
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
The basic integration kernel with a logarithmic singularity at the origin.
This class is the ABC (abstract base class) of GiNaC's class hierarchy.
Definition: basic.h:105
void ensure_if_modifiable() const
Ensure the object may be modified without hurting others, throws if this is not the case.
Definition: basic.cpp:894
virtual int compare_same_type(const basic &other) const
Returns order relation between two objects of same type.
Definition: basic.cpp:719
virtual ex evalf() const
Evaluate object numerically.
Definition: basic.cpp:425
Wrapper template for making GiNaC classes out of STL containers.
Definition: container.h:73
const_iterator end() const
Definition: container.h:240
const_iterator begin() const
Definition: container.h:239
size_t nops() const override
Number of operands/members.
Definition: container.h:118
ex op(size_t i) const override
Return operand/member at position i.
Definition: container.h:295
container & append(const ex &b)
Add element at back.
Definition: container.h:391
Lightweight wrapper for GiNaC's symbolic objects.
Definition: ex.h:72
ex evalf() const
Definition: ex.h:121
ex series(const ex &r, int order, unsigned options=0) const
Compute the truncated series expansion of an expression.
Definition: pseries.cpp:1259
ex subs(const exmap &m, unsigned options=0) const
Definition: ex.h:826
bool info(unsigned inf) const
Definition: ex.h:132
int compare(const ex &other) const
Definition: ex.h:322
ex lhs() const
Left hand side of relational expression.
Definition: ex.cpp:225
void print(const print_context &c, unsigned level=0) const
Print expression to stream.
Definition: ex.cpp:56
ex rhs() const
Right hand side of relational expression.
Definition: ex.cpp:233
ex coeff(const ex &s, int n=1) const
Definition: ex.h:175
The base class for integration kernels for iterated integrals.
void do_print(const print_context &c, unsigned level) const
ex get_numerical_value_impl(const ex &lambda, const ex &pre, int shift, int N_trunc) const
The actual implementation for computing a numerical value for the integrand.
A kernel corresponding to a polynomial in Eisenstein series.
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
ex series(const relational &r, int order, unsigned options=0) const override
The series method for this class returns the qbar-expansion of the modular form, without an additiona...
void do_print(const print_context &c, unsigned level) const
size_t nops() const override
Number of operands/members.
bool is_numeric(void) const override
This routine returns true, if the integration kernel can be evaluated numerically.
modular_form_kernel(const ex &k, const ex &P, const ex &C_norm=numeric(1))
ex op(size_t i) const override
Return operand/member at position i.
ex get_numerical_value(const ex &qbar, int N_trunc=0) const override
Returns the value of the modular form.
ex Laurent_series(const ex &qbar, int order) const override
Returns the Laurent series, starting possibly with the pole term.
bool uses_Laurent_series() const override
Returns true, if the coefficients are computed from the Laurent series (in which case the method Laur...
ex q_expansion_modular_form(const ex &q, int order) const
The integration kernel for multiple polylogarithms.
cln::cl_N series_coeff_impl(int i) const override
For only the coefficient of is non-zero.
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
size_t nops() const override
Number of operands/members.
void do_print(const print_context &c, unsigned level) const
ex op(size_t i) const override
Return operand/member at position i.
bool is_numeric(void) const override
This routine returns true, if the integration kernel can be evaluated numerically.
This class is a wrapper around CLN-numbers within the GiNaC class hierarchy.
Definition: numeric.h:82
bool is_odd() const
True if object is an exact odd integer.
Definition: numeric.cpp:1182
cln::cl_N to_cl_N() const
Returns a new CLN object of type cl_N, representing the value of *this.
Definition: numeric.cpp:1332
bool is_even() const
True if object is an exact even integer.
Definition: numeric.cpp:1175
int to_int() const
Converts numeric types to machine's int.
Definition: numeric.cpp:1303
Base class for print_contexts.
Definition: print.h:103
This class holds a relation consisting of two expressions and a logical relation between them.
Definition: relational.h:35
Basic CAS symbol.
Definition: symbol.h:39
A user-defined integration kernel.
user_defined_kernel(const ex &f, const ex &x)
ex Laurent_series(const ex &x, int order) const override
Returns the Laurent series, starting possibly with the pole term.
bool is_numeric(void) const override
This routine returns true, if the integration kernel can be evaluated numerically.
ex op(size_t i) const override
Return operand/member at position i.
bool uses_Laurent_series() const override
Returns true, if the coefficients are computed from the Laurent series (in which case the method Laur...
void do_print(const print_context &c, unsigned level) const
ex & let_op(size_t i) override
Return modifiable operand/member at position i.
size_t nops() const override
Number of operands/members.
Interface to GiNaC's constant types and some special constants.
unsigned options
Definition: factor.cpp:2480
vector< int > k
Definition: factor.cpp:1466
size_t n
Definition: factor.cpp:1463
size_t c
Definition: factor.cpp:770
size_t r
Definition: factor.cpp:770
umodpoly one
Definition: factor.cpp:1462
mvec m
Definition: factor.cpp:771
Interface to class of symbolic functions.
Interface to GiNaC's initially known functions.
static std::vector< ex > cache_vec
static symbol x
int order
Interface to GiNaC's integration kernels for iterated integrals.
Interface to GiNaC's products of expressions.
Definition: add.cpp:38
const numeric I
Imaginary unit.
Definition: numeric.cpp:1433
const numeric pow(const numeric &x, const numeric &y)
Definition: numeric.h:251
bool is_discriminant_of_quadratic_number_field(const numeric &n)
Returns true if the integer n is either one or the discriminant of a quadratic number field.
const numeric bernoulli(const numeric &nn)
Bernoulli number.
Definition: numeric.cpp:2166
const numeric mod(const numeric &a, const numeric &b)
Modulus (in positive representation).
Definition: numeric.cpp:2328
const numeric abs(const numeric &x)
Absolute value.
Definition: numeric.cpp:2315
const ex _ex1
Definition: utils.cpp:385
ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) and b(X) in Z[X].
Definition: normal.cpp:1432
ex series_to_poly(const ex &e)
Convert the pseries object embedded in an expression to an ordinary polynomial in the expansion varia...
Definition: pseries.h:136
const numeric irem(const numeric &a, const numeric &b)
Numeric integer remainder.
Definition: numeric.cpp:2363
ex series(const ex &thisex, const ex &r, int order, unsigned options=0)
Definition: ex.h:781
numeric dirichlet_character(const numeric &n, const numeric &a, const numeric &N)
Defines a Dirichlet character through the Kronecker symbol.
ex diff(const ex &thisex, const symbol &s, unsigned nth=1)
Definition: ex.h:778
const numeric exp(const numeric &x)
Exponential function.
Definition: numeric.cpp:1439
numeric kronecker_symbol(const numeric &a, const numeric &n)
Returns the Kronecker symbol a: integer n: integer.
ex Bernoulli_polynomial(const numeric &k, const ex &x)
The Bernoulli polynomials.
const numeric factorial(const numeric &n)
Factorial combinatorial function.
Definition: numeric.cpp:2113
const numeric cos(const numeric &x)
Numeric cosine (trigonometric function).
Definition: numeric.cpp:1470
bool is_even(const numeric &x)
Definition: numeric.h:281
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
Definition: numeric.cpp:2341
const constant Pi("Pi", PiEvalf, "\\pi", domain::positive)
Pi.
Definition: constant.h:82
print_func< print_context >(&varidx::do_print). print_func< print_latex >(&varidx
Definition: idx.cpp:45
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(add, expairseq, print_func< print_context >(&add::do_print). print_func< print_latex >(&add::do_print_latex). print_func< print_csrc >(&add::do_print_csrc). print_func< print_tree >(&add::do_print_tree). print_func< print_python_repr >(&add::do_print_python_repr)) add
Definition: add.cpp:40
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
Definition: numeric.cpp:2404
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT_T(lst, basic, print_func< print_context >(&lst::do_print). print_func< print_tree >(&lst::do_print_tree)) template<> bool lst GINAC_BIND_UNARCHIVER(lst)
Specialization of container::info() for lst.
Definition: lst.cpp:42
_numeric_digits Digits
Accuracy in decimal digits.
Definition: numeric.cpp:2586
const numeric sin(const numeric &x)
Numeric sine (trigonometric function).
Definition: numeric.cpp:1461
ex normal(const ex &thisex)
Definition: ex.h:754
numeric primitive_dirichlet_character(const numeric &n, const numeric &a)
Defines a primitive Dirichlet character through the Kronecker symbol.
ex ifactor(const numeric &n)
Returns the decomposition of the positive integer n into prime numbers in the form lst( lst(p1,...
numeric generalised_Bernoulli_number(const numeric &k, const numeric &b)
The generalised Bernoulli number.
const ex _ex0
Definition: utils.cpp:369
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to class for extended truncated power series.
#define GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(classname, supername, options)
Macro for inclusion in the implementation of each registered class.
Definition: registrar.h:185
Interface to relations between expressions.
Interface to GiNaC's symbolic objects.
Interface to several small and furry utilities needed within GiNaC but not of any interest to the use...

This page is part of the GiNaC developer's reference. It was generated automatically by doxygen. For an introduction, see the tutorial.