static ex tanh_eval(const ex & x)
{
	if (x.info(info_flags::numeric)) {

		// tanh(0) -> 0
		if (x.is_zero())
			return _ex0;

		// tanh(float) -> float
		if (!x.info(info_flags::crational))
			return tanh(ex_to<numeric>(x));

		// tanh() is odd
		if (x.info(info_flags::negative))
			return -tanh(-x);
	}
	
	if ((x/Pi).info(info_flags::numeric) &&
		ex_to<numeric>(x/Pi).real().is_zero())  // tanh(I*x) -> I*tan(x);
		return I*tan(x/I);
	
	if (is_exactly_a<function>(x)) {
		const ex &t = x.op(0);

		// tanh(atanh(x)) -> x
		if (is_ex_the_function(x, atanh))
			return t;

		// tanh(asinh(x)) -> x/sqrt(1+x^2)
		if (is_ex_the_function(x, asinh))
			return t*power(_ex1+power(t,_ex2),_ex_1_2);

		// tanh(acosh(x)) -> sqrt(x-1)*sqrt(x+1)/x
		if (is_ex_the_function(x, acosh))
			return sqrt(t-_ex1)*sqrt(t+_ex1)*power(t,_ex_1);
	}
	
	return tanh(x).hold();
}

static ex cosh_eval(const ex & x)
{
	if (x.info(info_flags::numeric)) {

		// cosh(0) -> 1
		if (x.is_zero())
			return _ex1;

		// cosh(float) -> float
		if (!x.info(info_flags::crational))
			return cosh(ex_to<numeric>(x));

		// cosh() is even
		if (x.info(info_flags::negative))
			return cosh(-x);
	}
	
	if ((x/Pi).info(info_flags::numeric) &&
		ex_to<numeric>(x/Pi).real().is_zero())  // cosh(I*x) -> cos(x)
		return cos(x/I);
	
	if (is_exactly_a<function>(x)) {
		const ex &t = x.op(0);

		// cosh(acosh(x)) -> x
		if (is_ex_the_function(x, acosh))
			return t;

		// cosh(asinh(x)) -> sqrt(1+x^2)
		if (is_ex_the_function(x, asinh))
			return sqrt(_ex1+power(t,_ex2));

		// cosh(atanh(x)) -> 1/sqrt(1-x^2)
		if (is_ex_the_function(x, atanh))
			return power(_ex1-power(t,_ex2),_ex_1_2);
	}
	
	return cosh(x).hold();
}

//
//  On return, lags is a GiNac::lst, where each element is a GiNac::lst
//  of length 4 containing {lagsym, variable_index + 1, var, lag_time}
//
void VectorField::convert_delay_to_lagvalue(ex& f, lst &lags)
{
    symbol t(IndependentVariable);
    exset dlist;
    f.find(delay(wild(1),wild(2)),dlist);
    for (exset::const_iterator iter = dlist.begin(); iter != dlist.end(); ++iter) {
        ex delayfunc = *iter;
        ex delayexpr = delayfunc.op(0);
        lst vars = FindVarsInEx(delayexpr);
        ex del = delayfunc.op(1);
        for (lst::const_iterator iter = vars.begin(); iter != vars.end(); ++iter) {
            ostringstream os;
            lst tmp;
            os << lags.nops() + 1;
            symbol lagsym("lag" + os.str());
            int vindex = FindVar(ex_to<symbol>(*iter));
            delayexpr = delayexpr.subs(*iter == lagsym);
            tmp = {lagsym, vindex + 1, *iter, del};
            lags.append(tmp);
        }
        f = f.subs(delayfunc == delayexpr);
    }
}

bealab::function<double(double)> makefun( const ex& fun, const symbol& x )
{
	FUNCP_CUBA fp;
	ex fun1 = fun.subs( Pi==pi );
	compile_ex( lst(fun1), lst(x), fp );
	return [fp]( double x )
	{
		int ndim  = 1;
		int ncomp = 1;
		double res;
		fp( &ndim, &x, &ncomp, &res );
		return res;
	};
}

static ex tan_series(const ex &x,
                     const relational &rel,
                     int order,
                     unsigned options)
{
	GINAC_ASSERT(is_a<symbol>(rel.lhs()));
	// method:
	// Taylor series where there is no pole falls back to tan_deriv.
	// On a pole simply expand sin(x)/cos(x).
	const ex x_pt = x.subs(rel, subs_options::no_pattern);
	if (!(2*x_pt/Pi).info(info_flags::odd))
		throw do_taylor();  // caught by function::series()
	// if we got here we have to care for a simple pole
	return (sin(x)/cos(x)).series(rel, order, options);
}

