integration_kernel.cpp

Go to the documentation of this file.

 
/*
 *  GiNaC Copyright (C) 1999-2026 Johannes Gutenberg University Mainz, Germany
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
#include "integration_kernel.h"
#include "add.h"
#include "mul.h"
#include "operators.h"
#include "power.h"
#include "relational.h"
#include "symbol.h"
#include "constant.h"
#include "numeric.h"
#include "function.h"
#include "pseries.h"
#include "utils.h"
#include "inifcns.h"
 
#include <iostream>
#include <stdexcept>
#include <cln/cln.h>
 
 
namespace GiNaC {
 
// anonymous namespace for helper function
namespace {

numeric kronecker_symbol_prime(const numeric & a, const numeric & n)
{
    if ( n == 1 ) {
        return 1;
    }
    else if ( n == -1 ) {
        if ( a >= 0 ) {
            return 1;
        }
        else {
            return -1;
        }
    }
    else if ( n == 2 ) {
        if ( GiNaC::smod(a,8) == 1 ) {
            return 1;
        }
        else if ( GiNaC::smod(a,8) == -1 ) {
            return 1;
        }
        else if  ( GiNaC::smod(a,8) == 3 ) {
            return -1;
        }
        else if ( GiNaC::smod(a,8) == -3 ) {
            return -1;
        }
        else {
            return 0;
        }
    }
 
    // n is an odd prime number
    return GiNaC::smod( pow(a,(n-1)/2), n);
}

numeric divisor_function(const numeric & n, const numeric & a, const numeric & b, const numeric & k)
{
    ex res = 0;
 
    for (numeric i1=1; i1<=n; i1++) {
        if ( irem(n,i1) == 0 ) {
            numeric ratio = n/i1;
            res += primitive_dirichlet_character(ratio,a) * primitive_dirichlet_character(i1,b) * pow(i1,k-1);
        }
    }
 
    return ex_to<numeric>(res);
}

numeric coefficient_a0(const numeric & k, const numeric & a, const numeric & b)
{
    ex conductor = abs(a);
 
    numeric a0;
    if ( conductor == 1 ) {
        a0 = -numeric(1,2)/k*generalised_Bernoulli_number(k,b);
    }
    else {
        a0 = 0;
    }
 
    return a0;
}

ex eisenstein_series(const numeric & k, const ex & q, const numeric & a, const numeric & b, const numeric & N)
{
    ex res = coefficient_a0(k,a,b);
 
    for (numeric i1=1; i1<N; i1++) {
        res += divisor_function(i1,a,b,k) * pow(q,i1);
    }
 
    return res;
}

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)
{
    int N_order_int = N_order.to_int();
 
    ex res = eisenstein_series(k,pow(q,K),a,b,iquo(N_order,K));
 
    res += Order(pow(q,N_order_int));
    res = res.series(q,N_order_int);
 
    return res;
}

ex B_eisenstein_series(const ex & q, const numeric & N_level, const numeric & K, const numeric & N_order)
{
    int N_order_int = N_order.to_int();
 
    ex res = eisenstein_series(2,q,1,1,N_order) - K*eisenstein_series(2,pow(q,K),1,1,iquo(N_order,K));
 
    res += Order(pow(q,N_order_int));
    res = res.series(q,N_order_int);
 
    return res;
}

struct subs_q_expansion : public map_function
{
    subs_q_expansion(const ex & arg_qbar, int arg_order) : qbar(arg_qbar), order(arg_order)
        {}
 
    ex operator()(const ex & e)
        {
            if ( is_a<Eisenstein_kernel>(e) || is_a<Eisenstein_h_kernel>(e) ) {
                return series_to_poly(e.series(qbar,order));
            }
            else {
                return e.map(*this);
            }
        }
 
    ex qbar;
    int order;
};

class Li_negative
{
    // ctors
public:
    Li_negative();
 
    // non-virtual functions 
public:
    ex get_symbolic_value(int n, const ex & x_val);
    ex get_numerical_value(int n, const ex & x_val);
 
    // member variables :
protected:
    static std::vector<ex> cache_vec;
    static symbol x;
};
 
 
Li_negative::Li_negative() {}
 
ex Li_negative::get_symbolic_value(int n, const ex & x_val)
{
    int n_cache = cache_vec.size();
 
    if ( n >= n_cache ) {
        for (int j=n_cache; j<=n; j++) {
            ex f;
            if ( j == 0 ) {
                f = x/(1-x);
            }
            else {
                f = normal( x*diff(cache_vec[j-1],x));
            }
            cache_vec.push_back( f );
        }
    }
 
    return cache_vec[n].subs(x==x_val);
}
 
ex Li_negative::get_numerical_value(int n, const ex & x_val)
{
    symbol x_symb("x_symb");
 
    ex f = this->get_symbolic_value(n,x_symb);
 
    ex res = f.subs(x_symb==x_val).evalf();
 
    return res;
}
 
// initialise static data members
std::vector<ex> Li_negative::cache_vec;
symbol Li_negative::x = symbol("x");
 
 
} // end of anonymous namespace

ex ifactor(const numeric & n)
{
    if ( !n.is_pos_integer() ) throw (std::runtime_error("ifactor(): argument not a positive integer"));
 
    lst p_lst, exp_lst;
 
    // implementation for small integers
    numeric n_temp = n;
    for (numeric p=2; p<=n; p++) {
        if ( p.info(info_flags::prime) ) {
            numeric exp_temp = 0;
            while ( irem(n_temp, p) == 0 ) {
                n_temp = n_temp/p;
                exp_temp++;
            }
            if ( exp_temp>0 ) {
                p_lst.append(p);
                exp_lst.append(exp_temp);
            }
        }
        if ( n_temp == 1 ) break;
    }
 
    if ( n_temp != 1 ) throw (std::runtime_error("ifactor(): probabilistic primality test failed"));
 
    lst res = {p_lst,exp_lst};
 
    return res;
}

bool is_discriminant_of_quadratic_number_field(const numeric & n)
{
    if ( n == 0 ) {
        return false;
    }
 
    if ( n == 1 ) {
        return true;
    }
 
    lst prime_factorisation = ex_to<lst>(ifactor(abs(n)));
    lst p_lst = ex_to<lst>(prime_factorisation.op(0));
    lst e_lst = ex_to<lst>(prime_factorisation.op(1));
 
    size_t n_primes = p_lst.nops();
 
    if ( n_primes > 0 ) {
        // take the last prime
        numeric p = ex_to<numeric>(p_lst.op(n_primes-1));
    
        if ( p.is_odd() ) {
            if ( e_lst.op(n_primes-1) != 1 ) {
                return false;
            }
 
            numeric pstar = p;
            if ( mod(p,4) == 3 ) {
                pstar = -p;
            }
            return is_discriminant_of_quadratic_number_field(n/pstar);
        }
    }
    // power of two now
    if ( (n==-4) || (n==-8) || (n==8) || (n==-32) || (n==32) || (n==-64) || (n==128) ) {
        return true; 
    }
 
    return false;
}

numeric kronecker_symbol(const numeric & a, const numeric & n)
{
    // case n=0 first, include kronecker_symbol(0,0)=0
    if ( n == 0 ) {
        if ( (a == 1) || (a == -1) ) {
            return 1;
        }
        else {
            return 0;
        }
    }
 
    numeric unit = 1;
    numeric n_pos = n;
    if ( n_pos<0 ) {
        unit = -1;
        n_pos = -n;
    }
 
    ex res = kronecker_symbol_prime(a,unit);
 
    numeric n_odd = n_pos;
    numeric alpha = 0;
    while ( n_odd.is_even() ) {
        alpha++;
        n_odd = n_odd/2;
    }
    if ( alpha>0 ) {
        res *= pow(kronecker_symbol_prime(a,2),alpha);
    }
 
    lst temp_lst = ex_to<lst>(ifactor(n_odd));
    lst prime_lst = ex_to<lst>(temp_lst.op(0));
    lst expo_lst = ex_to<lst>(temp_lst.op(1));
 
    for (auto it_p = prime_lst.begin(), it_e = expo_lst.begin(); it_p != prime_lst.end(); it_p++, it_e++) {
        res *= pow(kronecker_symbol_prime(a,ex_to<numeric>(*it_p)),ex_to<numeric>(*it_e));
    }
 
    return ex_to<numeric>(res);
}

numeric primitive_dirichlet_character(const numeric & n, const numeric & a)
{
    return kronecker_symbol(a,n);
}

numeric dirichlet_character(const numeric & n, const numeric & a, const numeric & N)
{
    if ( gcd(n,N) == 1 ) {
        return primitive_dirichlet_character(n,a);
    }
 
    return 0;
}

numeric generalised_Bernoulli_number(const numeric & k, const numeric & b)
{
    int k_int = k.to_int();
 
    symbol x("x");
 
    numeric conductor = abs(b);
 
    ex gen_fct = 0;
    for (numeric i1=1; i1<=conductor; i1++) {
        gen_fct += primitive_dirichlet_character(i1,b) * x*exp(i1*x)/(exp(conductor*x)-1);
    }
 
    gen_fct = series_to_poly(gen_fct.series(x,k_int+1));
 
    ex B = factorial(k) * gen_fct.coeff(x,k_int);
 
    return ex_to<numeric>(B);
}

ex Bernoulli_polynomial(const numeric & k, const ex & x)
{
    int k_int = k.to_int();
 
    symbol t("t");
 
    ex gen_fct = t*exp(x*t)/(exp(t)-1);
 
    gen_fct = series_to_poly(gen_fct.series(t,k_int+1));
 
    ex B = factorial(k) * gen_fct.coeff(t,k_int);
 
    return B;
}
 
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integration_kernel, basic,
  print_func<print_context>(&integration_kernel::do_print))
 
integration_kernel::integration_kernel() : inherited(), cache_step_size(100), series_vec()
{
}
 
int integration_kernel::compare_same_type(const basic &other) const
{
    return 0;
}
 
ex integration_kernel::series(const relational & r, int order, unsigned options) const
{
    if ( r.rhs() != 0 ) {
        throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
    }
 
    return Laurent_series(r.lhs(),order);
}

bool integration_kernel::has_trailing_zero(void) const
{
    if ( cln::zerop( series_coeff(0) ) ) {
        return false;
    }
 
    return true;
}

bool integration_kernel::is_numeric(void) const
{
    return true;
}

cln::cl_N integration_kernel::series_coeff(int i) const
{
    int n_vec = series_vec.size();
 
    if ( i >= n_vec ) {
        int N = cache_step_size*(i/cache_step_size+1);
 
        if ( uses_Laurent_series() ) {
            symbol x("x");
            // series_vec[0] gives coefficient of 1/z, series_vec[N-1] coefficient of z^(N-2),
            // thus expansion up to order (N-1) is required
            ex temp = Laurent_series(x, N-1);
            for (int j=n_vec; j<N; j++) {
                series_vec.push_back( ex_to<numeric>(temp.coeff(x,j-1).evalf()).to_cl_N() );
            }
        }
        else {
            for (int j=n_vec; j<N; j++) {
                series_vec.push_back( series_coeff_impl(j) );
            }
        }
    }
 
    return series_vec[i];
}

cln::cl_N integration_kernel::series_coeff_impl(int i) const
{
    if ( i == 1 ) {
        return 1;
    }
 
    return 0;
}

ex integration_kernel::Laurent_series(const ex & x, int order) const
{
    ex res = 0;
    for (int n=-1; n<order; n++) {
        res += numeric(series_coeff(n+1)) * pow(x,n);
    }
    res += Order(pow(x,order));
    res = res.series(x,order);
 
    return res;
}

ex  integration_kernel::get_numerical_value(const ex & lambda, int N_trunc) const
{
    return get_numerical_value_impl(lambda, 1, 0, N_trunc);
}

bool integration_kernel::uses_Laurent_series() const
{
    return false;
}

size_t integration_kernel::get_cache_size(void) const
{
    return series_vec.size();
}

void integration_kernel::set_cache_step(int cache_steps) const
{
    cache_step_size = cache_steps;
}

ex integration_kernel::get_series_coeff(int i) const
{
    return numeric(series_coeff(i));
}

ex  integration_kernel::get_numerical_value_impl(const ex & lambda, const ex & pre, int shift, int N_trunc) const
{
    cln::cl_N lambda_cln = ex_to<numeric>(lambda.evalf()).to_cl_N();
    cln::cl_N pre_cln = ex_to<numeric>(pre.evalf()).to_cl_N();
 
    cln::cl_F one = cln::cl_float(1, cln::float_format(Digits));
 
    cln::cl_N res = 0;
    cln::cl_N resbuf;
    cln::cl_N subexpr;
 
    if ( N_trunc == 0 ) {
        // sum until precision is reached
        bool flag_accidental_zero = false;
 
        int N = 0;
 
        do {
            resbuf = res;
     
            subexpr = series_coeff(N);
 
            res += pre_cln * subexpr * cln::expt(lambda_cln,N-1+shift);
 
            flag_accidental_zero = cln::zerop(subexpr);
 
            N++;
        } while ( (res != resbuf) || flag_accidental_zero );
    }
    else {
        // N_trunc > 0, sum up the first N_trunc terms
        for (int N=0; N<N_trunc; N++) {
            subexpr = series_coeff(N);
 
            res += pre_cln * subexpr * cln::expt(lambda_cln,N-1+shift);
        }
    }
 
    return numeric(res);
}
 
void integration_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "integration_kernel()";
}
 
GINAC_BIND_UNARCHIVER(integration_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(basic_log_kernel, integration_kernel,
                     print_func<print_context>(&basic_log_kernel::do_print))
 
basic_log_kernel::basic_log_kernel() : inherited()
{ 
}
 
int basic_log_kernel::compare_same_type(const basic &other) const
{
    return 0;
}
 
cln::cl_N basic_log_kernel::series_coeff_impl(int i) const
{
    if ( i == 0 ) {
        return 1;
    }
 
    return 0;
}
 
void basic_log_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "basic_log_kernel()";
}
 
GINAC_BIND_UNARCHIVER(basic_log_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(multiple_polylog_kernel, integration_kernel,
                     print_func<print_context>(&multiple_polylog_kernel::do_print))
 
multiple_polylog_kernel::multiple_polylog_kernel() : inherited(), z(_ex1)
{ 
}
 
multiple_polylog_kernel::multiple_polylog_kernel(const ex & arg_z) : inherited(), z(arg_z)
{
}
 
int multiple_polylog_kernel::compare_same_type(const basic &other) const
{
    const multiple_polylog_kernel &o = static_cast<const multiple_polylog_kernel &>(other);
 
    return z.compare(o.z);
}
 
size_t multiple_polylog_kernel::nops() const
{
    return 1;
}
 
ex multiple_polylog_kernel::op(size_t i) const
{
    if ( i != 0 ) {
        throw(std::range_error("multiple_polylog_kernel::op(): out of range"));
    }
 
    return z;
}
 
ex & multiple_polylog_kernel::let_op(size_t i)
{
    ensure_if_modifiable();
 
    if ( i != 0 ) {
        throw(std::range_error("multiple_polylog_kernel::let_op(): out of range"));
    }
 
    return z;
}
 
bool multiple_polylog_kernel::is_numeric(void) const
{
    return z.evalf().info(info_flags::numeric);
}
 
cln::cl_N multiple_polylog_kernel::series_coeff_impl(int i) const
{
    if ( i == 0 ) {
        return 0;
    }
 
    return -cln::expt(ex_to<numeric>(z.evalf()).to_cl_N(),-i);
}
 
void multiple_polylog_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "multiple_polylog_kernel(";
    z.print(c);
    c.s << ")";
}
 
GINAC_BIND_UNARCHIVER(multiple_polylog_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(ELi_kernel, integration_kernel,
                     print_func<print_context>(&ELi_kernel::do_print))
 
ELi_kernel::ELi_kernel() : inherited(), n(_ex0), m(_ex0), x(_ex0), y(_ex0)
{ 
}
 
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)
{
}
 
int ELi_kernel::compare_same_type(const basic &other) const
{
    const ELi_kernel &o = static_cast<const ELi_kernel &>(other);
    int cmpval;
 
    cmpval = n.compare(o.n);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = m.compare(o.m);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = x.compare(o.x);
    if ( cmpval) {
        return cmpval;
    }
 
    return y.compare(o.y);
}
 
size_t ELi_kernel::nops() const
{
    return 4;
}
 
ex ELi_kernel::op(size_t i) const
{
    switch (i) {
    case 0:
        return n;
    case 1:
        return m;
    case 2:
        return x;
    case 3:
        return y;
    default:
        throw (std::out_of_range("ELi_kernel::op() out of range"));
    }
}
 
ex & ELi_kernel::let_op(size_t i)
{
    ensure_if_modifiable();
 
    switch (i) {
    case 0:
        return n;
    case 1:
        return m;
    case 2:
        return x;
    case 3:
        return y;
    default:
        throw (std::out_of_range("ELi_kernel::let_op() out of range"));
    }
}
 
bool ELi_kernel::is_numeric(void) const
{
    return (n.info(info_flags::nonnegint) && m.info(info_flags::numeric) && x.evalf().info(info_flags::numeric) && y.evalf().info(info_flags::numeric));
}
 
cln::cl_N ELi_kernel::series_coeff_impl(int i) const
{
    if ( i == 0 ) {
        return 0;
    }
 
    int n_int = ex_to<numeric>(n).to_int();
    int m_int = ex_to<numeric>(m).to_int();
 
    cln::cl_N x_cln = ex_to<numeric>(x.evalf()).to_cl_N();
    cln::cl_N y_cln = ex_to<numeric>(y.evalf()).to_cl_N();
 
    cln::cl_N res_cln = 0;
 
    for (int j=1; j<=i; j++) {
        if ( (i % j) == 0 ) {
            int k = i/j;
 
            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);
        }
    }
 
    return res_cln;
}

ex  ELi_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
{
    return get_numerical_value_impl(qbar, 1, 1, N_trunc);
}
 
void ELi_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "ELi_kernel(";
    n.print(c);
    c.s << ",";
    m.print(c);
    c.s << ",";
    x.print(c);
    c.s << ",";
    y.print(c);
    c.s << ")";
}
 
GINAC_BIND_UNARCHIVER(ELi_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Ebar_kernel, integration_kernel,
                     print_func<print_context>(&Ebar_kernel::do_print))
 
Ebar_kernel::Ebar_kernel() : inherited(), n(_ex0), m(_ex0), x(_ex0), y(_ex0)
{ 
}
 
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)
{
}
 
int Ebar_kernel::compare_same_type(const basic &other) const
{
    const Ebar_kernel &o = static_cast<const Ebar_kernel &>(other);
    int cmpval;
 
    cmpval = n.compare(o.n);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = m.compare(o.m);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = x.compare(o.x);
    if ( cmpval) {
        return cmpval;
    }
 
    return y.compare(o.y);
}
 
size_t Ebar_kernel::nops() const
{
    return 4;
}
 
ex Ebar_kernel::op(size_t i) const
{
    switch (i) {
    case 0:
        return n;
    case 1:
        return m;
    case 2:
        return x;
    case 3:
        return y;
    default:
        throw (std::out_of_range("Ebar_kernel::op() out of range"));
    }
}
 
ex & Ebar_kernel::let_op(size_t i)
{
    ensure_if_modifiable();
 
    switch (i) {
    case 0:
        return n;
    case 1:
        return m;
    case 2:
        return x;
    case 3:
        return y;
    default:
        throw (std::out_of_range("Ebar_kernel::let_op() out of range"));
    }
}
 
bool Ebar_kernel::is_numeric(void) const
{
    return (n.info(info_flags::nonnegint) && m.info(info_flags::numeric) && x.evalf().info(info_flags::numeric) && y.evalf().info(info_flags::numeric));
}
 
cln::cl_N Ebar_kernel::series_coeff_impl(int i) const
{
    if ( i == 0 ) {
        return 0;
    }
 
    int n_int = ex_to<numeric>(n).to_int();
    int m_int = ex_to<numeric>(m).to_int();
 
    cln::cl_N x_cln = ex_to<numeric>(x.evalf()).to_cl_N();
    cln::cl_N y_cln = ex_to<numeric>(y.evalf()).to_cl_N();
 
    cln::cl_N res_cln = 0;
 
    for (int j=1; j<=i; j++) {
        if ( (i % j) == 0 ) {
            int k = i/j;
 
            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);
        }
    }
 
    return res_cln;
}

ex  Ebar_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
{
    return get_numerical_value_impl(qbar, 1, 1, N_trunc);
}
 
void Ebar_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "Ebar_kernel(";
    n.print(c);
    c.s << ",";
    m.print(c);
    c.s << ",";
    x.print(c);
    c.s << ",";
    y.print(c);
    c.s << ")";
}
 
GINAC_BIND_UNARCHIVER(Ebar_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Kronecker_dtau_kernel, integration_kernel,
                     print_func<print_context>(&Kronecker_dtau_kernel::do_print))
 
Kronecker_dtau_kernel::Kronecker_dtau_kernel() : inherited(), n(_ex0), z(_ex0), K(_ex1), C_norm(_ex1)
{ 
}
 
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)
{
}
 
int Kronecker_dtau_kernel::compare_same_type(const basic &other) const
{
    const Kronecker_dtau_kernel &o = static_cast<const Kronecker_dtau_kernel &>(other);
    int cmpval;
 
    cmpval = n.compare(o.n);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = z.compare(o.z);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = K.compare(o.K);
    if ( cmpval) {
        return cmpval;
    }
 
    return C_norm.compare(o.C_norm);
}
 
size_t Kronecker_dtau_kernel::nops() const
{
    return 4;
}
 
ex Kronecker_dtau_kernel::op(size_t i) const
{
    switch (i) {
    case 0:
        return n;
    case 1:
        return z;
    case 2:
        return K;
    case 3:
        return C_norm;
    default:
        throw (std::out_of_range("Kronecker_dtau_kernel::op() out of range"));
    }
}
 
ex & Kronecker_dtau_kernel::let_op(size_t i)
{
    ensure_if_modifiable();
 
    switch (i) {
    case 0:
        return n;
    case 1:
        return z;
    case 2:
        return K;
    case 3:
        return C_norm;
    default:
        throw (std::out_of_range("Kronecker_dtau_kernel::let_op() out of range"));
    }
}
 
bool Kronecker_dtau_kernel::is_numeric(void) const
{
    return (n.info(info_flags::nonnegint) && z.evalf().info(info_flags::numeric) && K.info(info_flags::posint) && C_norm.evalf().info(info_flags::numeric));
}
 
cln::cl_N Kronecker_dtau_kernel::series_coeff_impl(int i) const
{
    numeric n_num = ex_to<numeric>(n);
    int n_int = n_num.to_int();
 
    // case n=0
    if ( n_num == 0 ) {
        if ( i == 0 ) {
            ex res = -C_norm*K;
 
            return ex_to<numeric>(res.evalf()).to_cl_N();
        }
 
        return 0;
    }
 
    // case n=1
    if ( n_num == 1 ) {
        return 0;
    }
 
    // case n>1
    if ( i == 0 ) {
        ex res = C_norm*K / factorial(n_num-2) * bernoulli(n_num)/n_num;
 
        return ex_to<numeric>(res.evalf()).to_cl_N();
    }
 
    // n>1, i>0
 
    // if K>1 the variable i needs to be a multiple of K
    int K_int = ex_to<numeric>(K).to_int();
 
    if ( (i % K_int) != 0 ) {
        return 0;
    }
    int i_local = i/K_int;
 
    ex w = exp(ex_to<numeric>((2*Pi*I*z).evalf()));
    cln::cl_N w_cln = ex_to<numeric>(w).to_cl_N();
    cln::cl_N res_cln = 0;
    for (int j=1; j<=i_local; j++) {
        if ( (i_local % j) == 0 ) {
            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); 
        }
    }
    ex pre = -C_norm*K/factorial(n_num-2);
 
    return ex_to<numeric>(pre.evalf()).to_cl_N() * res_cln;
}

ex  Kronecker_dtau_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
{
    numeric n_num = ex_to<numeric>(n);
 
    if ( n_num == 0 ) {
        return 1;
    }
 
    // use the direct formula here
    if ( n_num == 1 ) {
        ex wbar = exp(ex_to<numeric>((2*Pi*I*z).evalf()));
        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));
 
        return ex_to<numeric>(res.evalf());
    }
 
    ex pre = pow(2*Pi*I,n_num)/C_norm/K/(n_num-1);
 
    return get_numerical_value_impl(qbar, pre, 1, N_trunc);
}
 
void Kronecker_dtau_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "Kronecker_dtau_kernel(";
    n.print(c);
    c.s << ",";
    z.print(c);
    c.s << ",";
    K.print(c);
    c.s << ",";
    C_norm.print(c);
    c.s << ")";
}
 
GINAC_BIND_UNARCHIVER(Kronecker_dtau_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Kronecker_dz_kernel, integration_kernel,
                     print_func<print_context>(&Kronecker_dz_kernel::do_print))
 
Kronecker_dz_kernel::Kronecker_dz_kernel() : inherited(), n(_ex0), z_j(_ex0), tau(_ex0), K(_ex1), C_norm(_ex1)
{ 
}
 
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)
{
}
 
int Kronecker_dz_kernel::compare_same_type(const basic &other) const
{
    const Kronecker_dz_kernel &o = static_cast<const Kronecker_dz_kernel &>(other);
    int cmpval;
 
    cmpval = n.compare(o.n);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = z_j.compare(o.z_j);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = tau.compare(o.tau);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = K.compare(o.K);
    if ( cmpval) {
        return cmpval;
    }
 
    return C_norm.compare(o.C_norm);
}
 
size_t Kronecker_dz_kernel::nops() const
{
    return 5;
}
 
ex Kronecker_dz_kernel::op(size_t i) const
{
    switch (i) {
    case 0:
        return n;
    case 1:
        return z_j;
    case 2:
        return tau;
    case 3:
        return K;
    case 4:
        return C_norm;
    default:
        throw (std::out_of_range("Kronecker_dz_kernel::op() out of range"));
    }
}
 
ex & Kronecker_dz_kernel::let_op(size_t i)
{
    ensure_if_modifiable();
 
    switch (i) {
    case 0:
        return n;
    case 1:
        return z_j;
    case 2:
        return tau;
    case 3:
        return K;
    case 4:
        return C_norm;
    default:
        throw (std::out_of_range("Kronecker_dz_kernel::let_op() out of range"));
    }
}
 
bool Kronecker_dz_kernel::is_numeric(void) const
{
    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));
}
 
cln::cl_N Kronecker_dz_kernel::series_coeff_impl(int i) const
{
    numeric n_num = ex_to<numeric>(n);
 
    ex w_j_inv = exp(ex_to<numeric>((-2*Pi*I*z_j).evalf()));
    cln::cl_N w_j_inv_cln = ex_to<numeric>(w_j_inv).to_cl_N();
 
    ex qbar = exp(ex_to<numeric>((2*Pi*I*K*tau).evalf()));
 
    // case n=0
    if ( n_num == 0 ) {
        return 0;
    }
 
    // case n=1
    if ( n_num == 1 ) {
        if ( i == 1 ) {
            return ex_to<numeric>((C_norm * 2*Pi*I).evalf()).to_cl_N();
        }
 
        return 0;
    }
 
    // case n=2
    if ( n_num == 2 ) {
        if ( ex_to<numeric>(z_j.evalf()).is_zero() ) {
            if ( i == 0 ) {
                return ex_to<numeric>((C_norm).evalf()).to_cl_N();
            }
            else if ( i == 1 ) {
                return 0;
            }
            else {
                ex res = -bernoulli(i)/numeric(i);
                if ( numeric(i).is_even() ) {
                    Ebar_kernel Ebar = Ebar_kernel( 1-i, 0, numeric(1), numeric(1) );
                    res += Ebar.get_numerical_value(qbar);
                }
 
                res *= -pow(2*Pi*I,i)*C_norm/factorial(i-1);
 
                return ex_to<numeric>(res.evalf()).to_cl_N();
            }
        }
        else {
            // z_j is not zero
            if ( i == 0 ) {
                return 0;
            }
            else {
                Li_negative my_Li_negative;
 
                ex res = 0;
                if ( i == 1 ) {
                    res = numeric(1,2);
                }
 
                Ebar_kernel Ebar = Ebar_kernel( 1-i, 0, w_j_inv, numeric(1) );
 
                res += my_Li_negative.get_numerical_value(i-1,w_j_inv) + Ebar.get_numerical_value(qbar);
 
                res *= -pow(2*Pi*I,i)*C_norm/factorial(i-1);
 
                return ex_to<numeric>(res.evalf()).to_cl_N();
            }
        }
    }
 
    // case n>2
    ex res = 0;
    if ( i == 1 ) {
        res += - bernoulli(n_num-1)/(n_num-1);
    }
    if ( i > 0 ) {
        if ( ex_to<numeric>(z_j.evalf()).is_zero() ) {
            if ( (numeric(i)+n_num).is_even() ) {
                Ebar_kernel Ebar = Ebar_kernel( 1-i, 2-n_num, numeric(1), numeric(1) );
 
                res += pow(2*Pi*I,i-1)/factorial(i-1) * Ebar.get_numerical_value(qbar);
            }
        }
        else {
            // z_j is not zero
            Ebar_kernel Ebar = Ebar_kernel( 1-i, 2-n_num, w_j_inv, numeric(1) );
 
            res += pow(2*Pi*I,i-1)/factorial(i-1) * Ebar.get_numerical_value(qbar);
        }
    }
 
    res *= - C_norm * 2*Pi*I/factorial(n_num-2);
 
    return ex_to<numeric>(res.evalf()).to_cl_N();
}

ex  Kronecker_dz_kernel::get_numerical_value(const ex & z, int N_trunc) const
{
    numeric n_num = ex_to<numeric>(n);
 
    if ( n_num == 1 ) {
        return 1;
    }
 
    ex pre = pow(2*Pi*I,n-2)/C_norm;
 
    return get_numerical_value_impl(z, pre, 0, N_trunc);
}
 
void Kronecker_dz_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "Kronecker_dz_kernel(";
    n.print(c);
    c.s << ",";
    z_j.print(c);
    c.s << ",";
    tau.print(c);
    c.s << ",";
    K.print(c);
    c.s << ",";
    C_norm.print(c);
    c.s << ")";
}
 
GINAC_BIND_UNARCHIVER(Kronecker_dz_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eisenstein_kernel, integration_kernel,
                     print_func<print_context>(&Eisenstein_kernel::do_print))
 
Eisenstein_kernel::Eisenstein_kernel() : inherited(), k(_ex0), N(_ex0), a(_ex0), b(_ex0), K(_ex0), C_norm(_ex1)
{ 
}
 
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)
{ 
}
 
int Eisenstein_kernel::compare_same_type(const basic &other) const
{
    const Eisenstein_kernel &o = static_cast<const Eisenstein_kernel &>(other);
    int cmpval;
 
    cmpval = k.compare(o.k);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = N.compare(o.N);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = a.compare(o.a);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = b.compare(o.b);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = K.compare(o.K);
    if ( cmpval) {
        return cmpval;
    }
 
    return C_norm.compare(o.C_norm);
}

ex Eisenstein_kernel::series(const relational & r, int order, unsigned options) const
{
    if ( r.rhs() != 0 ) {
        throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
    }
 
    ex qbar = r.lhs();
    ex res = q_expansion_modular_form(qbar, order);
    res = res.series(qbar,order);
 
    return res;
}
 
size_t Eisenstein_kernel::nops() const
{
    return 6;
}
 
ex Eisenstein_kernel::op(size_t i) const
{
    switch (i) {
    case 0:
        return k;
    case 1:
        return N;
    case 2:
        return a;
    case 3:
        return b;
    case 4:
        return K;
    case 5:
        return C_norm;
    default:
        throw (std::out_of_range("Eisenstein_kernel::op() out of range"));
    }
}
 
ex & Eisenstein_kernel::let_op(size_t i)
{
    ensure_if_modifiable();
 
    switch (i) {
    case 0:
        return k;
    case 1:
        return N;
    case 2:
        return a;
    case 3:
        return b;
    case 4:
        return K;
    case 5:
        return C_norm;
    default:
        throw (std::out_of_range("Eisenstein_kernel::let_op() out of range"));
    }
}
 
bool Eisenstein_kernel::is_numeric(void) const
{
    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));
}
 
ex Eisenstein_kernel::Laurent_series(const ex & x, int order) const
{
    ex res = C_norm * q_expansion_modular_form(x, order+1)/x;
    res = res.series(x,order);
 
    return res;
}

ex  Eisenstein_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
{
    ex pre = numeric(1)/C_norm;
 
    return get_numerical_value_impl(qbar, pre, 1, N_trunc);
}
 
bool Eisenstein_kernel::uses_Laurent_series() const
{
    return true;
}
 
ex Eisenstein_kernel::q_expansion_modular_form(const ex & q, int order) const
{
    numeric k_num = ex_to<numeric>(k);
    numeric N_num = ex_to<numeric>(N);
    numeric a_num = ex_to<numeric>(a);
    numeric b_num = ex_to<numeric>(b);
    numeric K_num = ex_to<numeric>(K);
 
    if ( (k==2) && (a==1) && (b==1) ) {
        return B_eisenstein_series(q, N_num, K_num, order);
    }
 
    return E_eisenstein_series(q, k_num, N_num, a_num, b_num, K_num, order);
}
 
void Eisenstein_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "Eisenstein_kernel(";
    k.print(c);
    c.s << ",";
    N.print(c);
    c.s << ",";
    a.print(c);
    c.s << ",";
    b.print(c);
    c.s << ",";
    K.print(c);
    c.s << ",";
    C_norm.print(c);
    c.s << ")";
}
 
GINAC_BIND_UNARCHIVER(Eisenstein_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eisenstein_h_kernel, integration_kernel,
                     print_func<print_context>(&Eisenstein_h_kernel::do_print))
 
Eisenstein_h_kernel::Eisenstein_h_kernel() : inherited(), k(_ex0), N(_ex0), r(_ex0), s(_ex0), C_norm(_ex1)
{ 
}
 
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)
{ 
}
 
int Eisenstein_h_kernel::compare_same_type(const basic &other) const
{
    const Eisenstein_h_kernel &o = static_cast<const Eisenstein_h_kernel &>(other);
    int cmpval;
 
    cmpval = k.compare(o.k);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = N.compare(o.N);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = r.compare(o.r);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = s.compare(o.s);
    if ( cmpval) {
        return cmpval;
    }
 
    return C_norm.compare(o.C_norm);
}

ex Eisenstein_h_kernel::series(const relational & r, int order, unsigned options) const
{
    if ( r.rhs() != 0 ) {
        throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
    }
 
    ex qbar = r.lhs();
    ex res = q_expansion_modular_form(qbar, order);
    res = res.series(qbar,order);
 
    return res;
}
 
size_t Eisenstein_h_kernel::nops() const
{
    return 5;
}
 
ex Eisenstein_h_kernel::op(size_t i) const
{
    switch (i) {
    case 0:
        return k;
    case 1:
        return N;
    case 2:
        return r;
    case 3:
        return s;
    case 4:
        return C_norm;
    default:
        throw (std::out_of_range("Eisenstein_h_kernel::op() out of range"));
    }
}
 
ex & Eisenstein_h_kernel::let_op(size_t i)
{
    ensure_if_modifiable();
 
    switch (i) {
    case 0:
        return k;
    case 1:
        return N;
    case 2:
        return r;
    case 3:
        return s;
    case 4:
        return C_norm;
    default:
        throw (std::out_of_range("Eisenstein_h_kernel::let_op() out of range"));
    }
}
 
bool Eisenstein_h_kernel::is_numeric(void) const
{
    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));
}
 
ex Eisenstein_h_kernel::Laurent_series(const ex & x, int order) const
{
    ex res = C_norm * q_expansion_modular_form(x, order+1)/x;
    res = res.series(x,order);
 
    return res;
}

ex  Eisenstein_h_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
{
    ex pre = numeric(1)/C_norm;
 
    return get_numerical_value_impl(qbar, pre, 1, N_trunc);
}
 
bool Eisenstein_h_kernel::uses_Laurent_series() const
{
    return true;
}

ex Eisenstein_h_kernel::coefficient_a0(const numeric & k, const numeric & r, const numeric & s, const numeric & N) const
{
    if ( k == 1 ) {
        if ( irem(s,N) != 0 ) {
            return numeric(1,4) - mod(s,N)/numeric(2)/N;
        }
        else if ( (irem(r,N)==0) && (irem(s,N)==0) ) {
            return 0;
        }
        else {
            return I*numeric(1,4)*cos(Pi*mod(r,N)/N)/sin(Pi*mod(r,N)/N);
        }
    }
 
    // case k > 1
    return -Bernoulli_polynomial(k,mod(s,N)/N)/numeric(2)/k;
}

ex Eisenstein_h_kernel::coefficient_an(const numeric & n, const numeric & k, const numeric & r, const numeric & s, const numeric & N) const
{
    ex res = 0;
 
    for (numeric m=1; m<=n; m++) {
        if ( irem(n,m) == 0 ) {
            for (numeric c1=0; c1<N; c1++) {
                numeric c2 = n/m;
                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));
            }
        }
    }
 
    return res/numeric(2)/pow(N,k);
}
 
ex Eisenstein_h_kernel::q_expansion_modular_form(const ex & q, int N_order) const
{
    numeric N_order_num = numeric(N_order);
 
    numeric k_num = ex_to<numeric>(k);
    numeric r_num = ex_to<numeric>(r);
    numeric s_num = ex_to<numeric>(s);
    numeric N_num = ex_to<numeric>(N);
 
    ex res = coefficient_a0(k_num,r_num,s_num,N_num);
 
    for (numeric i1=1; i1<N_order_num; i1++) {
        res += coefficient_an(i1,k_num,r_num,s_num,N_num) * pow(q,i1);
    }
 
    res += Order(pow(q,N_order));
    res = res.series(q,N_order);
 
    return res;
}
 
void Eisenstein_h_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "Eisenstein_h_kernel(";
    k.print(c);
    c.s << ",";
    N.print(c);
    c.s << ",";
    r.print(c);
    c.s << ",";
    s.print(c);
    c.s << ",";
    C_norm.print(c);
    c.s << ")";
}
 
GINAC_BIND_UNARCHIVER(Eisenstein_h_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(modular_form_kernel, integration_kernel,
                     print_func<print_context>(&modular_form_kernel::do_print))
 
modular_form_kernel::modular_form_kernel() : inherited(), k(_ex0), P(_ex0), C_norm(_ex1)
{ 
}
 
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)
{ 
}
 
int modular_form_kernel::compare_same_type(const basic &other) const
{
    const modular_form_kernel &o = static_cast<const modular_form_kernel &>(other);
    int cmpval;
 
    cmpval = k.compare(o.k);
    if ( cmpval) {
        return cmpval;
    }
 
    cmpval = P.compare(o.P);
    if ( cmpval) {
        return cmpval;
    }
 
    return C_norm.compare(o.C_norm);
}

ex modular_form_kernel::series(const relational & r, int order, unsigned options) const
{
    if ( r.rhs() != 0 ) {
        throw (std::runtime_error("integration_kernel::series: non-zero expansion point not implemented"));
    }
 
    ex qbar = r.lhs();
 
    subs_q_expansion do_subs_q_expansion(qbar, order);
 
    ex res = do_subs_q_expansion(P).series(qbar,order);
    res += Order(pow(qbar,order));
    res = res.series(qbar,order);
 
    return res;
}
 
size_t modular_form_kernel::nops() const
{
    return 3;
}
 
ex modular_form_kernel::op(size_t i) const
{
    switch (i) {
    case 0:
        return k;
    case 1:
        return P;
    case 2:
        return C_norm;
    default:
        throw (std::out_of_range("modular_form_kernel::op() out of range"));
    }
}
 
ex & modular_form_kernel::let_op(size_t i)
{
    ensure_if_modifiable();
 
    switch (i) {
    case 0:
        return k;
    case 1:
        return P;
    case 2:
        return C_norm;
    default:
        throw (std::out_of_range("modular_form_kernel::let_op() out of range"));
    }
}
 
bool modular_form_kernel::is_numeric(void) const
{
    bool flag = (k.info(info_flags::nonnegint) && C_norm.evalf().info(info_flags::numeric));
    if ( !flag ) {
        return false;
    }
 
    symbol qbar("qbar");
 
    // test with a random number and random expansion
    return series_to_poly(q_expansion_modular_form(qbar,18)).subs(qbar==numeric(1,937)).evalf().info(info_flags::numeric);
}
 
ex modular_form_kernel::Laurent_series(const ex & qbar, int order) const
{
    ex res = series_to_poly(q_expansion_modular_form(qbar,order+1));
    res = C_norm * res/qbar;
    res = res.series(qbar,order);
    return res;
}

ex  modular_form_kernel::get_numerical_value(const ex & qbar, int N_trunc) const
{
    ex pre = numeric(1)/C_norm;
 
    return get_numerical_value_impl(qbar, pre, 1, N_trunc);
}
 
bool modular_form_kernel::uses_Laurent_series() const
{
    return true;
}
 
ex modular_form_kernel::q_expansion_modular_form(const ex & q, int N_order) const
{
    return this->series(q==0,N_order);
}
 
void modular_form_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "modular_form_kernel(";
    k.print(c);
    c.s << ",";
    P.print(c);
    c.s << ",";
    C_norm.print(c);
    c.s << ")";
}
 
GINAC_BIND_UNARCHIVER(modular_form_kernel);
 
 
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(user_defined_kernel, integration_kernel,
                     print_func<print_context>(&user_defined_kernel::do_print))
 
user_defined_kernel::user_defined_kernel() : inherited(), f(_ex0), x(_ex0)
{ 
}
 
user_defined_kernel::user_defined_kernel(const ex & arg_f, const ex & arg_x) : inherited(), f(arg_f), x(arg_x)
{ 
}
 
int user_defined_kernel::compare_same_type(const basic &other) const
{
    const user_defined_kernel &o = static_cast<const user_defined_kernel &>(other);
    int cmpval;
 
    cmpval = f.compare(o.f);
    if ( cmpval) {
        return cmpval;
    }
 
    return x.compare(o.x);
}
 
size_t user_defined_kernel::nops() const
{
    return 2;
}
 
ex user_defined_kernel::op(size_t i) const
{
    switch (i) {
    case 0:
        return f;
    case 1:
        return x;
    default:
        throw (std::out_of_range("user_defined_kernel::op() out of range"));
    }
}
 
ex & user_defined_kernel::let_op(size_t i)
{
    ensure_if_modifiable();
 
    switch (i) {
    case 0:
        return f;
    case 1:
        return x;
    default:
        throw (std::out_of_range("user_defined_kernel::let_op() out of range"));
    }
}
 
bool user_defined_kernel::is_numeric(void) const
{
    // test with a random number
    return f.subs(x==numeric(1,937)).evalf().info(info_flags::numeric);
}
 
ex user_defined_kernel::Laurent_series(const ex & x_up, int order) const
{
    ex res = f.series(x,order).subs(x==x_up);
 
    return res;
}
 
bool user_defined_kernel::uses_Laurent_series() const
{
    return true;
}
 
void user_defined_kernel::do_print(const print_context & c, unsigned level) const
{
    c.s << "user_defined_kernel(";
    f.print(c);
    c.s << ",";
    x.print(c);
    c.s << ")";
}
 
GINAC_BIND_UNARCHIVER(user_defined_kernel);
 
} // namespace GiNaC