* Implementation of GiNaC's initially known functions. */
/*
- * GiNaC Copyright (C) 1999-2018 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2024 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
return Order(x).hold();
}
+static ex Order_power(const ex & x, const ex & e)
+{
+ // Order(x)^e -> Order(x^e) for positive integer e
+ if (is_exactly_a<numeric>(e) && e.info(info_flags::posint))
+ return Order(pow(x, e));
+ // NB: For negative exponents, the above could be wrong.
+ // This is because series() produces Order(x^n) to denote the order where
+ // it gave up. So, Order(x^n) can also be an x^(n+1) term if the x^n term
+ // vanishes. In this situation, 1/Order(x^n) can also be a x^(-n-1) term.
+ // Transforming it to Order(x^-n) would miss that.
+
+ return power(Order(x), e).hold();
+}
+
static ex Order_expl_derivative(const ex & arg, const symbol & s)
{
return Order(arg.diff(s));
expl_derivative_func(Order_expl_derivative).
conjugate_func(Order_conjugate).
real_part_func(Order_real_part).
- imag_part_func(Order_imag_part));
+ imag_part_func(Order_imag_part).
+ power_func(Order_power));
//////////
// Solve linear system
//////////
-static void insert_symbols(exset &es, const ex &e)
-{
- if (is_a<symbol>(e)) {
- es.insert(e);
- } else {
- for (const ex &sube : e) {
- insert_symbols(es, sube);
+class symbolset {
+ exset s;
+ void insert_symbols(const ex &e)
+ {
+ if (is_a<symbol>(e)) {
+ s.insert(e);
+ } else {
+ for (const ex &sube : e) {
+ insert_symbols(sube);
+ }
}
}
-}
-
-static exset symbolset(const ex &e)
-{
- exset s;
- insert_symbols(s, e);
- return s;
-}
+public:
+ explicit symbolset(const ex &e)
+ {
+ insert_symbols(e);
+ }
+ bool has(const ex &e) const
+ {
+ return s.find(e) != s.end();
+ }
+};
ex lsolve(const ex &eqns, const ex &symbols, unsigned options)
{
for (size_t r=0; r<eqns.nops(); r++) {
const ex eq = eqns.op(r).op(0)-eqns.op(r).op(1); // lhs-rhs==0
- const exset syms = symbolset(eq);
+ const symbolset syms(eq);
ex linpart = eq;
for (size_t c=0; c<symbols.nops(); c++) {
- if (syms.count(symbols.op(c)) == 0) continue;
+ if (!syms.has(symbols.op(c)))
+ continue;
const ex co = eq.coeff(ex_to<symbol>(symbols.op(c)),1);
linpart -= co*symbols.op(c);
sys(r,c) = co;
}
// test if system is linear and fill vars matrix
- const exset sys_syms = symbolset(sys);
- const exset rhs_syms = symbolset(rhs);
+ const symbolset sys_syms(sys);
+ const symbolset rhs_syms(rhs);
for (size_t i=0; i<symbols.nops(); i++) {
vars(i,0) = symbols.op(i);
- if (sys_syms.count(symbols.op(i)) != 0)
+ if (sys_syms.has(symbols.op(i)))
throw(std::logic_error("lsolve: system is not linear"));
- if (rhs_syms.count(symbols.op(i)) != 0)
+ if (rhs_syms.has(symbols.op(i)))
throw(std::logic_error("lsolve: system is not linear"));
}