static ex eta_eval(const ex &x, const ex &y)
{
	// trivial:  eta(x,c) -> 0  if c is real and positive
	if (x.info(info_flags::positive) || y.info(info_flags::positive))
		return _ex0;

	if (x.info(info_flags::numeric) &&	y.info(info_flags::numeric)) {
		// don't call eta_evalf here because it would call Pi.evalf()!
		const numeric nx = ex_to<numeric>(x);
		const numeric ny = ex_to<numeric>(y);
		const numeric nxy = ex_to<numeric>(x*y);
		int cut = 0;
		if (nx.is_real() && nx.is_negative())
			cut -= 4;
		if (ny.is_real() && ny.is_negative())
			cut -= 4;
		if (nxy.is_real() && nxy.is_negative())
			cut += 4;
		return (I/4)*Pi*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)-
		                 (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut);
	}
	
	return eta(x,y).hold();
}

static ex eta_evalf(const ex &x, const ex &y)
{
	// It seems like we basically have to replicate the eval function here,
	// since the expression might not be fully evaluated yet.
	if (x.info(info_flags::positive) || y.info(info_flags::positive))
		return _ex0;

	if (x.info(info_flags::numeric) &&	y.info(info_flags::numeric)) {
		const numeric nx = ex_to<numeric>(x);
		const numeric ny = ex_to<numeric>(y);
		const numeric nxy = ex_to<numeric>(x*y);
		int cut = 0;
		if (nx.is_real() && nx.is_negative())
			cut -= 4;
		if (ny.is_real() && ny.is_negative())
			cut -= 4;
		if (nxy.is_real() && nxy.is_negative())
			cut += 4;
		return evalf(I/4*Pi)*((csgn(-imag(nx))+1)*(csgn(-imag(ny))+1)*(csgn(imag(nxy))+1)-
		                      (csgn(imag(nx))+1)*(csgn(imag(ny))+1)*(csgn(-imag(nxy))+1)+cut);
	}

	return eta(x,y).hold();
}

static ex zeta1_series(const ex& m, const relational& rel, int order, unsigned options)
{
	// use taylor expansion everywhere except at the singularity at 1
	const numeric val = ex_to<numeric>(m.subs(rel, subs_options::no_pattern));
	if (val != 1)
		throw do_taylor();  // caught by function::series()
	// at 1, use the expansion with the stieltjes-constants
	ex ser = 1/(m-1);
	numeric fac = 1;
	for (int n = 0; n <= order; ++n) {
	  fac = fac.mul(n+1);
	  ser += pow(-1, n) * stieltjes(n) * pow(m-1, n) * fac.inverse();
	}
	return ser.series(rel, order, options);
}

static ex zeta1_evalf(const ex& x, PyObject* parent)
{
        /*
	if (is_exactly_a<lst>(x) && (x.nops()>1)) {

		// multiple zeta value
		const int count = x.nops();
		const lst& xlst = ex_to<lst>(x);
		std::vector<int> r(count);

		// check parameters and convert them
		auto it1 = xlst.begin();
		auto it2 = r.begin();
		do {
			if (!(*it1).info(info_flags::posint)) {
				return zeta(x).hold();
			}
			*it2 = ex_to<numeric>(*it1).to_int();
			it1++;
			it2++;
		} while (it2 != r.end());

		// check for divergence
		if (r[0] == 1) {
			return zeta(x).hold();
		}

		// decide on summation algorithm
		// this is still a bit clumsy
		int limit = 10;
		if ((r[0] < limit) || ((count > 3) && (r[1] < limit/2))) {
			return numeric(zeta_do_sum_Crandall(r));
		} else {
			return numeric(zeta_do_sum_simple(r));
		}
	}*/

	// single zeta value
	if (x == 1) {
		return UnsignedInfinity;
	} else	if (is_exactly_a<numeric>(x)) {
		try {
			return zeta(ex_to<numeric>(x.evalf(0, parent)));
		} catch (const dunno &e) { }
	}

	return zeta(x).hold();
}

static ex Order_eval(const ex & x)
{
	if (is_exactly_a<numeric>(x)) {
		// O(c) -> O(1) or 0
		if (!x.is_zero())
			return Order(_ex1).hold();
		else
			return _ex0;
	} else if (is_exactly_a<mul>(x)) {
		const mul &m = ex_to<mul>(x);
		// O(c*expr) -> O(expr)
		if (is_exactly_a<numeric>(m.op(m.nops() - 1)))
			return Order(x / m.op(m.nops() - 1)).hold();
	}
	return Order(x).hold();
}

