C++ ex类代码示例

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));
}