void
Model_block::add_objective(const ex &obj, const ex &obj_eq, const ex &lambda, int l)
{
    m_obj_var = obj;
    m_obj_eq = obj_eq;
    m_obj_eq_in = obj_eq;
    m_obj_lm_in = lambda;
    if (m_obj_lm_in) {
        m_obj_lm = m_obj_lm_in;
    } else {
        m_obj_lm = ex("lambda__" + m_name + "_" + obj.str(DROP_IDX | DROP_T), 0);
        m_obj_lm = add_idx(m_obj_lm, m_i1);
        m_obj_lm = add_idx(m_obj_lm, m_i2);
        m_redlm.insert(m_obj_lm);
    }
    m_obj_line = l;
}

static ex zeta2_eval(const ex& m, const ex& s_)
{
	if (is_exactly_a<lst>(s_)) {
		const lst& s = ex_to<lst>(s_);
		for (const auto & elem : s) {
			if ((elem).info(info_flags::positive)) {
				continue;
			}
			return zeta(m, s_).hold();
		}
		return zeta(m);
	} else if (s_.info(info_flags::positive)) {
		return zeta(m);
	}

	return zeta(m, s_).hold();
}

static void zeta1_print_latex(const ex& m_, const print_context& c)
{
	c.s << "\\zeta(";
	if (is_a<lst>(m_)) {
		const lst& m = ex_to<lst>(m_);
		auto it = m.begin();
		(*it).print(c);
		it++;
		for (; it != m.end(); it++) {
			c.s << ",";
			(*it).print(c);
		}
	} else {
		m_.print(c);
	}
	c.s << ")";
}

static ex binomial_sym(const ex & x, const numeric & y)
{
	if (y.is_integer()) {
		if (y.is_nonneg_integer()) {
			const unsigned N = y.to_int();
			if (N == 0) return _ex1;
			if (N == 1) return x;
			ex t = x.expand();
			for (unsigned i = 2; i <= N; ++i)
				t = (t * (x + i - y - 1)).expand() / i;
			return t;
		} else
			return _ex0;
	}

	return binomial(x, y).hold();
}

static void
collect_term(ex_collect_priv_t& ec, const ex& e, const exvector& vars)
{
	if (e.is_zero())
		return;
	static const ex ex1(1);
	exp_vector_t key(vars.size());
	ex pre_coeff = e;
	for (std::size_t i = 0; i < vars.size(); ++i) {
		const int var_i_pow = pre_coeff.degree(vars[i]);
		key[i] = var_i_pow;
		pre_coeff = pre_coeff.coeff(vars[i], var_i_pow);
	}
	auto i = ec.find(key);
	if (i != ec.end())
		i->second += pre_coeff;
	else
		ec.insert(ex_collect_priv_t::value_type(key, pre_coeff));
}

void VectorField::CheckForDelay(const ex& f)
    {
    lst occurrences;
    if (f.find(delay(wild(1),wild(2)),occurrences))
        {
        IsDelay = true;
        for (lst::const_iterator iter = occurrences.begin(); iter != occurrences.end(); ++iter)
            {
            ex del = iter->op(1);
            AddDelay(del);
            if (del.has(IndVar))
                HasNonconstantDelay = true;  // time-dependent delay
            for (lst::const_iterator viter = varname_list.begin(); viter != varname_list.end(); ++viter)
                {
                if (del.has(*viter))
                    HasNonconstantDelay = true; // state-dependent delay
                }
            }
        }
    }

static ex log_expand(const ex & arg, unsigned options)
{
	if ((options & expand_options::expand_transcendental)
		&& is_exactly_a<mul>(arg) && !arg.info(info_flags::indefinite)) {
		exvector sumseq;
		exvector prodseq;
		sumseq.reserve(arg.nops());
		prodseq.reserve(arg.nops());
		bool possign=true;

		// searching for positive/negative factors
		for (const_iterator i = arg.begin(); i != arg.end(); ++i) {
			ex e;
			if (options & expand_options::expand_function_args)
				e=i->expand(options);
			else
				e=*i;
			if (e.info(info_flags::positive))
				sumseq.push_back(log(e));
			else if (e.info(info_flags::negative)) {
				sumseq.push_back(log(-e));
				possign = !possign;
			} else
				prodseq.push_back(e);
		}

		if (sumseq.size() > 0) {
			ex newarg;
			if (options & expand_options::expand_function_args)
				newarg=((possign?_ex1:_ex_1)*mul(prodseq)).expand(options);
			else {
				newarg=(possign?_ex1:_ex_1)*mul(prodseq);
				ex_to<basic>(newarg).setflag(status_flags::purely_indefinite);
			}
			return add(sumseq)+log(newarg);
		} else {
			if (!(options & expand_options::expand_function_args))
				ex_to<basic>(arg).setflag(status_flags::purely_indefinite);
		}
	}

	if (options & expand_options::expand_function_args)
		return log(arg.expand(options)).hold();
	else
		return log(arg).hold();
}

static ex atan_series(const ex &arg,
                      const relational &rel,
                      int order,
                      unsigned options)
{
	GINAC_ASSERT(is_a<symbol>(rel.lhs()));
	// method:
	// Taylor series where there is no pole or cut falls back to atan_deriv.
	// There are two branch cuts, one runnig from I up the imaginary axis and
	// one running from -I down the imaginary axis.  The points I and -I are
	// poles.
	// On the branch cuts and the poles series expand
	//     (log(1+I*x)-log(1-I*x))/(2*I)
	// instead.
	const ex arg_pt = arg.subs(rel, subs_options::no_pattern);
	if (!(I*arg_pt).info(info_flags::real))
		throw do_taylor();     // Re(x) != 0
	if ((I*arg_pt).info(info_flags::real) && abs(I*arg_pt)<_ex1)
		throw do_taylor();     // Re(x) == 0, but abs(x)<1
	// care for the poles, using the defining formula for atan()...
	if (arg_pt.is_equal(I) || arg_pt.is_equal(-I))
		return ((log(1+I*arg)-log(1-I*arg))/(2*I)).series(rel, order, options);
	if (!(options & series_options::suppress_branchcut)) {
		// method:
		// This is the branch cut: assemble the primitive series manually and
		// then add the corresponding complex step function.
		const symbol &s = ex_to<symbol>(rel.lhs());
		const ex &point = rel.rhs();
		const symbol foo;
		const ex replarg = series(atan(arg), s==foo, order).subs(foo==point, subs_options::no_pattern);
		ex Order0correction = replarg.op(0)+csgn(arg)*Pi*_ex_1_2;
		if ((I*arg_pt)<_ex0)
			Order0correction += log((I*arg_pt+_ex_1)/(I*arg_pt+_ex1))*I*_ex_1_2;
		else
			Order0correction += log((I*arg_pt+_ex1)/(I*arg_pt+_ex_1))*I*_ex1_2;
		epvector seq;
		if (order > 0) {
			seq.reserve(2);
			seq.push_back(expair(Order0correction, _ex0));
		}
		seq.push_back(expair(Order(_ex1), order));
		return series(replarg - pseries(rel, std::move(seq)), rel, order);
	}
	throw do_taylor();
}

void compile_ex(const ex& expr, const symbol& sym1, const symbol& sym2, FUNCP_2P& fp, const std::string filename)
{
	symbol x("x"), y("y");
	ex expr_with_xy = expr.subs(lst(sym1==x, sym2==y));

	std::ofstream ofs;
	std::string unique_filename = filename;
	global_excompiler.create_src_file(unique_filename, ofs);

	ofs << "double compiled_ex(double x, double y)" << std::endl;
	ofs << "{" << std::endl;
	ofs << "double res = ";
	expr_with_xy.print(GiNaC::print_csrc_double(ofs));
	ofs << ";" << std::endl;
	ofs << "return(res); " << std::endl;
	ofs << "}" << std::endl;

	ofs.close();

	global_excompiler.compile_src_file(unique_filename, filename.empty());
	// This is not standard compliant! ... no conversion between
	// pointer-to-functions and pointer-to-objects ...
	fp = (FUNCP_2P) global_excompiler.link_so_file(unique_filename+".so", filename.empty());
}

static ex atanh_eval(const ex & x)
{
	if (x.info(info_flags::numeric)) {

		// atanh(0) -> 0
		if (x.is_zero())
			return _ex0;

		// atanh({+|-}1) -> throw
		if (x.is_equal(_ex1) || x.is_equal(_ex_1))
			throw (pole_error("atanh_eval(): logarithmic pole",0));

		// atanh(float) -> float
		if (!x.info(info_flags::crational))
			return atanh(ex_to<numeric>(x));

		// atanh() is odd
		if (x.info(info_flags::negative))
			return -atanh(-x);
	}
	
	return atanh(x).hold();
}